{"id":4410,"date":"2021-02-26T12:00:12","date_gmt":"2021-02-26T17:00:12","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/?p=4410"},"modified":"2021-03-06T09:29:10","modified_gmt":"2021-03-06T14:29:10","slug":"fft-spectral-leakage","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2021\/02\/26\/fft-spectral-leakage\/","title":{"rendered":"FFT Spectral Leakage"},"content":{"rendered":"
A MATLAB user recently contacted MathWorks tech support to ask why the output of <\/span>fft<\/span> did not meet their expectations, and tech support asked the MATLAB Math Team for assistance. Fellow Georgia Tech graduate Chris Turnes wrote a detailed response that I enjoyed reading. I thought it would be worth adapting the case for a blog post. Although the case is about the 1-D FFT, the underlying issues show up in image processing, too, and I have written about them in the past. (See my <\/span>Fourier transform category<\/span><\/a>.)<\/span><\/div>
Let's set up the case by starting with a sine wave at 50 Hz, sampled at 1 kHz, with 4080 samples.<\/span><\/div>
Fs = 1000;<\/span><\/span><\/div><\/div>
N = 4080;<\/span><\/span><\/div><\/div>
t = (0:(N-1))\/Fs;<\/span><\/span><\/div><\/div>
y = sin(2*pi*50*t);<\/span><\/span><\/div><\/div>
plot(t,y)<\/span><\/span><\/div><\/div>
axis([0 0.1 -1.2 1.2])<\/span><\/span><\/div><\/div>
xlabel(<\/span>'t (sec)'<\/span>)<\/span><\/span><\/div><\/div>
ylabel(<\/span>'y(t)'<\/span>)<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div><\/div>
title(<\/span>'The first 0.1 seconds of y(t)'<\/span>)<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
Next, let's compute and display the FFT, scaling the frequency axis so that it is in Hz, and scaling the magnitude by the square root of the FFT length.<\/span><\/div>
Y = fft(y);<\/span><\/span><\/div><\/div>
f = (Fs\/N) * (0:N-1);<\/span><\/span><\/div><\/div>
plot(f,abs(Y)\/sqrt(N))<\/span><\/span><\/div><\/div>
xlabel(<\/span>'Frequency (Hz)'<\/span>)<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
Let's zoom around the peak at 50 Hz.<\/span><\/div>
xlim([45 55])<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
This all seems reasonable, so far. Now, let's try the procedure using only the first 4096 samples of <\/span>$ y(t) $<\/span>. (That's approximately 4.1 seconds.)<\/span><\/div>
N2 = 4096;<\/span><\/span><\/div><\/div>
t2 = (0:(N2-1))\/Fs;<\/span><\/span><\/div><\/div>
y2 = sin(2*pi*50*t2);<\/span><\/span><\/div><\/div>
Y2 = fft(y2);<\/span><\/span><\/div><\/div>
f2 = (Fs\/N2) * (0:N2-1);<\/span><\/span><\/div><\/div>
plot(f2,abs(Y2)\/sqrt(N2))<\/span><\/span><\/div><\/div>
xlabel(<\/span>'Frequency (Hz)'<\/span>)<\/span><\/span><\/div><\/div>
xlim([45 55])<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
That looks a lot different. Let's plot them together.<\/span><\/div>
hold <\/span>on<\/span><\/span><\/div><\/div>
plot(f,abs(Y)\/sqrt(N))<\/span><\/span><\/div><\/div>
hold <\/span>off<\/span><\/span><\/div><\/div>
legend([<\/span>\"|Y_2|\"<\/span>,<\/span>\"|Y|\"<\/span>])<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
Hmm. The FFT of the signal with 4080 points has a clean, sharp peak at 50 Hz, whereas the FFT of the signal with 4096 points is much more spread out around that peak.<\/span><\/div>
This was the tech support question: What is the explanation for the difference? Is there something wrong with the <\/span>fft<\/span> function for <\/span>$ N = 4096 $<\/span>?<\/span><\/div>
Chris immediately provided a nice, conceptual explanation of what's going on here, and I'll get to that in a just bit. First, though, I want to go through a thought process I sometimes use when faced with similar questions. Is there an independent way to verify whether the function in question is returning the correct answer? Well, in this case, I know that the <\/span>fft<\/span> function computes something called the <\/span>discrete Fourier transform<\/span>, or DFT. For a discrete-time signal <\/span>$ y[n] $<\/span>, defined over <\/span>$ 0 \\leq n < N $<\/span>, the DFT is given by:<\/span><\/div>
$ Y[k] = \\sum_{n=0}^{N-1} y[n] e^{-j 2 \\pi k n \/ N}, k = 0, 1, \\ldots, N-1 $<\/span><\/div>
The <\/span>fft<\/span> function doesn't compute the DFT using this formula because it is slow, but we can use the formula to double-check its output.<\/span><\/div>
k = (0:N2-1);<\/span><\/span><\/div><\/div>
n = k';<\/span><\/span><\/div><\/div>
Y2p = sum(y2' .* exp(-1i * 2 * pi * k .* n \/ N2));<\/span><\/span><\/div><\/div><\/div>
Y2<\/span> and <\/span>Y2p<\/span> are the same to within about 12 decimal digits of relative precision:<\/span><\/div>
norm(Y2p - Y2)\/norm(Y2)<\/span><\/span><\/div>
ans = 7.9563e-13<\/div><\/div><\/div><\/div>
That is well within the difference I would expect to see from ordinary floating-point round-off differences. But let's plot them together, just to be sure.<\/span><\/div>
plot(f2,abs(Y2)\/sqrt(N2),f2,abs(Y2p)\/sqrt(N2))<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div><\/div>
xlim([45 55])<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
The two results are visually indistinguishable.<\/span><\/div>
OK, at this point, then, I would be confident that <\/span>fft<\/span> is returning the correct answer. We can turn to the question of why the output appears to be so different just because we changed the total number of samples slightly. Let me hand it over to Chris, who provided following explanation:<\/span><\/div>
This is just <\/span>spectral leakage<\/span><\/a>; it's a well-known and <\/span>well-described<\/span><\/a> problem in Fourier analysis. The discrete spectral frequencies sampled by the DFT of length <\/span>N<\/span> are <\/span>(0:(N-1))\/N<\/span>. In the customer's case, the discrete frequency of his pure tone is not an integer multiple of <\/span>$ 1\/N $<\/span>. As a result, the frequency content \"leaks\" into other bins.<\/span><\/div>
Another way of seeing this is to recall that the DFT\/FFT inherently assumes the input is a periodic signal. However, neither of the signals above, <\/span>y<\/span> and <\/span>y2<\/span>, contains an integer number of periods. You can see this if you plot the end of the signals instead of the beginning.<\/span><\/div>
plot(t((end-50):end),y((end-50):end))<\/span><\/span><\/div><\/div>
ylim([-1.2 1.2])<\/span><\/span><\/div><\/div>
title(<\/span>'N = 4080'<\/span>)<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
plot(t2((end-50):end),y2((end-50):end))<\/span><\/span><\/div><\/div>
ylim([-1.2 1.2])<\/span><\/span><\/div><\/div>
title(<\/span>'N = 4096'<\/span>)<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
<\/span><\/div>
With <\/span>$ N = 4080 $<\/span>, the signal has just about completed its period at the end, so the last sample is almost equal to the first sample (which is 0). Therefore, the artifacts from the assumed periodicity aren't that strong. With 4096 samples, though, the signal ends with a quarter of its period still to go. There's a big jump from the last sample value (-1) to the first sample value (0). This big jump is like a discontinuity in the signal, and that spreads the spectral peak, resulting in more \"leakage.\" <\/span><\/div>
If everything lines up, then life becomes beautiful; with any integer <\/span>$ k\\geq 2 $<\/span>, using <\/span>$ N=50k $<\/span> will result in the ideal spike we'd all like to see. For other lengths, however, where <\/span>$ 50\/1000 $<\/span> can't be written as <\/span>$ q\/N $<\/span> for an integer <\/span>q<\/span>, you'll get the spectral leakage phenomenon. The amplitude will dip and the peak will spread out. The result is not incorrect; it is the mathematically expected result, even if it does not match our initial intuition.<\/span><\/div>
Thanks for letting me use your explanation, Chris!<\/span><\/div><\/div>