Hi all!
I hope this is the right subcategory to ask in. I’m playing around with the FFTW package to perform some Fourier transforms. I got how to define a signal and transform it, but I can’t get how (and why) can I shift the frequency domain to reflect the actual frequencies of my signal. My Google-fu made me understand I need to call fftshift() and fftfreq(), but apparently I can’t understand how and why.
As an example, I have this script:
using FFTW
using Plots
x = 0:100
y = sin.(x)
# y = [abs(x) <= 10 ? 1 : 0 for x in x]
transform = fftshift(fft(y))
f = fftshift(fftfreq(length(x)) * length(x))
plot(f, abs.(transform))
However, I got this:
Inspecting it with plotlyjs() I can see the spike on the right is at x ~= 16. Why, since the original signal had frequency = 1?
Is this the right way to go? If not, what it is? And why? Thanks 
1 Like
The original signal repeats how many times per 100 samples? It’s 100/(2pi) ~16
Time is measured by the duration of the sample… Frequency is measured by cycles per duration of sample
If you want to calculate with physical frequencies you will need to explicitly tell fftfreq the sample rate. Like if you are sampling every millisecond then sample rate is 1/0.001
4 Likes
If the original angular frequency is 1, to get it we have to multiply fftfreq() by 2π. See annotated code below:
using FFTW, Measures, Plots; gr()
x = 0:100
N = length(x)
y = sin.(x) # ω0*x = x; with angular frequency ω0 = 1
Y = fft(y) # Fundamental DFT period T = N*ΔT = N; with ΔT = 1
ω = 2π*fftfreq(N, 1) # 2π-scaler to get angular frequency. Sampling rate fs = 1/ΔT = 1
p1 = plot(x, y, st=:line, m=:o, ms=2, xlabel="Sample number");
p2 = plot(ω, abs.(Y), st=:stem, m=:o, ms=2, xlabel="Radians per sample");
plot(p1, p2, margins=5mm, size=(800,350))
2 Likes
Uhm, I think I see the problem but I can’t understand the solution: what frequencies does the DFT outputs? Radians? If this is so, shouldn’t I get the “right” frequencies in the output (and graph) without calling fftfreq (since sin.(x) also expects radians)?
Here is a recent thread with similar discussion.
Cycles per window length.
1 Like
Give the sample freqiency as the second argument to fftfreq, and do not multiply by length(x). I usually do something like this to get it right:
using FFTW
using UnicodePlots
let
fs = 1000 # sample freq
duration = 1
timestamps = range(0, duration, step=1 / fs) # step = sample period = 1 / fs
f = 250
A = 2
s = @. A * sin(f * 2π * timestamps)
plt1 = scatterplot(timestamps, s)
transform = fft(s)|>fftshift
freqs = fftfreq(length(s), fs)|>fftshift
plt2 = scatterplot(freqs, transform .|> abs)
display(plt1)
display(plt2)
end
┌────────────────────────────────────────┐
2 │⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶⠶│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
-2 │⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀│
└────────────────────────────────────────┘
⠀0⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀1⠀
┌────────────────────────────────────────┐
1000 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⢈⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⢠⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀│
0 │⣀⣀⣀⣀⣀⣀⣀⣀⣀⡼⢇⣀⣀⣀⣀⣀⣀⣀⣀⣀⣇⣀⣀⣀⣀⣀⣀⣀⣀⡸⢧⣀⣀⣀⣀⣀⣀⣀⣀⣀│
└────────────────────────────────────────┘
⠀-500⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀500⠀
EDIT :
f=250 makes the time domain look funny, lower frequencies make the dot harder to spot in the frequency domain, and may reault in leakage. I would reccomend copying the code and varying e.g. the signal frequency yourself. You can also try adding another signal, to get two distrinct frequency peaks.
The DFT assumes the input data is periodic and it outputs amplitudes of the sinusoids whose frequencies are defined as a multiple of a fundamental frequency that is equal to the inverse of the input data duration (called also fundamental period T = N*dT; assume dT=1 in your example). The first sample is at 0, the second sample corresponds to 1 cycle per fundamental period (i.e., think of one sine wave with a period equal to the input window length), the third sample corresponds to 2 cycles per fundamental period, etc. (after N/2 the negative frequencies wrap around). If we want to convert cycles per fundamental period to radians per fundamental period (i.e. angular frequencies instead of regular frequencies), we need to multiply by 2pi.
1 Like
I think I’m slowly getting now. I think I misunderstood some quirks of the DFT, and the problems that arise in the sampling process. For instance, I now get that the DFT has no way of teling what the frequencies were, if you don’t take into account the sampling frequency. That iiuc it’s the reason why fftfreq() asks for it!
I am trying to fill my gaps using this tutorial book.
Anyway, a comment on the line
fftshift(fftfreq(length(x)) * length(x))
I took it from another tutorial, but iiuc it’s not idiomatic: there’s no point in multiplying for lenght(x) (the number of samples), the usual way to go is to pass the sampling time to fftfreq(), right? With this, frequencies are already good, I think?
1 Like
There is no need to define physical units in the mathematical definition of FFT, as hinted in my previous post. However, the main applications of FFT are in engineering where physical units of time or space are commonly used.
Examples have been posted multiple times in Discourse, including by TheLateKronos above. FWIW find here below another time domain example, annotated.
using FFTW, Measures, Plots; gr()
Δt = 2e-3 # sampling period [s]
fs = 1/Δt # sampling frequency [Hz]
t = 0:Δt:0.24 # times sampled in [s]
N = length(t)
f0 = 25.0 # sinusoid frequency [Hz]
yt = sin.(2π*f0*t) # sample values at times t
Yf = fft(yt)
f = fftfreq(N, fs) # frequencies of the Yf amplitudes
p1 = plot(t, yt, m=:o, ms=2, xlabel="Time [s]")
p2 = plot(f, abs.(Yf), st=:stem, lw=0.4, m=:o, ms=2, xlabel="Frequency [Hz]")
plot(p1, p2, size=(800,350), margins=5mm)
1 Like
Many thanks! I think I’m getting it now. 
One more thing: in Julia, which is is the coefficient of the DFT? The unitary one (1/\sqrt(N))? And why? (I can open a new topic if this is too OT)
