Following our last post about the FT Theory, now we will practice in Matlab with code exercises and solutions. The first thing you need to do is to download the code for this tutorial here: https://www.mathworks.com/matlabcentral/fileexchange/56443-discrete-fourier-transform-function-to-avoid-sampling-errors
EXERCISE 1: Calculate the FFT of a sinusoidal signal and analyse it
In this exercise, first, we will generate 64 samples of a sinusoidal signal (using the function sine) with frequency f=20 Hz and sampling frequency, fs=128 Hz. We will then calculate its DFT by suing the 64 points of the signal, we will represent its module and its phase. We are going to use the stem function to represent its module. For example, write the following code:
x=sine(start, L, f, fs);
title('DFT module representation (omegas-axis)')
title('DFT module representation (frequencies-axis)')
As this is a pure sinusoidal signal, we will obtain an unique significant value at the frequency of 20 Hz. We will named the obtained vector Y and it will represent the samples of the signal’s continuous spectrum. Observe that we will use the frequency interval 0≤Ѡ≤2π, which corresponds to the real frequencies 0≤f≤fs:
Figure 1. Sine DFT omega-axis and frequency-axis representation
As you can see, the sine DFT is a two deltas function. However, notice that in this exercise, the DFT has been calculated between 0≤2π. In next exercises, we will how to calculate between -π and π by using some Matlab functions.
Remember that the sine DFT is a complex number and what we have represented before it was its module. Now, the phase will look like this:
Figure 2. Sine DFT phase omega-axis and frequency-axis representation
EXERCISE 2: Leakage effect.
Now we are going to repeat the same exercise as before but we will use a signal frequency of 19 Hz. Why do you think this will cause a leakage effect? 🙂
If we represent the DFT with 64 points, being the sampling frequency fs=128 Hz, the leakage effect will happen because the sample frequency is not a integer multiple of the number of samples that we are taking (however, in the previous exercises, it was 128/64 = 2)
Therefore, we will observe that the DFT Matlab representation is not accurate: we are losing signal information in the frequency domain because the spectral components are spread around other frequencies. This error is the result of sampling the signal over a non-integer number of periods.
Let’s check this in the module of the DFT for the sine signal:
Figure 3. Sine DFT (leakage effect) omega-axis and frequency-axis representation
Therefore, instead of obtaining the two corresponding deltas at f=19 Hz and its replica, we have obtained an approximation of the spectrum in the frequencies around.
Regarding to the argument, we know that it is define as:
Also, we know that in Matlab, -π≤ATAN2(Y,X)≤π. Therefore, when representing the phase by using the function stem, we won’t observe the leakage effect very pronounced but we will still can see it because the interval between -π and π is not clearly defined:
Figure 4. Sine DFT phase (leakage effect) omega-axis and frequency-axis representation
Now, we suggest go a bit further and see what happens if you represent the DFT of 20 Hz sine signal but using a bigger number of points, a power of 2, such as 128 and 256. You should still observe the two deltas in -fc and fc but also some reflections in the nearest frequencies: this is one of the limitations of the Fast Fourier Transform algorithm; when taking more points for the signal spectrum, we will observe more spectral components due to the inaccuracy of this algorithm.
Now, if you try the 19 Hz sine’s DFT with 128, you will observe that there 2 deltas around -fc and fc, but as we saw before, you will observe the leakage effect as well. When using 256 points, you will still observe the leakage effect but, this time, the 2 deltas will be closer to -fc and fc.
EXERCISE 3: Representing the spectrum centered in 0 Hz
In the first exercise, we commented that the FFT Matlab algorithm, represents the spectrum in 0≤Ѡ≤2π, which corresponds to the real frequencies 0≤f≤fs. Also, remember that the Fourier transform is symmetric in the interval π≤Ѡ≤2π and this spectrum is equivalent to the one in the interval -π≤Ѡ≤0.
In the available code, you will see that we have created a DFT function that takes an input signal of period N and sampling frequency fs. The output is a vector that contains the module of its Fourier Transform in the interval -π≤Ѡ≤π or equivalently, -fs/2≤f≤fs/2. To do this, we have used the Matlab functions commented in our previous post: fft and fftshift. Also, in this function, we have restricted the signal period to be a multiple of 2 and greater than the sequence length. We have used the Matlab functions error and length to detect and throw errors in case this restrictions aren’t met.
We have taken into account the following considerations:
In Matlab, it is not possible to compute the continuous Fourier Transform, because the computer just works with a finite number of discrete or quantified values; therefore, the signal must be sampled and that’s why we use the Discrete Fourier Transform.
The previous point means that the Fourier Transform integrals are replaced by samples summations and the period, instead of being T (a real number) is N, and it has to be an integer number. In the case of periodic signals, which are the ones we study in this post, the DFT expression is as follows:
The Matlab DFT algorithm takes the following number of samples:
no. of operations=(no of samples)log2(no of samples)
Therefore, the number of samples needs to be a power of 2. If not, the number of samples will be segmented in summations of powers of 2 which will result in a less efficient algorithm.
In addition, with our function we will avoid the DFT downsampling because this will cause a lost of information. In order to get the original input signal of L samples, we need to do a DFT of, at least, L coefficients. In other words, if N is the fundamental period of the signal, the condition, L≤N needs to be met. If this requirement is not met, there will be overlapping and we will observe errors in the spectrum amplitudes, as we observed in previous posts.
Now, you can use that function and try different cases to understand why we apply these considerations. For example, start by generating a sine of 64 samples with a frequency of 20 Hz and a sampling frequency of 128 Hz. If you try to generate its DFT by using 65 points, you will get an error because the fundamental period needs to be a power of 2. Also, if you do N=4, you will get an error explaining that the fundamental period needs to be greater than the signal length.
On the other hand, if a sine of 64 samples with a frequencies of 30, 40 and 50 Hz and a sampling frequency of 128 Hz. If you try to generate its DFT using 64 points, you will be able to represent the spectrum (with less accuracy, as explained before) correctly.
EXERCISE 4: Downsampling and Aliasing
Using the function explained in the previous exercise, generate a sine of 70 Hz, sample it at 128 Hz and calculate its DFT. What do you observe?
In this case, we are downsampling the signal, so we will observe the aliasing effect. Notice that f=182Hz
Figure 5. Sine DFT (aliasing and leakage effects) frequency-axis representation
As you can observe, we have the same spectrum but for a sine of lower frequency (59 Hz).
EXERCISE 5: Inverse Fourier Transform (IFFT)
In this exercise, we are going to use the Matlab Function IFFT. We are going to compare the 20 Hz signal generated in the first exercise and, after calculating its DFT (using the function we have created above), we going to apply the IFFT (we will get x'[n]) and compare with the original sine, x[n]. Then, do x[n]-x'[n] and you will see that this is not zero, because of the error of the IFFT algorithm:
Figure 6. Sines IFFT comparison, n-axis representation
Now, do the same for a signal of 70 Hz, sampled at 128 Hz. As you already know, this will cause aliasing and the deltas won’t be at -fc and fc, but they will be at -59 and 59 Hz:
Figure 7. Sine DFT (aliasing effect), frequency-axis representation
Now, check the leakage effect by changing the sampling frequency to 1280 Hz. The DFT will be correct but inaccurate:
Figure 8. Sine DFT (leakage effect), frequency-axis representation
Finally, calculate the IFFT of the original signal and the aliased one and subtract them:
Figure 9. Sines IFFT comparison, n-axis representation
Again, the aliased signal reconstruction is not perfect and therefore, the subtraction is not zero.
In this tutorial, we have studied several concepts related to the DFT. We have used sinusoidal signals obtained by sampling a continuous signal in the time domain. By studying their spectra using Matlab, we have learnt that the DFT is calculated with the FFT function which represents this spectrum between 0 and 2π. Also, we saw that by using the function fftshift, we represent the spectrum between -π and π.
Remember that the conditions that need to be met to sample correctly are:
Matlab limits us to consider the number of DFT points taken, N, needs to be a power of 2.
In order to avoid the overlap of samples, the length of the sequence needs to be smaller or equal to N.
The sampling frequency needs to be at least, twice the bandwidth of the signal in order to avoid aliasing.
Finally, we also studied the leakage effect and see that this happens when the signal frequency is not an integer multiple of fs/N.
Therefore, if we don’t meet any of the previous conditions, there will be a lost of information in the signal.
We hope you enjoyed this post and it helps you to understand some concepts about the DFT and FFT. If you still have more questions, don’t hesitate in asking by leaving an comment below. Your request for future posts will be also very welcomed 🙂