Implement FFT with Matlab and What is digital filter

Introduction

Let us understand FFT. It is Fast Fourier Transform, an algorithm to calculate DFT or discrete fourier transform in fast and efficient way. The first question that arises seeing the title is what the hell a tutorial on FFT doing in the new article section of code project in the year 2012 when the algorithm is about 50 years old. Great Question. By the way, you can refer following links for understanding the same in a complicated way.

http://en.wikipedia.org/wiki/Discrete_Fourier_transform

http://en.wikipedia.org/wiki/Fast_Fourier_transform

Background

One of the reasons, I thought I must write this tutorial is because I have seen many engineers finding it difficult to understand what exactly is a frequency domain, what exactly is filtering. How can we filter a signal in simple way. What is this symmetrical features of FFT and more importantly understanding the core facts about FFT. So if you are a programming Pro and a killer programmer having written your own DSP logics and algorithms, then this article may not give you any new or interesting things. But if you have just picked signal processing, you might find it quite helpful (Hopefully).

Using the Code

First Let us construct a simple sinusoidal signal of 50Hz with amplitude=5. I am considering a sampling frequency of 100 times the message frequency(it should it least be twice as par Nequist rate), which means I will collect 1000 samples from the actual analog sinusoidal signal. So if the frequency is 50, that means you will get one sample at 1/50 or .02 seconds. If you want 10 samples, you have to gather the samples for .2 seconds. You will sample the signal at every 1/5000 seconds.

Code MATLAB :

 

f=50;

A=5;

Fs=f*100;

Ts=1/Fs;

t=0:Ts:10/f;

x=A*sin(2*pi*f*t);

plot(x),grid on

Figure 1: Sinusoidal Signal of frequency 50Hz. Plot of Variable x.So the figure shows 10 cycles of sine wave of 50 Hz. You can change the plot to stem type to see the samples.

So lets take FFT.

Code MATLAB :

 

F=fft(x);

You must have studied about N point FFT. If you do not specify N, then by default N is length of message signal. In this case it is 1001.

Very important thing: FFT divides your Sampling frequency into N equal parts and returns the strength of the signal at each of these frequency levels. What it means is you are dividing frequencies from 0 to 5000 into 1001 equal parts. So first point in fft is 5Hz, next represents 10 Hz and so on. As the actual frequency of your message is 50 Hz, you must get highest peak at that level. So let us plot FFT.

Code MATLAB :

 

plot(real(F)),grid on

Figure 2: Plot of real part of the FFT output of signal x with highest peak marked. Plot of real part of Variable F.What you see? The highest value is located at 11th position, which is nothing but 10th sample as matlab starts from 1.

10*5?=50Hz

It is so simple. But wait why do we get the imaginary part also? There is absolutely nothing called Imaginary values, they reflect past values.

Your signal is Sinusoidal, so you get negative cycle also right. Check this.

Code MATLAB :

plot(real(imag(F))),grid on

Figure 3: Plot of imaginary part of the FFT output of signal x with highest peak marked. Plot of imaginary part Variable F.

what you see? 

Yes the same thing with one negative amplitude peak also. I hope you get the idea.

Now let us understand how you can have some filtering with FFT. What I will do is mix couple of more signal with different frequency with same amplitude and same number of samples with your signal x.

Code MATLAB :

x1=A*sin(2*pi*(f+50)*t);

x2=A*sin(2*pi*(f+250)*t);

x=x+x1+x2; 

 

Figure 4: Plot of hybrid signal x containing 50Hz,100Hz,300Hz frequency components. Plot of Variable x.

Our signal now contains three components, a 50Hz,100Hz and 300Hz.

Its fft plot is the proof. Three frequencies at three levels, exactly at 50,100 and 300Hz.

Figure 5: Plot of real part of the FFT output of Hybrid signal x.

I want to recover my original signal. So if I keep the sample at 11th sample and leave the rest and then take ifft, I must get back my signal right. But the values of both 10th and 11th sample must be kept. It is called windowing.Figure 6: Desired frequency band marked as a rectangular window in FFT plot of hybrid signal x .See the red color window I have put over the FFT response of our Hybrid signal.

Code MATLAB :

F2=zeros(length(F),1);

F2(10:11)=F(10:11);

xr=ifft(F2);

plot(real(xr)),grid on

Figure 7: Recovery of Original signal of Frequency component 50Hz from hybrid signal shown in figure 4, using rectangular window.

What you see? It is your original signal. But with less amplitude. This is because you choose a rectangular window and the loss is due to negative filter gain. You can now read a bit of digital filters and apply that.

Here is the complete code which may not be needed, but still I am posting it, in case you decide to spend some time and play with it.

Code MATLAB :

clear all;

clc;

close all

f=50;

A=5;

Fs=f*100;

Ts=1/Fs;

t=0:Ts:10/f;

x=A*sin(2*pi*f*t);

x1=A*sin(2*pi*(f+50)*t);

x2=A*sin(2*pi*(f+250)*t);

x=x+x1+x2;

plot(x)

F=fft(x);

figure

N=Fs/length(F);

baxis=(1:N:N*(length(x)/2-1));

plot(baxis,real(F(1:length(F)/2)))