How come the frequency vector starts at 0 for even length FFT ?
34 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Marguerite Marnat
el 1 de Abr. de 2016
Editada: Derek Alley
el 3 de Ag. de 2021
In many examples including the documentation about fft, the one-sided frequency vector is defined as: f = Fs*(0:(L/2))/L where L=length of the FFT. However it seems to me that for even-length FFT, the DC component is missing and so on. So you would plot at 0 the component computed for (+delta f /2 ) which is fairly ok because it is averaged of the 2 points around 0, but also the component plotted at frequency (+delta f) is the one computed for (+3.delta f/2), which is false. How do you explain that ?
And what would be the definition of a two-sided frequency vector for both even and odd lengths ? None of the answers provided on the forum are satisfying... Thank you for your help.
0 comentarios
Respuesta aceptada
Rick Rosson
el 2 de Abr. de 2016
Editada: Rick Rosson
el 2 de Abr. de 2016
Regardless of whether L is even or odd, the first element of the vector returned by fft is always the DC component. No matter what.
So, using the full frequency span, the frequency vector is computed as
L = length(X);
dF = Fs/L;
f = dF*(0:L-1)';
This code is correct for any length L, whether even or odd.
For the one-sided span, the code becomes:
f = dF*(0:floor(L/2)-1)';
For the two-sided span centered on DC, the code becomes:
f = (-Fs/2:dF:Fs/2-dF)' + mod(L,2)*dF/2;
Again, this code is correct for both even and odd L.
3 comentarios
Fernando Zigunov
el 13 de Ag. de 2020
I think his answer is correct, but perhaps this is a better way of interpreting.
The actual spectrum is always one-sided:
f = dF*(0:(L-1))';
But since it exceeds the Nyquist frequency (dF*L/2), the frequencies higher than (dF*L/2) will alias and become negative. In general:
f(f>(dF*L/2))=f(f>(dF*L/2))-(dF*L);
We're subtracting dF*L from all over-Nyquist frequencies, making them negative. This helps when doing some PDE solving operations, where the negative frequencies need to be used with their original sign.
Please correct me if I'm mistaken!
Derek Alley
el 3 de Ag. de 2021
Editada: Derek Alley
el 3 de Ag. de 2021
Hello,
I agree with Gladir, the one-sided solution by Rick is discarding the highest frequency bin for odd sized vectors.
Consider this example:
Lodd=7; Leven=8;
FFTodd = linspace(0,Lodd-1,Lodd); % fabricated FFT output for clarity
FFTeven = linspace(0,Leven-1,Leven); % fabricated FFT output for clarity
% the first element of the FFTs is the DC component
% use fftshift to put the DC component in the 'middle'
FFTodd_shift = fftshift(FFTodd);
FFTeven_shift = fftshift(FFTeven);
The FFTs and their shifted counterparts look like this:
FFTodd = 0 1 2 3 4 5 6
FFTodd_shift = 4 5 6 0 1 2 3
FFTeven = 0 1 2 3 4 5 6 7
FFTeven_shift = 4 5 6 7 0 1 2 3
In both cases, the positive freq portion (including DC) of the FFT contains 4 elements
Using:
0:floor(L/2)-1
to define the number of elements in the one-sided positive frequency axis will truncate the axis for odd length FFTs. For example, if L = 7, then floor(L/2)-1 = 2, and the positive frequency axis would only have 3 elements instead of the 4 that are actually there. Maybe that last frequency bin doesn't matter much in practical problems, but for correctness, and to match the output size of the positive frequency portion of fftshift(), ceil() would be a good alternative to floor().
For the one sided example, the formula becomes:
f = dF*(0:ceil(L/2)-1)'
On the other hand, Rick's formula for the 2-sided frequency axis is mostly correct, but I have found that the resulting frequency vector's DC value is not exactly 0, possibly due to floating point calculation errors. The two-sided formula provided by Fernando, which is basically the same as using fftshift on the frequency vector given by:
f = dF*(0:(L-1))'
and then subtracting Fs from only the negative frequencies, has always returned a value of 0 for the DC frequency component.
Más respuestas (0)
Ver también
Categorías
Más información sobre Spectral Measurements en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!