Steve on Image ProcessingConcepts, algorithms & MATLAB

This is machine translation

Translated by
Mouseover text to see original. Click the button below to return to the English version of the page.

Sinusoids and FFT frequency bins13

Posted by Steve Eddins,

Contents

Note: I really didn't expect to be writing about the FFT two weeks in a row. But ready or not, here we go!

My friend Joe walked into my office the other day and announced, "Steve, I've got a question about the FFT."

(Joe was one of my MathWorks interviewers 21 years ago this month. His question for me then was, "So what do you bring to the party?" I don't remember how I answered.)

Joe is one of the biggest MATLAB experimenters I know. He's always trying to use MATLAB for some new kind of analysis or to automate something dubious (such as changing the hue of the lights in his house. Hmm.)

Anyway, lately he's been into automatic analysis of music. And something was bothering him. If you have a pure sinusoid signal and take the FFT of it, why does it spread out in the frequency domain? How do you interpret the output of fft when the frequency is in between two adjacent frequency bins of the FFT? And is that situation fundamentally different somehow that when the sinusoid's frequency is exactly one of the FFT bin frequencies?

Oh boy. I told Joe, "Well, you don't really have a sinusoid. A sinusoid is an infinitely long signal. You have a finite portion of a sinusoid, which I think of as a rectangularly windowed sinusoid. Go away, let me think more about how to explain all that, and then I'll come find you tomorrow."

Because of my background in digital signal processing, I tend to think about this in terms of the relationships between three kinds of Fourier transforms:

• Fourier transform (continuous time - continuous frequency)
• discrete-time Fourier transform (DTFT) (discrete time - continuous frequency
• discrete Fourier transform (DFT) (discrete time - discrete frequency)

I have written before (23-Nov-2009) about the various kinds of Fourier transforms. It is this last kind, the DFT, that is computed by the MATLAB fft function.

It's hard to tell exactly where to start in any discussion about Fourier transform properties because the use of terminology and the mathematical convensions vary so widely. Let me just dive in and walk through Joe's example. If you have questions about what's going on, please post a comment.

Note: In what follows, I am making a very big simplifying assumption. Specifically, I'm assuming that the sampling frequency is high enough for us to be able to discount the possibility of aliasing in the sampling of the sinusoid. This isn't a completely valid assumption because, unlike a real sinusoid, the truncated or windowed sinusoid doesn't have a finite bandwidth.

Here some sample parameters that Joe gave me.

Sinusoid frequency in Hz:

f0 = 307;


Sampling frequency in Hz:

Fs = 44100;


Sampling time:

Ts = 1/Fs;


How many samples of the sinusoid do we have? From a theoretical perspective, think of this as the length of a rectangular window that multiplies our infinitely long sinusoid sequence.

M = 1500;


Construct the signal

n = 0:(M-1);
y = cos(2*pi*f0*Ts*n);


Compute M-point discrete Fourier transform (DFT) of f. This is what the function fft computes.

Y = fft(y);


Plot the DFT. The DFT itself doesn't have physical frequency units associated with it, but we can plot it on a Hz frequency scale by incorporating our knowledge of the sampling rate. (I'm plotting the DFT as disconnected points to emphasize its discrete-frequency nature. Also, I'll be plotting magnitudes of the Fourier transforms throughout.)

f = (0:(M-1))*Fs/M;
plot(f,abs(Y),'*')
xlim([0 Fs/2])
title('Magnitude of DFT')
xlabel('Frequency (Hz)')


Exploring the "spread" around the peak

Zoom in on the peak.

p = 2*Fs/(M+1);
f1 = f0 - 10*p;
f2 = f0 + 10*p;
xlim([f1 f2])


What's with the "spread" around the peak? The spread is caused by the fact that we do not have a pure sinusoid. We have only a finite portion of a pure sinusoid.

To explore this further, let's relate the DFT to what some might call the "real" Fourier transform of our windowed sinusoid. If we can discount aliasing, then the discrete-time Fourier transform (DTFT) is the same as the "real" Fourier transform up to half the sampling frequency. And, for a signal that is nonzero only in the interval 0 <= n < N, an M-point DFT exactly samples the DTFT as long as M >= N. The DTFT of an N-point windowed sinuoid has the form a bunch of copies of $sin(\omega (N+1)/2) / sin(\omega/2)$, where the copies are spaced apart by multiples of the sampling frequency.

Can we get a better idea of the true shape of the DTFT? Yes, by zero-padding the input to the DFT. This has the effect of sampling the DTFT more closely. We can use this relationship to plot a pretty good approximation of the DTFT and then superimpose the DFT (the output of fft).

M10 = length(y)*10;
Y10 = fft(y,M10);
f10 = (0:(M10-1))*Fs/M10;

plot(f10,abs(Y10));
hold on
plot(f,abs(Y),'*')
hold off

xlim([f1 f2])

legend({'DTFT','DFT'})
xlabel('Frequency (Hz)')


Moving the sinusoid frequency to line up with a bin

The DFT samples the Fourier transform at a spacing of Fs/M. Let's pick a frequency that lines up exactly on a bin.

f0 = 294;
w0 = 2*pi*f0/Fs;
n = 0:(M-1);
y = cos(w0*n);
Y = fft(y);
p = 2*Fs/(M+1);
f1 = f0 - 10*p;
f2 = f0 + 10*p;
plot(f,abs(Y),'*')
xlim([f1 f2])
title('DFT')
xlabel('Frequency (Hz)')


The DFT samples (output of fft) appears to be perfect above. It has just one nonzero value at the sinusoid frequency.

But the appearance of perfection is just a kind of illusion. You can see what's really happening by superimposing the DFT samples on top of the DTFT, which we compute as before using the zero-padding trick.

Y10 = fft(y,M10);
plot(f10,abs(Y10));
hold on
plot(f,abs(Y),'*')
hold off
xlim([f1 f2])
title('DTFT with DFT superimposed')
xlabel('Frequency (Hz)')


You can see that most of the DFT samples happen to line up precisely at the zeros of the DTFT.

Using a longer signal

Finally, let's look at what happens when you use a longer signal. (Joe said, "Can't I make the spread go away by using a sufficiently large number of samples of the sinusoid?" The answer, as we'll see demonstrated, is that you can make the spread get smaller, but you can't make it go away entirely as long as you have a finite number of samples.) Let's go up a factor of 10, to 15,000 samples of the sinusoid. (And let's go back to the original frequency of 307 Hz.)

f0 = 307;
M = 15000;
w0 = 2*pi*f0/Fs;
n = 0:(M-1);
y = cos(w0*n);
Y = fft(y);

M10 = length(y)*10;
Y10 = fft(y,M10);
f10 = (0:(M10-1))*Fs/M10;
f = (0:(M-1))*Fs/M;

plot(f10,abs(Y10));
hold on
plot(f,abs(Y),'*')
hold off

xlim([f1 f2])

title('DTFT with DFT superimposed')
xlabel('Frequency (Hz)')


The plot above is shown at the same scale as the one for M = 1500 further up. You can see that the peak is much narrower, but it's hard to see the details. Let's zoom in closer.

p = 2*Fs/(M+1);
f1 = f0 - 10*p;
f2 = f0 + 10*p;
xlim([f1 f2])


Now you can see that the DTFT shape and DFT samples look very much like what we got for the 1500-point windowed sinusoid. The width of the peak, as well as the spacing of the sidelobes, has been reduced by a factor of 10 (15000 / 1500). And the spacing (in Hz) of the DFT samples has been reduced by the same amount.

Summary

In summary, the Fourier transform of a finite portion of a sinusoid has a characteristic shape: a central peak (the "main lobe") with a bunch of "side lobes" in either direction from the peak. The DFT simply samples this shape. Depending on the exact frequency, the output of the DFT might look a little different, depending on where the frequency happens to line up with the DFT frequency bins.

Get the MATLAB code

Published with MATLAB® R2014a

Note

Cris Luengo replied on : 1 of 13

Dear Steve,

You are usually spot-on in your blog posts, but here you seem to make a mistake.

The DFT is not just samples from the DTFT. In the DTFT one assumes an infinitely long signal, but the DFT takes a finite signal, M samples. It is then assumed that these samples are repeated to infinity (sample M is equal to sample 0, M+1 equal to 1, etc). Same as sampling the spatial domain causes the frequency domain to be finite and repeated, making the spatial domain finite and repeated causes the frequency domain to be sampled. Simply put, in the discrete signal all periods are integer divisions of M: M, M/2, M/3, M/4, M/5, etc. You cannot create a signal with a frequency of 0.99 M, because of the forced repetition with period M. In the signal in your example, if you look at what happens at the transition between sample M-1 and M (or -1 and 0), that is at the boundary where we repeat the samples, you will see a sharp transition. This transition has many frequencies, and hence the single peak at f0 looks spread out.

When you pad the DFT with zeros and do the IFFT, you end up with a signal interpolated to the new size, with a sinc interplating kernel. In the same way, if you pad your signal with zeros, you interpolate the DFT with a sinc. This is where the lobes come from that you show. This is an artefact of the padding, and doesn't represent frequency content of the original signal. If you compute the FFT of repmat(y,1,10) you will get a "higher resolution" DFT.

Steve Eddins replied on : 2 of 13
Cris—I did not clearly express how I was thinking about Joe's question. I told him, "You don't have a sinusoid. You have a finite portion of a sinusoid, or a rectangularly windowed sinusoid." So, as you know (but as I didn't say well enough), for a signal that is nonzero only in the interval 0 <= n < N, an M-point DFT exactly samples the DTFT of the windowed sinusoid as long as M >= N. In turn, the DTFT is the convolution of an impulse train with sin(\omega (M+1)/2) / sin(\omega / 2). I'll look over my post again and edit it to try to make this point of view more clear.
Steve Eddins replied on : 3 of 13
Cris—I made several edits. I included my comment to Joe that he didn't really have a true sinusoid, which is infinitely long, but instead a rectangularly windowed sinusoid. I noted that my assumption about avoiding aliasing wasn't completely valid because a windowed sinusoid does not have finite bandwidth. I indicated the conditions under which the DFT exactly samples the DTFT. And I explained the form of the DTFT for a windowed sinusoid.
Bjorn Gustavsson replied on : 4 of 13
Chris, The repmat would give you a signal with discontinuities at the "seams" for all cases except when the period of your signal matches the original sampling window. This would also not give you a flawlesser Fourier transform...
Eric replied on : 5 of 13
Great post as usual Steve :) Maybe you could follow it up with a treatment of the same topic, but using the Goertzel algorithm and discussing what advantages it has? I've found from talking to 'regular folks' who need to do DSP stuff that it's an under-appreciated tool, compared to the more widely known FFT/DFT, but it has big advantages in some not-uncommon scenarios. Maybe at the same time I can plug a request to add the Generalized Goertzel alg to the SPT :D
BolderMan replied on : 6 of 13
ItsWasAlwaysHardForMeToMakeThroughThisKindOfMathematicalFunctionsAndManyNumbers!!!
Steve Eddins replied on : 7 of 13
Eric—Thanks for the suggestion! I have forwarded your enhancement request along to the Signal Processing Toolbox development team.
Eric Sampson replied on : 8 of 13
Thanks Steve!
Yusuf replied on : 9 of 13
Its been more than 6 years since I joined this blog. I have not logged on since 3 or 4 years from today. I see its a lot of change. I finished my master thesis in 2010 in image processing and I eager no to continue my phd in the same manner, the dust has filled my shelves I guess but i will try my best to achieve my goal with the aid of this great blog that it has not cut its connection with me even I have not participated nor respond to its mail feeds since years, THANKs Yusuf
Greg replied on : 10 of 13
Choose M = 44100.
Jim replied on : 11 of 13
Dear Steve, Excellent post! Tell Joe he will be able to control the spread (and the main lobe) by multiplying the data with a window function (e.g. Hanning, Hamming). Best, Jim
Kevin replied on : 12 of 13
I think Chris' comment is excellent, I would just like to add another (perhaps more intuitive) way of thinking about it. Say we have some true signal sin(w*t) that persists for infinite time - this has a sharp peak in the FT. When we sample only a limited time of it, and then perform the DFT, you can think of as we are multiplying it by a square function in time and Fourier transforming that. In the frequency space then, this is a convolution of our delta function (FT of sin(w*t)) with a sinc function (FT of the square 'window'). Then what you're seeing is the broadening by the sinc function. As you sample longer ,the sinc function gets narrower. Another way to look at this is the time-frequency uncertainty inherent in the DFT - if you sample for shorter times, you have more uncertainty in the frequency - hence a broader spectral peak. As you show with the "DTFT" this doesn't go away if the frequency happens to line up with a bin, rather it samples it such that it just appears to be a delta function.
muna shehan replied on : 13 of 13
Hi; Pleas can anyone explain how i can determine sinusoid frequency if I have a step response where the data form is two vectors one for the acceleration and the other is timeز