6.2 Using FFT to Obtain Simple Spectral Analysis Plots

How can you correctly scale the output of the FFT function to obtain a meaningful magnitude versus frequency plot?

Solution

Assume "x" is a vector containing your data. A sample vector used with this technical note is

Fs=1000;                              % sampling frequency
Fn=Fs/2;                              % Nyquist frequency
t=0:1/Fs:1;                           % time vector sampled at Fs Hz,
                                      % length of 1 second
x = sin(2*pi*t*200);                  % sine wave of 200 Hz.

First, you need to call the fft function. For the fastest possible ffts, you will want to pad your data with enough zeros to make its length a power of two. The fft function does this for you automatically if you give a second argument specifying the overall length of the fft, as demonstrated below:

NFFT=2.^(ceil(log(length(x))/log(2)));% Next highest power of 2
                                      % greater than or equal to
                                      % length(x)
FFTX=fft(x,NFFT);                     % Take fft, padding with zeros.
                                      % length(FFTX)==NFFT

If NFFT is even (which it will be, if the above two commands are used), then magnitude of the fft will be symmetric such that the first (1+NFFT/2) points are unique, and the rest are symmetrically redundant. FFTX(1) is the DC component of x, and fftx(1+NFFT/2) is the Nyquist frequency component of x.If NFFT is odd, however, the Nyquist frequency component is not evaluated, and the number of unique points is (NFFT+1)/2. This can be generalized for both cases to ceil((NFFT+1)/2).

NumUniquePts = ceil((NFFT+1)/2);
FFTX=FFTX(1:NumUniquePts);            % fft is symmetric, throw away
                                      % second half

Next, we calculate the magnitude of the fft:

MX=abs(FFTX);                         % Take magnitude of X
MX=MX*2;                              % Multiply by 2 to take into
                                      % account the fact that we
                                      % threw out second half of
                                      % FFTX above

Note that simply multiplying MX by 2 here is the quick-and-dirty solution. Technically, the DC component and the Nyquist component (if it exists, see above) are unique and should not be multiplied by 2. To be accurate, the following lines can be added:

MX(1)=MX(1)/2;
if ~rem(NFFT,2),
     MX(length(MX))=MX(length(MX))/2;
end

Next, we need to take into account the fact that MATLAB doesn't scale the output of fft by the length of the input:

MX=MX/length(x);                      % Scale the FFT so that it is
                                      % not a function of the length
                                      % of x

Now we create a frequency vector:

f=(0:NumUniquePts-1)*2/NFFT;          % This is an evenly spaced
                                      % frequency vector with
                                      % NumUniquePts points.
f=f*Fn;                               % Multiply this by the Nyquist 
                                      % frequency (Fn ==1/2 sample
                                      % freq.)

The above two lines can be combined to

f=(0:NumUniquePts-1)*2*Fn/NFFT;
Finally, we can generate the plot
plot(f,MX);

Bringing this all together, we get the following M-file:

Fs=1000;                              % sampling frequency
Fn=Fs/2;                              % Nyquist frequency
t=0:1/Fs:1;                           % time vector sampled at Fs Hz,
                                      % length of 1 second
x = sin(2*pi*t*200);                  % sine wave of 200 Hz.
NFFT=2.^(ceil(log(length(x))/log(2)));% Next highest power of 2
                                      % greater than length(x).
FFTX=fft(x,NFFT);                     % Take FFT, padding with zeros.
                                      % length(FFTX)==NFFT
NumUniquePts = ceil((NFFT+1)/2);
FFTX=FFTX(1:NumUniquePts);            % FFT is symmetric, throw away
                                      % second half
MX=abs(FFTX);                         % Take magnitude of X
MX=MX*2;                              % Multiply by 2 to take into
                                      % account the fact that we
                                      % threw out second half of 
                                      %FFTX above
MX(1)=MX(1)/2;                        % Account for endpoint
                                      % uniqueness
MX(length(MX))=MX(length(MX))/2;      % We know NFFT is even
MX=MX/length(x);                      % Scale the FFT so that it is
                                      % not a function of the length
                                      % of x.
f=(0:NumUniquePts-1)*2*Fn/NFFT;
plot(f,MX);