Fast Fourier Transform in MATLAB

[member 428835]

Hi PF!

I'm following a tutorial in MATLAB, shown here

t = 0:.001:.25;
x = sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t));
Y = fft(y,251);
Pyy = Y.*conj(Y)/251;
f = 1000/251*(0:127);
plot(f,Pyy(1:128))
title('Power spectral density')
xlabel('Frequency (Hz)')

I read the documentation but a few things are still unclear, which I'll list for clarity:
Line 4: why use the 251 part of the transform?
Line 5: Why this rather than, say Norm(Y)?
Line 6: so the 1000 is evidently from line 1, the time step? But where is the 251 coming from; looks like it's the time step * end time + 1.

Thanks!

 

[jedishrfu]

Heres what Mathworks has to say

https://www.mathworks.com/help/matlab/ref/fft.html
You should be able search this yourself. This page has some code to try. Matlab is a math lab for its users allowing to explore data numerically.

When we did this in school we created our data from the sum of two sin functions and then applied the fft to it to verify that the fft found the frequencies.

 

[member 428835]

Heres what Mathworks has to say

https://www.mathworks.com/help/matlab/ref/fft.html
You should be able search this yourself. This page has some code to try. Matlab is a math lab for its users allowing to explore data numerically.

When we did this in school we created our data from the sum of two sin functions and then applied the fft to it to verify that the fft found the frequencies.

Yea, I looked at this too. But I could not figure out what was meant by "n-point DFT" and when to specify n.

Also, can you tell me why they plot abs(Y/L) rather than simply abs(Y)?

 

[jedishrfu]

This explains it better with Matlab too

https://www.csun.edu/~skatz/ece460/matlab_tut_two.pdf

 

[RPinPA]

But I could not figure out what was meant by "n-point DFT"

A digital Fourier transform (DFT) is a transformation of n points. The equations are a sum from 1 to $n$. $n$ is the number of points you are transforming. The number of time values. And you will get out the same number of frequency values. If you do a 32-point transform, transforming 32 time samples, you'll get 32 frequency samples.

The fast implementation, the FFT, was originally designed for $n$ being a power of 2. The algorithm has been modified over the years so that's no longer required, but I believe it's still most efficient if $n$ is a power of 2.

As to your original questions:

Line 4: why use the 251 part of the transform?

Since $t$ in this example looks like it has 251 points, that argument is redundant. It will be a 251-point transform by default. But if you changed line 1 to have more or less points, the effect would be to always use a transform of the same length, which will give the same frequency resolution and would allow comparison.

Perhaps the tutorial experimented with changes in line 1 while keeping the rest the same?

Line 5: Why this rather than, say Norm(Y)?

Because norm(Y) would calculate $\sum_{i=1}^n |Y_i|^2$. This is calculating the vector of power values $|Y_i|^2$, not their sum.

Edit to add: abs(Y).^2 would have the same effect as Y .* conj(Y). So it's really a stylistic thing. Y . * conj(Y) makes it clear to some people (to me for instance, and apparently the author) exactly the nature of this calculation. No strong reason for it though. Sometimes you just write equations in a certain way as a note to yourself.

The division by $n$ has to do with the scaling factors that the FFT introduces. It's converting the power to some particular units. I don't know the details of what scale factors they're accounting for.

Line 6: so the 1000 is evidently from line 1, the time step? But where is the 251 coming from; looks like it's the time step * end time + 1.

The length of the original vector is $T = 251/1000$ seconds. That is the time-span represented by the original time sample. This causes the frequency spacing, the $\Delta f$ between frequency points, to be $1/T = (1000/251)$ Hz.

And the reason for the (0:127) is that when you DFT a real signal, the upper half is the mirror image of the lower half (well, actually its complex conjugate). This line is plotting 128 points of the magnitude of the FFT, just the lower half. Since it was a 251-point transform, that seems to me to be going a couple of points past the halfway mark.

Try changing it to (0:250) to see what happens. (Of course you need to modify line 7 so you use the same number of points of Pyy).