Using FFT

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:

  1. 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)
  2. 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.
  3. FFT gives data in mirror image form. So FFTSHIFT is also used.
  4. 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.
  5. 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.
  6. To know the actual frequency, abscissa frequency must be divided by  2*pi .
  7. 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.

sqrharm