DSP Basics: Fast Fourier Transform (FFT)
Fast Fourier Transform (FFT, FFT — Fast Fourier Transform) is an algorithm for fast computing of the discrete Fourier transform (DFT), which is widely used in digital signal processing (DSP) for analysis and processing in the frequency domain. In DSP, it is often necessary to transfer signals from the time domain to the frequency domain (and vice versa) in order to:
- Analyze the signal spectrum
- Filter the signals
- Compress data
- Implement modulation algorithms (OFDM in Wi-Fi, 5G)
- Process radar and medical signals
Let's look at how to use the FFT function from the plug-in library. FFTW.jl for the task of signal spectrum analysis:
using FFTW
A simple signal and its spectrum
Let's get acquainted with the features of using the function fft to analyze the spectrum using the example of a real sinusoidal signal with a duration of 0.2 seconds, with a sampling frequency of 2 kHz, an amplitude of 0.5 and a repetition frequency of 100 Hz:
fs = 2000;
dt = 1/fs;
t = 0:dt:0.2-dt;
sinewave = 0.5 * sin.(2*pi*t*100);
plot(t,sinewave)
It is worth noting the length of our discrete signal sinewave - 400 counts. This is important for further analysis.
nsamp = length(t)
Let's apply the function to it fft and visualize the result:
A = fft(sinewave);
plot(A)
The output of the FFT function is a complex vector A with the same number of samples as the original signal in the time domain.
To estimate the amplitude spectrum of the signal, it is necessary to take samples of the FFT output modulo:
plot(abs.(A))
We observe two peaks corresponding to the positive and negative frequency ranges. But this representation does not give us reliable information about either the frequency or the amplitude of these peaks.
For a more explicit display of the spectrum, we use the function fftshift to shift the output vector of the FFT, and also set the frequency vector, which is shifted along the X-axis. We observe the signal spectrum in the so-called first Nyquist zone, which extends from -fs/2 to +fs/2. We can find out the frequency resolution (RBW, resolution bandwidth), or frequency grid pitch, by dividing the entire frequency range under consideration (either from -fs/2 to +fs/2, or from 0 to fs) by the number of FFT points.:
B = fftshift(A);
println("Разрешение по частоте (RBW) = ", fs/nsamp, " Гц")
df = fs/nsamp;
freq_vec = -fs/2:df:fs/2-df;
plot(freq_vec, abs.(B))
Now we see two peaks at frequencies of -100 Hz and +100 Hz. But their amplitudes do not correspond to the amplitude of the sinusoid in question.
To accurately estimate the energy of a signal on its spectrum, we normalize the output vector of the FFT by the number of samples, that is, divide by nsamp. But since half of the energy of our actual signal has "gone" into the negative frequency range, we need to multiply the resulting spectrum by 2. Well, we'll use a bar graph to display (stem) with round markers in the range of positive frequencies from 0 to 400 Hz:
S = 2*B/nsamp;
println("Разрешение по частоте (RBW) = ", fs/nsamp, " Гц")
plot(freq_vec, abs.(S),
line=:stem,
marker=:circle,
xlim = (0, 400),
title = "Амплитудный спектр синусоиды")
We observe a peak at a frequency of 100 Hz with an amplitude of 0.5, which exactly corresponds to the parameters of the original sinusoidal signal.
Changing the frequency resolution
Now let's try to change the number of samples of the input signal (the "window" of the FFT). Let's apply the functions to the input fft the first 80 samples of the sine wave. In this case, the output of the FFT algorithm will also have 80 complex samples, which means that the step of the frequency grid will change. The signal has become 5 times shorter, therefore, the frequency resolution has become 25 Hz.:
winsize = 80;
rbw80 = fs/winsize;
freq80 = 0:rbw80:fs - rbw80;
S80 = (2*fft(sinewave[1:winsize]))./winsize;
println("Разрешение по частоте (RBW) = ", rbw80, " Гц")
plot(freq80, abs.(S80),
line=:stem,
marker=:circle,
xlim = (0, 400))
Changing the frequency grid did not affect the estimation of the amplitude and frequency of the signal. A 25 Hz discount allows us to "accurately hit" the spectral component at 100 Hz. But what happens if the pitch of the spectrum is not a multiple of the frequency of the estimated sine wave?
winsize = 110;
rbw110 = fs/winsize;
freq110 = 0:rbw110:fs - rbw110;
S110 = (2*fft(sinewave[1:winsize]))./winsize;
println("Разрешение по частоте (RBW) = ", rbw110, " Гц")
plot(freq110, abs.(S110),
line=:stem,
marker=:circle,
xlim = (0, 400))
We observe two peaks at frequencies of 91 Hz and 109 Hz. This is the phenomenon of spectral leakage: the energy of the signal at 100 Hz was captured by the main and side lobes (see DSP theory) of the window function on neighboring samples of the frequency grid.
Addition of zeros in FFT
Zero padding in FFT is a technique that involves adding zeros to the original signal to the desired length before performing the FFT. Adding zeros increases the number of FFT points, reducing the frequency grid pitch and making the spectrum visually smoother. But this does not add new information about the signal — it just interpolates the existing spectrum.
Let's add zeros to our sinusoid to the length of the vector in 2000 samples:
winsize = 2000;
rbw2000 = fs/winsize;
freq2000 = 0:rbw2000:fs - rbw2000;
numofzeros = winsize - length(sinewave);
println("Добавлено нулей = ", numofzeros)
S2000 = (2*fft([sinewave; zeros(numofzeros)]))./winsize;
println("Разрешение по частоте (RBW) = ", rbw2000, " Гц")
plot(freq2000, abs.(S2000),
line=:stem,
marker=:circle,
xlim = (0, 400))
You can see that the peak of the spectrum is at 0.1, not 0.5. This is due to the fact that the signal has become five times longer due to the zeros, and its energy, respectively, has also been weakened five times.
The classical algorithm works more efficiently when the number of FFT points is a power of 2, and adding zeros is useful for length alignment and visualization, but it does not increase the true frequency resolution. For a real improvement in spectral analysis, it is necessary to increase the recording time of the signal.
Interactive spectrum analysis
In conclusion, we will apply spectrum estimation techniques for a signal that is the sum of three sinusoids with frequencies of 80 Hz, 120 Hz, and 245 Hz, and amplitudes of 10, 5, and 2.4, respectively. Let's display the shape of the resulting vector in the time domain:
three_sines = [10 5 2.4] .* sin.(2*pi*t*[80 120 245]);
sig = sum(three_sines, dims = 2);
sig = vec(sig);
plot(t,sig)
And now let's use the functionality of the scripts. Engee to create a code cell mask. We can change the frequency resolution value (RBW) using the tool slider.
In the cell code, the number of FFT output points is calculated, compared with the number of discrete signal points in time, and if one exceeds the other, either the addition of zeros or the truncation of the input "window" of the algorithm is performed. Then the amplitude spectrum is calculated and drawn.
RBW = 10 # @param {type:"slider",min:1,max:20,step:1}
nfft = Int64(round(fs/RBW));
println("Точек БПФ: ", nfft,
" Отсчётов сигнала: ", nsamp);
if nsamp >= nfft
win = sig[1:nfft];
Y = 2*fft(win)./nfft;
println("Усечение входного окна")
else
win = [sig; zeros(nfft - nsamp)];
Y = 2*fft(win)./nsamp;
println("Дополнение нулями: интерполяция спектра")
end
ff = (0:nfft-1)*fs/nfft;
plot(ff, abs.(Y),
line=:stem,
xlim = (0, 400),
marker =:circle,
lab = "Сумма синусоид",
title = "Амплитудный спектр сигнала")
Try experimenting with different values. RBW, from 1 to 20 Hz, and look at the accuracy of estimating the frequency and amplitude of individual spectral components of the signal.
Conclusion
We got acquainted with the features of using the function fft libraries FFTW.jl to evaluate the spectrum of discrete signals, and also learned how to apply the maxima of code cells and perform digital signal processing tasks interactively.