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;

Answer 1

[Paul]

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.

[David Goodmanson]

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]

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.

Answer 2

[dpb]

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:

Also note for <the output variables for fft> says

"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)"