Effacer les filtres
Effacer les filtres

Add zero padding to fft for accurate estimation

22 vues (au cours des 30 derniers jours)
Niusha
Niusha le 6 Sep 2023
Commenté : Paul le 7 Sep 2023
Hi,
I wonder why in this line of code(bolded), there is the length divided by 2+1?
Fs = 1e3;
t = 0:0.001:1-0.001;
x = cos(2*pi*100*t)+sin(2*pi*202.5*t);
xdft = fft(x);
xdft = xdft(1:length(x)/2+1);
xdft = xdft/length(x);
xdft(2:end-1) = 2*xdft(2:end-1);
freq = 0:Fs/length(x):Fs/2;

Réponse acceptée

Paul
Paul le 6 Sep 2023
Modifié(e) : Paul le 6 Sep 2023
Hi Niusha
Here is the first part of the code:
Fs = 1e3;
t = 0:0.001:1-0.001;
x = cos(2*pi*100*t)+sin(2*pi*202.5*t);
xdft = fft(x);
Here is the frequency vector that corresponds to all of the elements fo xdft
N = numel(xdft)
N = 1000
tempfreq = (0:N-1)/N*Fs;
Because N is even, the value of tempfreq that corresponds to length(x)/2 + 1 is the Nyquist frequency (Fs/2)
tempfreq(length(x)/2+1)
ans = 500
So the goal for the remainder of the code is to only use the elements of xdft that correspond to frequencies from 0 to Fs/2
xdft = xdft(1:length(x)/2+1);
xdft = xdft/length(x);
xdft(2:end-1) = 2*xdft(2:end-1);
freq = 0:Fs/length(x):Fs/2;
freq(end)
ans = 500
However, if x has an odd number of points, then this code needs to be modified
Fs = 1e3;
t = 0:0.001:1;
x = cos(2*pi*100*t)+sin(2*pi*202.5*t);
xdft = fft(x);
N = numel(xdft)
N = 1001
tempfreq = (0:N-1)/N*Fs;
try
tempfreq(length(x)/2+1)
catch ME
ME.message
end
ans = 'Array indices must be positive integers or logical values.'
That same indexing won't work for xdft either.
For N odd, it's common to use the frequencies from 0 to
tempfreq((length(x)-1)/2 + 1)
ans = 499.5005
And then
xdft = xdft(1:(length(x)-1)/2 + 1);
xdft = xdft/length(x);
Because the last point does not correspond to Fs/2, it should be muliplied by 2
xdft(2:end) = 2*xdft(2:end);
And the frequency vector is most easily computed as
%freq = 0:Fs/length(x):Fs/2;
freq = (0:(length(x)-1)/2)/N*Fs;
freq(end)
ans = 499.5005
Basically, if you're going for the "one-sided" DFT, you have to be careful if the length of xdft is even or odd.
Also, consider using numel instead of length
Also, the title of the Question asks about zero padding, but there is no zero-padding in the code.
  2 commentaires
David Goodmanson
David Goodmanson le 7 Sep 2023
Hi Niusha,
let y = fft(x), n = length(y) sampling frequency = Fs
Then to obtain values of y corresponding to positive (and 0) frequencies, the following expressions work for both n even and n odd.
frequencies: (0:floor(n/2))*(Fs/n) % safer with integer increment than
% floating point increment
y values: y(1:floor(n/2)+1)
After that, as Paul pointed out, for n even, all but the first and last (0 frequency and nyquist) points get multiplied by 2. For n odd, all but the first (zero frequency) point gets multiplied by 2. For odd n, the nyquist frequency is not present in the frequency grid.
soapbox time: It's ok for plotting purposes I suppose, but I am not a big fan of the multiply by 2, positive frequency aproach, and definitely not if the lower half of the y array is discarded. That makes it inconvenient to do an inverse transform and increases the possibility of mistakes when doing so. If the absolute value of y is taken and the phase is discarded, then one can't do the inverse transform at all. I think it's much preferable to keep the entire array around and access parts of it as needed.
Paul
Paul le 7 Sep 2023
Is there room for two on your soapbox?
Another problem with the "one-sided" approach is that it's most common, at least on this forum, to plot the spectrum with a plain-old line plot. Because the first point is not multiplied by two, the interval of the plot line between the first and second points of the spectrum is not correct.
Also, the one-sided-multiply-by-two approach is a view on magnitude of the time-domain signal, but loses the energy relationship between the time domain and frequency domain.
Also, the one-sided-multiply-by-two approach does not work at all for complex input signals.
I guess it all depends on what one is trying to show and how one is trying to show it.

Connectez-vous pour commenter.

Plus de réponses (1)

dpb
dpb le 6 Sep 2023
Modifié(e) : dpb le 6 Sep 2023
Returning the one-sided frequency spectrum instead of two-sided (the zero frequency/DC component is in the middle).
Fs = 1e3;
t = 0:0.001:1-0.001;
x = cos(2*pi*100*t)+sin(2*pi*202.5*t);
xdft = fft(x);
subplot(2,1,1)
plot(abs(xdft))
tmp=xdft;
title('Two Sided DFT Magnitude vs Index')
xdft = xdft(1:length(x)/2+1);
subplot(2,1,2)
plot(abs(xdft))
title('One Sided DFT Magnitude vs Index')
Nota Bene: the second of the two peaks is roughly half in magnitude of the first (not exact because binning isn't exact for 202.5 Hz so energy is spread over several bins as the shape shows). By happenstance of the sampling rate and length chosen, the frequency and bin number are equivalent for the first half.
You can see the symmetry also if look at the values
L=length(tmp); M=L/2+1;
fprintf('%0.8f ',real(tmp(M-4:M+4)))
-0.73879176 -0.73871287 -0.73865653 -0.73862273 -0.73861146 -0.73862273 -0.73865653 -0.73871287 -0.73879176
You can observe the values are symmetric about the midpoint, M.
ADDENDUM:
"Y - Frequency domain representation returned as a vector, matrix, or multidimensional array.
...
If X is real, then Y is conjugate symmetric, and the number of unique points in Y is ceil((n+1)/2)"

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by