原创 Lesson 3 : MATLAB - Introductory FFT Tutorial

http://blinkdagger.com/matlab/matlab-introductory-fft-tutorial

 

Introduction

In this tutorial, we will discuss how to use the fft (Fast Fourier Transform) command within MATLAB. The fft command is in itself pretty simple, but takes a little bit of getting used to in order to be used effectively.

When we represent a signal within matlab, we usually use two vectors, one for the x data, and one for the y data. The fft command only operates on the y-data (converting the y-data from the time domain into the frequency domain), so it’s up to the user to determine what the x-data in the frequency domain will be! This tutorial will show you how to define your x-axis so that your fft results are meaningful. In addition, it will show you how to obtain a two-sided spectrum as well as a positive frequency spectrum for a given signal.

A Simple Example

1. Let’s start off with a simple cosine wave, written in the following manner:
   

   

2. Next, let’s generate this curve within matlab using the following commands:
   
   
   fo = 4; %frequency of the sine wave
   Fs = 100; %sampling rate
   Ts = 1/Fs; %sampling time interval
   t = 0:Ts:1-Ts; %sampling period
   n = length(t); %number of samples
   y = 2*sin(2*pi*fo*t); %the sine curve
   
   %plot the cosine curve in the time domain
   sinePlot = figure; plot(t,y)
   xlabel('time (seconds)')
   ylabel('y(t)')
   title('Sample Sine Wave')
   grid

Here’s what we get:

1. When we take the fft of this curve, we would ideally expect to get the following spectrum in the frequency domain (based on fourier theory, we expect to see one peak of amplitude 1 at -4 Hz, and another peak of amplitude 1 at +4 Hz):

   

2. There is also a phase component, but we’ll discuss that in a future tutorial.

Using Matlab’s FFT Command

So now that we know what to expect, let’s use MATLAB’s built in fft command to try to recreate the frequency spectrum:

%plot the frequency spectrum using the MATLAB fft command
matlabFFT = figure; %create a new figure
YfreqDomain = fft(y); %take the fft of our sin wave, y(t)
stem(abs(YfreqDomain)); % use abs command to get the magnitude
% similary, we would use angle command to get the phase plot!
% we'll discuss phase in another post though!

xlabel('Sample Number')
ylabel('Amplitude')
title('Using the Matlab fft command')
grid

axis([0,100,0,120])

This doesn’t quite look like what we predicted above. If you notice, there are a couple of things that are missing.

• The x-axis gives us no information on the frequency. How can we tell that the peaks are in the right place?
• The amplitude is all the way up to 100
• The spectrum is not centered around zero

A Custom Function for fft to Obtain a Two-Sided Spectrum

Here is a helpful function I learned at Harvey Mudd College which will simplify the process of plotting a two-sided spectrum. Copy this code into an m-file and save it.

function [X,freq]=centeredFFT(x,Fs) % this is a custom function that helps in plotting the two-sided spectrum
% x is the signal that is to be transformed
% Fs is the sampling rate
N=length(x); %this part of the code generates that frequency axis

if mod(N,2)==0 
    k=-N/2:N/2-1; % N even
else
    k=-(N-1)/2N-1)/2; % N odd
end
T=N/Fs;
freq=k/T; %the frequency axis

%takes the fft of the signal, and adjusts the amplitude accordingly
X=fft(x)/N; % normalize the data
X=fftshift(X); %shifts the fft data so that it is centered

This is a relatively simple function to use. The function outputs the correct frequency range and the transformed signal. It takes in as input the signal to be transformed, and the sampling rate.

Let’s use the sine wave from above and do a quick example (Remember to set the Matlab directory to the location where you saved the previous m-file). Now, copy and paste these commands into the Matlab command prompt.

[YfreqDomain,frequencyRange] = centeredFFT(y,Fs);
centeredFFT = figure; %remember to take the abs of YfreqDomain to get the magnitude!
stem(frequencyRange,abs(YfreqDomain));
xlabel('Freq (Hz)')
ylabel('Amplitude')
title('Using the centeredFFT function')
grid
axis([-6,6,0,1.5])

Here’s what you should see:

As you can see, this plot is basically identical to what we would expect! We get peaks at both -4 Hz and +4 Hz, and the amplitude of the peaks are 1.

Redundant Information in the FFT

As you can see from the plots above, the information within the frequency spectrum is entirely symmetric. Thus, we only need one side of the spectrum. In general, the positive side of the spectrum is used, while the negative side is ignored. So let’s adjust our function above so that we only get the positive frequencies.

A Custom Function for fft to Obtain only the Positive Frequencies

The following function is a modification of the above function, and will help you plot only the positive frequencies of the spectrum.

function [X,freq]=positiveFFT(x,Fs) N=length(x); %get the number of points
k=0:N-1; %create a vector from 0 to N-1
T=N/Fs; %get the frequency interval
freq=k/T; %create the frequency range
X=fft(x)/N; % normalize the data

%only want the first half of the FFT, since it is redundant
cutOff = ceil(N/2);

%take only the first half of the spectrum
X = X(1:cutOff);
freq = freq(1:cutOff);

 

Once again, let’s use the same sine wave and put it through this function. Copy and paste the following code into the Matlab command prompt.

[YfreqDomain,frequencyRange] = positiveFFT(y,Fs);
positiveFFT = figure;
stem(frequencyRange,abs(YfreqDomain));
set(positiveFFT,'Position',[500,500,500,300])
xlabel('Freq (Hz)')
ylabel('Amplitude')
title('Using the positiveFFT function')
grid
axis([0,20,0,1.5])

 

Here’s what you should get:

These two functions are very useful, and I still use them all the time!

Power of 2

The fft command within Matlab allows you to specify how many data points are in the transform. The Matlab documentation recommends that a power of 2 be used for optimal computation time. In my experience, there really isn’t a need to specify N as a power of 2. By using the next greatest power of 2, the fft command pads the original signal with zeros and proceeds to do a FFT on the signal. I’ve done some quick runs using fft with N as a power of 2 and N not as a power of 2, and the speed difference was neglible. In some cases, a 120 point FFT took LESS time than a 128 point FFT in some of my runs.

I don’t know exactly how the Matlab fft works, but I believe that it internally pads the signal with zeroes to the next greatest power of 2, performs the fft , then spits out an answer without the padded zeros.

I highly encourage anyone with greater knowledge to shed some light on this topic!

Recap and Future Topics

In this post, we talked primarily about how to use the fft command to create a frequency spectrum from a given signal. Two important things we did were to appropriately define the frequency axis, adjust the amplitude, and to view the spectrum as a two-sided, or one-sided spectrum.

In this post, we used a very simple example signal that was very well behaved. Unfortunately, this will rarely be the case. In the next upcoming posts, we will discuss zero-padding, windowing, filtering, and other techniques that will help you interpret and analyze the frequency spectrum of various signals.

Another thing that was not covered in this post, but is of imperative importance, is the phase. Stay tuned, and be ready for more material in the future!

Download source files and Further Reading

You can download the source files here.

Brush up your Fourier by reading about the theory and background at these links:

http://www.complextoreal.com/chapters/fft1.pdf
http://www.dspguide.com/ch8.htm

This is the end of the tutorial.