FFT (fast fourier transform) is a method to find discrete fourier transform (DFT). To see the behaviour of a time domain signal in the frequency domain, we use FFT.
Here I find harmonic contents of a square wave of 50 Hz. The square wave amplitude is from -1 to +1.
#finding harmonics of the square wave clear clc N=4096; t=linspace(0,2*pi,N); # 2*pi*f = omega square wave equation of 50 hz f=square(2*pi*50*t); ##to plot the freq on x axis we decide the freq # interval is dt=t(2) i.e. diff between 2nd and 1st sample of "fun" dt=t(2); # nyquist's criterion applied n=[-N/2:1:N/2]; # last sample deleted to keep dimensions equal to those of "f" n(N+1)=; ## note these steps freq=2*pi*n./(N*dt); fun=dt*(fftshift(fft(f))); # signal strength matrix. Notice that ## the ABSOLUTE VALUE is considered strength=abs(fun); ## following code shows raw fft plot subplot(121) plot(freq,strength,"linewidth",2,"color","red") xlabel("freq","fontsize",18); ylabel("strength","fontsize",18); title("Raw FFT","fontsize",22) grid minor # find index of 0 value zeroindex=find(strength==0); # remove negative frequencies and the corrosponding strengths freq=freq(zeroindex:end)./(2*pi); strength=strength(zeroindex:end); ## The following code shows the fft output without -ve frequencies subplot(122) plot(freq,strength,"linewidth",2); title("Negative freq Removed","fontsize",22); xlabel("Actual frequency","fontsize",18); ylabel("strength","fontsize",18); grid minor print("sqrharm.png","-dpng")
While finding the harmonic contents of a square wave using “fft”, some important points were noted:
- GNU Octave has an inbuilt function “square” producing a square wave. To get a square wave of 50 Hz, the command is: f=square (2*pi*50*t)
- To study more harmonics, we need more samples. I considered 2^12=4096 samples. I could get harmonics up to 350 Hz (i.e. 7th harmonic of 50 Hz) using these samples.
- FFT gives data in mirror image form. So FFTSHIFT is also used.
- Outcome of fftshift(fft(function)) is a vector of complex numbers. So its absolute value is found out. In addition, multiplying the output (of fftshift and fft ) by “dt” scales the output. This enables us to know relative strengths of the frequency signals with respect to each other.
- In the above Octave script, I wanted to see only the right half of the original plot. This requires removing all negative frequencies. This is done by finding the index of zero of the vector and then manipulating the vector.
- To know the actual frequency, abscissa frequency must be divided by 2*pi .
- The plot on the right below shows the peaks from left to right at harmonic frequencies. They match the theoretical values derived from the fourier analysis of the square wave.