Add zero padding to fft for accurate estimation
显示 更早的评论
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;
采纳的回答
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)
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)
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)
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)
tempfreq = (0:N-1)/N*Fs;
try
tempfreq(length(x)/2+1)
catch ME
ME.message
end
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)
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)
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, the title of the Question asks about zero padding, but there is no zero-padding in the code.
2 个评论
David Goodmanson
2023-9-7
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
2023-9-7
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.
更多回答（1 个）
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)))
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)"
类别
在 帮助中心 和 File Exchange 中查找有关 Single-Rate Filters 的更多信息
Translated by 
