Zero-Padding Position for FFT

Question

Jason 2025-7-2

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2178202-zero-padding-position-for-fft

评论： Jason 2025-7-7

采纳的回答： Paul

在 MATLAB Online 中打开

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!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Paul 2025-7-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2178202-zero-padding-position-for-fft#answer_1567464

在 MATLAB Online 中打开

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.

7 个评论
显示 5更早的评论隐藏 5更早的评论

Paul 2025-7-3

编辑：Paul 2025-7-3

在 MATLAB Online 中打开

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.

Paul 2025-7-6

编辑：Paul 2025-7-6

在 MATLAB Online 中打开

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 2025-7-7

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

请先登录，再进行评论。

Answer 2

Matt J 2025-7-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2178202-zero-padding-position-for-fft#answer_1567469

编辑：Matt J 2025-7-2

在 MATLAB Online 中打开

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

2 个评论
显示无隐藏无

Matt J 2025-7-2

在 MATLAB Online 中打开

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 2025-7-2

在 MATLAB Online 中打开

Or just

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

so as to not corrupt the original data vector.

请先登录，再进行评论。

Zero-Padding Position for FFT

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

7 个评论
显示 5更早的评论隐藏 5更早的评论

更多回答（1 个）

2 个评论
显示无隐藏无

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Zero-Padding Position for FFT

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

7 个评论 显示 5更早的评论隐藏 5更早的评论

更多回答（1 个）

2 个评论 显示 无隐藏 无

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

7 个评论
显示 5更早的评论隐藏 5更早的评论

2 个评论
显示无隐藏无