Zero-Padding Position for FFT

Hi,

I have a time-domain signal x with 5000 samples with period T, and it spans from -0.5*T to +0.5*T.

I want to compute an 8000-point FFT and I am confused about where I should add the zero-padding. Which of the following two methods is correct?

First:

x = [x; zeros(1, 3000)];
X = fftshift( fft( fftshift(x), 8000) );

Second:

x = [zeros(1, 1500); x; zeros(1, 1500)];
X = fftshift( fft( fftshift(x), 8000) );

Thanks!

Answer 1

[Paul]

Hi Jason,

The data in x can't span exactly -0.5*T to 0.5*T if it has an even number of uniformly spaced samples.

Assuming the data actually starts at 0.5*T with 5000 samples, we have (for example)

T = 1;
N = 5000;
Ts = T/N;
n = (0:N-1) - N/2;
t = n*Ts;
format short e
t([1 end])  % end points are not exactly the same with N even

ans = 1×2

-5.0000e-01 4.9980e-01

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% example signal over (nearly) one period
x = sin(2*pi/T*t);

For this problem, want to zero pad the left and right edges as in your second option, but use ifftshift on the inside

Npad = 3000;
X = fftshift(fft(ifftshift([zeros(1,Npad/2),x,zeros(1,Npad/2)])));
% plot
Nfft = numel(X);
f = ceil(-Nfft/2:(Nfft/2-1))/Nfft/Ts;
figure
stem(f,abs(X),'-o');

●

figure
stem(f,abs(X),'-o');
xlim([-10,10])

Warning: Hardware-accelerated graphics is unavailable. Displaying fewer markers to preserve interactivity.

[Jason]

Thank you very much for the clarification. But I am confused about the difference between my first and second method, could you help me understnad why padding zero in the two sides are correct?

[Paul]

Well, I guess it really depends on what we're trying to do. I was assuming that we wanted to get finer sampling of the amplitude and phase of the Discrete Time Fourier Transform (DTFT) of the the signal.

If the goal is to only capture the amplitude, then either of your approaches will work and, as @Matt J noted, there is no need to worry about shifting on input to fft.

T = 1;
N = 50;  % use smaller N and Npad to make the plots less cluttered
Ts = T/N;
n = (0:N-1) - N/2;
t = n*Ts;
x = sin(2*pi/T*t);
Npad = 30;
x1 = [x,zeros(1,Npad)];
x2 = [zeros(1,Npad/2),x,zeros(1,Npad/2)];
X1 = fftshift(fft(ifftshift(x1)));  % option 1
X2 = fftshift(fft(ifftshift(x2)));  % option 2
X3 = fftshift(fft(x1));             % no shifting 
Nfft = numel(X1);
f = ceil(-Nfft/2:(Nfft/2-1))/Nfft/Ts;

The amplitudes of all three are the same

figure
stem(f,abs(X1),'o');hold on
stem(f,abs(X2),'x')
stem(f,abs(X3),'d')

But the phases are not

figure
stem(f,angle(X1),'o'),hold on
stem(f,angle(X2),'x')
stem(f,angle(X3),'d')

The DTFT of x is (computed for 4096 frequency points)

[xDTFT,wDTFT] = freqz(x,1,linspace(-pi,pi,4096));
xDTFT = xDTFT.*exp(-1j*wDTFT*n(1)); % freqz assumes first element of x is at n = 0, so we have to adjust the phase 

Comparing its real and imaginary parts to X2 from above

figure
plot(wDTFT/Ts/2/pi,real(xDTFT)),grid % floating point trash
hold on
stem(f,real(X2))
xlim([-20 20]) 

figure
plot(wDTFT/Ts/2/pi,imag(xDTFT)),grid
hold on
stem(f,imag(X2))
xlim([-20,20])

Neither X1 nor X3 will match the real and imaginary parts of xDTFT, but they will both match its magnitude. So if magnitude is all that matters, either X1 or X2 or X3 will do.

[Jason]

Thank you very much!

[Matt J]

Thank you very much for the clarification. But I am confused about the difference between my first and second method, could you help me understnad why padding zero in the two sides are correct?

You do not have to pad on both sides, even if you do care about getting both amplitude and phase spectra. My answer shows an alternative that works with 1-sided padding. The important thing is that you properly shift the zero-padded signal so that the signal origin moves to the beginning array.

Note also that both of your methods will fail if you were to pad the array to an odd length. In that case, fftshift and ifftshift are not equivalent:

ifftshift(1:7)

ans = 1×7

4 5 6 7 1 2 3

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fftshift(1:7)

ans = 1×7

5 6 7 1 2 3 4

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Paul]

To illustrate that two-sided padding followed by the ifftshift* is the same as @Matt J's approach of one-sided padding followed by a circular shift ... They both result in the same input to fft.

*Throughout this discussion I've been assuming that the original sequence spans -N/2 to N/2-1 for N even, or -(N-1)/2 to (N-1)/2 for N odd, and that we are padding with an even number of zeros (as in the original question).

T = 1;
N = 50;  % use smaller N and Npad to make the plots less cluttered
Ts = T/N;
n = (0:N-1) - N/2;
t = n*Ts;
x = sin(2*pi/T*t);
Npad = 30;
x2 = [zeros(1,Npad/2),x,zeros(1,Npad/2)];

n=floor(numel(x)/2)+1; %location of the signal origin
x(end+1:80)=0; %zero pad to 80
xs=circshift(x, 1-n); %shift so that signal origin is at xs(1)
stem(subplot(211),ifftshift(x2));
stem(subplot(212),xs)

[Paul]

As stated by @Matt J,

"The important thing is that you properly shift the zero-padded signal so that the signal origin moves to the beginning [of the] array."

Zero-pad first, and then (correctly) shift.

This point is so important and well stated that I thought it useful to illustrate it in more depth.

Let's define a finite duration signal of length N = 9

x = @(n) interp1((-10:-2),[1,8,6,10,3,2,8,1,5],n,'nearest',0);

Plotting x[n] over -30 <= n <=30 graphically shows that x[n] = 0 for n < -10 and n > -2

figure
nv = -30:30;
stem(nv,x(nv)),xlabel('n');ylabel('x[n]')

The Discrete-Time Fourier Transform (DTFT) of x[n] is defined as

The frequency variable, w, has units of rad/sample and covers the entire real line. The DTFT is periodic, with period = 2*pi.

We can compute the DTFT directly at any frequency because we only need to sum over -10 <= n <= -2.

Define the frequencies over one-period of the DTFT from 0 <= w < 2*pi

w = (0:4095).'/4096*2*pi;

and compute the DTFT directly based on its definition.

n1 = -10:-2;
h1 = sum(exp(-1j*w.*n1).*x(n1),2);

We can also use freqz to compute the DTFT of x[n], where again we only need to use x[-10 <=n <= -2]. However, freqz() assumes that the first point in the input signal is x[0], whereas ours will be x[-10]. So we need to use use the time shifting property of the DTFT to correct the phase of the output from freqz().

h2 = freqz(x(n1),1,w).*exp(-1j*w*n1(1));

An option for computing frequency-domain samples of the DTFT is to use the Discrete Fourier Transform (DFT), which is implemented by fft. Like freqz(), fft() assumes that the first point of the input corresponds to n = 0. So we have to "trick" fft() into giving us what we want.

One way to understand the "trick" is to realize that the DFT outputs are also Discrete Fourier Series (DFS) coefficients of a periodic signal. Here, the periodic signal in question is formed by taking the N-periodic summation of x[n], which is basically appending an infinite number of replicates of x[n] to x[n]. Plotting a few periods of the periodic summation

N1 = numel(n1);
figure
stem(nv,x(nv+3*N1)+x(nv+2*N1)+x(nv+N1)+x(nv)+x(nv-N1)+x(nv-2*N1)+x(nv-3*N1)+x(nv-4*N1)+x(nv-5*N1))

and plotting the original signal to show where it lies within the periodic summation

hold on;
stem(nv,x(nv),'rx')

Just like any other Fourier series, DFS coefficients can be computed from any period of the periodic signal. However, fft() always computes the coefficients from the period that starts at n = 0 and ends at n = N-1. This period can be found using circshift, which is plotted in green

stem(0:N1-1,circshift(x(n1),n1(1)),'gx');
xlabel('n'),ylabel('x_P[n]')

Therfore, we can compute the DTFT samples by taking the DFT of the points in green.

h3 = fft(circshift(x(n1),n1(1)));

The DFT outputs from fft() are computed at frequencies

w3 = (0:N1-1)/N1*2*pi;

If we want finer sampling in the frequency domain we can zero-pad. We can zero-pad at either end or both ends, as long we we correctly shift after padding, as stated at the outset.

Let's zero-pad to a 31-point signal, with five zeros to the left and 17 zeros to the right.

n2 = -15:15; N2 = numel(n2);

Circular shift and then compute the DFT

h4 = fft(circshift(x(n2),n2(1)));w4 = (0:N2-1)/N2*2*pi;

Put all four responses on a single plot, and we see that both methods for computing the DTFT are the same and that both DFTs yield frequency domain samples of the DTFT.

figure
plot(subplot(211,'Nextplot','add'),w,abs(h1),w,abs(h2));stem(w3,abs(h3));stem(w4,abs(h4));
legend('DTFT','freqz','DFT-9','DFT-31','location','north');
ylabel('magnitude')
plot(subplot(212,'Nextplot','add'),w,angle(h1),w,angle(h2));stem(w3,angle(h3));stem(w4,angle(h4));
ylabel('phase (rad)');xlabel('\omega (rad/sample)')

Plot the circularly shifted signal for the 31-point case.

figure
stem(0:N2-1,circshift(x(n2),n2(1)))

Circular shifting is a general approach, but the same can be accomplished with ifftshift if the signal to be DFT'd is ifftshift-symmetric, i.e., a) has an equal number of points the left and right of n = 0 if N is odd, or b) one more point to the left of zero than to the right if N is even. Our 31-point case above is ifftshift-symmetric, and ifftshift() yields the same result as circshift above (in fact, ifftshift and fftshift are both implemented in terms of circshift() ).

hold on
stem(0:N2-1,ifftshift(x(n2)),'x')
xlabel('n');ylabel('x_{shifted}[n]');legend('circshift','iffshift','Location','northwest')

In summary: zero-pad first and then "properly shift the zero-padded signal so that the signal origin moves to the beginning [of the] array" if you want to correctly compute the phase in the frequency domain.

[Jason]

Thank you for the clear explanation and help! I understand it now!

Answer 2

[Matt J]

If you only care about the amplitude spectrum the shifts don't really matter, but otherwise, you would do,

n=floor(numel(x)/2)+1; %location of the signal origin

x(end+1:8000)=0; %zero pad to 8000

xs=circshift(x, 1-n); %shift so that signal origin is at xs(1)

X = fftshift( fft( xs ) ); %transform

[Matt J]

If you only care about the amplitude spectrum the shifts don't really matter

And in that case, it simplifies to,

x(end+1:8000)=0; %zero pad to 8000

X=fftshift( abs(fft(x)) );

[Paul]

Or just

X = fftshift(abs(fft(x,8000)))

so as to not corrupt the original data vector.