Making matlab's fft() function useful

I use matlab almost exclusively for data processing. Matlab has many built-in functions that really make signal processing very quick and efficient. However, their FFT function is not one of them. Matlab's FFT is what I tend to think of as a ``mathematician's FFT''. It computes the basic transform exactly as written in the math textbooks. But what we are after is more of an ``engineer's FFT''. The main shortfalls of matlab's FFT function are:

• The output is not scaled to engineering units. I.e. if my input vector has units of inches (or mm, etc), I expect that the output is also in inches. But it's not.
• The output is in real-imaginary components, whereas magnitude-phase is usually more useful.
• The output includes components at negative frequencies. Negative frequencies are an interesting mathematical concept, but rarely useful for practical vibration analysis.
• There is no frequency vector in the output. i.e. all you get is y=fft(x). We'd like to get [y, freq]=fft(x), so that we could immediately do plot(freq,y), without having to compute the frequency vector manually every time.
• Even worse, the default frequency vector is not what most people would expect. If you're recognized that there are negative frequencies in the output, most people expect that if they do y=fft(x) for an x of length N, that y(1) corresponds to the frequency -Fs/2, y(N/2) corresponds to 0 Hz, and y(N) correspondings to Fs/2. But that's not the way matlab does it. Actually y(1) corresponds to 0, y((N-1)/2) corresponds to Fs/2, y(N/2) corresponds to -Fs/2, and y(N) to Fs/N. To get the order you were expecting, you need to do y=fftshift(fft(x)).

This shortfalls are illustrated in Figure 12. The input is just the sum of two sines waves (1 Hz and 10 Hz), both with amplitude 1 (i.e. signal = sin(2 * pi * t) + sin(10 * 2 * pi * t + pi/2)) This is such a simple signal that we can do the transform in our heads: there should be two peaks, each of amplitude 1, at 1 Hz and 10 Hz. But look at the output of fft(signal) in the bottom frame. It's not very useful.

Figure 12: Example of using matlab's FFT function as-is. Top: the input signal is the sum of a 1 Hz sine wave and a 10 Hz sine wave, both with amplitude 1. Bottom: the output signal is complex (real in blue, imaginary in green), is not scaled to the same units as the input, has a two-sided spectrum (i.e. four peaks instead of the expected two), and no x-axis frequency vector is provided.

The function below is simple wrapper script that corrects all of these deficiencies and gives a proper ``engineer's FFT''. I recommend that you use a function like this all the time, and never use matlab's fft() directly. Note: this function does not do any windowing. This topic is discussed later in this section. The result of using this function is shown in 13. See how much more useful that is?

% function [frq, amp, phase] = simpleFFT( signal, ScanRate)

% Purpose: perform an FFT of a real-valued input signal, and generate the single-sided 

% output, in amplitude and phase, scaled to the same units as the input.

%inputs: 

%    signal: the signal to transform

%    ScanRate: the sampling frequency (in Hertz)

% outputs:

%    frq: a vector of frequency points (in Hertz)

%    amp: a vector of amplitudes (same units as the input signal)

%    phase: a vector of phases (in radians)

function [frq, amp, phase] = simpleFFT( signal, ScanRate)

n = length(signal); 

z = fft(signal, n); %do the actual work

%generate the vector of frequencies

halfn = floor(n / 2)+1;

deltaf = 1 / ( n / ScanRate);

frq = (0:(halfn-1)) * deltaf;

% convert from 2 sided spectrum to 1 sided

%(assuming that the input is a real signal)

amp(1) = abs(z(1)) ./ (n);

amp(2:(halfn-1)) = abs(z(2:(halfn-1))) ./ (n / 2); 

amp(halfn) = abs(z(halfn)) ./ (n); 

phase = angle(z(1:halfn));

Figure 13: The same data from Figure 12 using the simple_fft script.

" . $row['Name'] . " Posted on " . $row['DateTime']; echo "
"; echo $row['Comment']; echo "

"; } echo ""; ?> Leave a comment on this page:

Name: (optional)
To prove you are not a robot, what is 2+3?

This work by Daniel Kiracofe (daniel dot kiracofe at gmail dot com) is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License./' $I