MATLAB - Fast Fourier Transform

MATLAB - Practical Fast Fourier Transform (FFT)

The Fast Fourier Transform (FFT) is a crucial algorithm for signal processing. Its primary purpose is to decompose a signal from its original domain (usually time) into its constituent frequencies. In essence, the FFT tells you what frequencies are present in your signal and what their amplitudes are.

MATLAB’s fft() function provides an efficient implementation of this algorithm. The key to using it effectively is not just calling the function, but understanding its output and how to plot it correctly against a meaningful frequency axis.

A Complete, Practical Example: Analyzing a Noisy Signal

Let’s walk through a complete, realistic example. We will create a signal containing two distinct sine waves, add some random noise, and then use the FFT to identify the original frequencies.

Step 1: Create the Time-Domain Signal

% Signal Parameters
Fs = 1000;            % Sampling frequency (Hz)
T = 1/Fs;             % Sampling period (s)
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector

% Create a signal containing a 50 Hz and a 120 Hz sinusoid
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

% Add zero-mean random noise to the signal
X = S + 2*randn(size(t));

% To visualize our noisy signal, we can plot it:
% plot(1000*t(1:50), X(1:50))
% title('Noisy Time-Domain Signal (First 50ms)')
% xlabel('t (milliseconds)')
% ylabel('X(t)')

At this point, just by looking at the time-domain plot, it would be very difficult to identify the 50 Hz and 120 Hz components.

Step 2: Compute the FFT

Now, we apply the fft function to our noisy signal X.

Y = fft(X);

Step 3: Compute the Frequency Spectrum

The output Y of the FFT is a complex-valued vector. To get the amplitude, we take its absolute value. The result of the FFT is also symmetric for real-valued inputs. We are typically only interested in the first half of the spectrum (the single-sided spectrum).

% Compute the two-sided spectrum P2.
P2 = abs(Y/L);

% Compute the single-sided spectrum P1 based on P2 and the even-valued signal length L.
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

Here, we divide by L to normalize the amplitude. We then take the first half of the points and multiply all but the first (DC) and last (Nyquist) components by 2 to conserve the signal’s energy.

Step 4: Define the Frequency Axis and Plot

To make the plot meaningful, we need to create a frequency vector that corresponds to the points in our spectrum P1.

% Define the frequency domain f
f = Fs*(0:(L/2))/L;

% Plot the single-sided amplitude spectrum P1
plot(f, P1)
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')
axis([0 200 0 1.2]) % Focus the view on our frequencies of interest

Expected Output Description: The resulting plot will show a frequency axis from 0 Hz onwards. Despite the heavy noise in the original signal, you will see two clear, sharp peaks. One will be located at 50 Hz with an amplitude of approximately 0.7, and the other will be at 120 Hz with an amplitude of approximately 1.0, successfully identifying the components of our original signal S.

Other fft Syntaxes

FFT with a Specified Length (n)

Y = fft(X, n) computes the n-point FFT. If X has fewer than n points, it is padded with zeros. If it has more, it is truncated. Zero-padding is a common technique to increase the frequency resolution of the FFT, making it easier to pinpoint peak frequencies.

% Using our original signal X (length 1500)
% Compute a 2048-point FFT (zero-padding)
Y_padded = fft(X, 2048);

FFT along a Specific Dimension

Y = fft(X, n, dim) computes the FFT along the dimension dim. This is useful for processing multiple signals stored in the columns or rows of a matrix.

% Create a matrix with 2 signals as columns
signal1 = sin(2*pi*50*t)';
signal2 = sin(2*pi*150*t)';
X_matrix = [signal1, signal2];

% Compute the FFT of each column (dim=1 is the default)
Y_matrix = fft(X_matrix);

% To compute the FFT of each row, you would use:
% Y_matrix_rows = fft(X_matrix, [], 2);

Important Considerations

• Inverse FFT (ifft): To transform a signal from the frequency domain back to the time domain, use the ifft function.
• fftshift: The raw output of fft places the zero-frequency (DC) component at the beginning of the array. For certain types of analysis and visualization, it’s useful to have the DC component in the center. The fftshift function accomplishes this.
• Spectral Leakage and Windowing: When a signal is not periodic within the measurement window, it causes its energy to “leak” into adjacent frequency bins, smearing the peaks. To mitigate this, you can multiply your signal by a windowing function (e.g., hann or hamming) before performing the FFT. This is a critical step in professional spectral analysis.