DSP basics: Fast Fourier Transform (FFT)

The Fast Fourier Transform (FFT) is a fast discrete Fourier Transform (DFT) algorithm that is widely used in Digital Signal Processing (DSP) for analysis and processing in the frequency domain. In DSP, it is often necessary to translate signals from the time domain to the frequency domain (and back) in order to:

• Analyse the spectrum of a signal
• Filter signals
• Compress data
• Implement modulation algorithms (OFDM in Wi-Fi, 5G)
• Process radar and medical signals

Let's see how to apply the FFT function from the FFTW.jl plug-in library to the task of signal spectrum analysis:

using FFTW

A simple signal and its spectrum

Let's get acquainted with the peculiarities of using the function fft to analyse the spectrum on the example of a real sinusoidal signal with duration of 0.2 sec, sampling rate of 2 kHz, amplitude of 0.5 and repetition rate 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 samples. This is important for further analyses.

nsamp = length(t)

400

Let's apply the function fft to it and visualise 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 the samples of the FFT output modulo:

plot(abs.(A))

We observe two peaks corresponding to the positive and negative frequency region. But this representation does not give us reliable information about either the frequency or the amplitude of these peaks.

For a more explicit representation of the spectrum, we will use the function fftshift to shift the output vector of the FFT, and also set a vector of frequencies, which we will postpone along the X axis. We observe the spectrum of the signal in the so-called first Nyquist zone, which extends from -fs/2 to +fs/2. We can find the RBW (resolution bandwidth), or frequency grid step, by dividing the whole 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, " Гц")
freq_vec = -fs/2:fs/nsamp:fs/2-df;
plot(freq_vec, abs.(B))

Разрешение по частоте (RBW) = 5.0 Гц

Now we see two peaks at -100 Hz and +100 Hz. But their amplitudes do not match the amplitude of the sinusoid under consideration.

To accurately estimate the energy of the signal on its spectrum, we normalise the FFT output vector by the number of samples, i.e. divide by nsamp. But since half of the energy of our actual signal has "gone" to the region of negative frequencies, we need to multiply the resulting spectrum by 2. And for displaying we will use a line graph (stem) with round markers in the area 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 = "Амплитудный спектр синусоиды")

Разрешение по частоте (RBW) = 5.0 Гц

We observe a peak at 100 Hz with an amplitude of 0.5, which corresponds exactly to the parameters of the original sinusoidal signal.

Change in frequency resolution

Now let's try to change the number of samples of the input signal ("window" of the FFT). Let's input the first 80 samples of the sinusoid to the function fft. In this case, the output of the FFT algorithm will also be 80 complex samples, which means that the pitch of the frequency grid will change. The signal has become 5 times shorter, so 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))

Разрешение по частоте (RBW) = 25.0 Гц

The change in frequency grid did not affect the estimation of the amplitude and frequency of the signal. The 25 Hz discretisation allows us to "exactly hit" the spectral component at 100 Hz. But what happens if the spectrum step is not a multiple of the frequency of the estimated sinusoid?

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))

Разрешение по частоте (RBW) = 18.181818181818183 Гц

We observe two peaks at 91 Hz and 109 Hz. This is a spectral leakage phenomenon: the energy of the signal at 100 Hz has been captured by the main and side lobes (see DSP theory) of the window function on neighbouring samples of the frequency grid.

Addition of zeros at FFT

Zero padding in FFT is a technique that involves padding the original signal with zeros to the desired length before performing the FFT. Adding zeros increases the number of FFT points, reducing the pitch of the frequency grid and making the spectrum visually smoother. But it does not add new information about the signal - it simply interpolates the existing spectrum.

Let's augment our sine wave with zeros up to a vector length of 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))

Добавлено нулей = 1600
Разрешение по частоте (RBW) = 1.0 Гц

We can notice that the peak of the spectrum is at 0.1 instead of 0.5. This is due to the fact that the signal has become five times longer due to the zeros, and its energy, accordingly, has also been attenuated by a factor of five.

The classical algorithm works more efficiently when the number of FFT points is degree 2, and the addition of zeros is useful for length equalisation and visualisation, but it does not increase the true frequency resolution. To really improve spectral analysis, you need to increase the duration of the signal recording.

Interactive spectrum analysis

Finally, let us apply the spectrum estimation techniques to 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 us 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)

Let us now use the scripting functionality of Engee to create the code cell mask. We can change the value of frequency resolution (RBW) using the tool slider.

The code cell calculates the number of FFT output points, compares it with the number of points of the discrete signal 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 = "Амплитудный спектр сигнала")

Точек БПФ:  200    Отсчётов сигнала:  400
Усечение входного окна

Try experimenting with different values of RBW, from 1 to 20 Hz, and see the accuracy of frequency and amplitude estimation of individual spectral components of the signal.

Conclusion

We got acquainted with the peculiarities of using the function fft of the library FFTW.jl for estimating the spectrum of discrete signals, and also learnt how to apply the maxima of code cells and perform the tasks of digital signal processing in the interactive mode.