Effect of zero padding on FFT amplitude

조회 수: 220 (최근 30일)
Rahul Kale
Rahul Kale 2023년 1월 25일
편집: Paul 2023년 2월 11일
Many text books and other literature comment on the improvement in FFT resolution due to zero padding but I have not come across a text that comments on the effect of zero padding on FFT amplitude due to change in signal length (effectively the same energy is spread over longer time). In the real world, say for processing of vibration signal, engineers are interested in amplitude at a frequency.
Would zero padding not change original signal and thereby underestimate amplitudes of FFT?
Thanks
Rahul
  댓글 수: 10
Paul
Paul 2023년 1월 28일
1. Yes, x1 in the first case are the first six elements in x3 in the second case. I found a typo in that coment and fixed it, it now says: "x3 = [x2 zeros(1,Nfft-3)]." Maybe that typo confused things. I don't understand the second question here that starts "Or, are you saying ...."
2. I've reread this question many times, and not sure I understand it. We've been saying and showing in this thread that zero-padding doesn't affect the amplitude spectrum of the signal. As I've said and MattJ showed, zero-padding changes the location of the frequency domain samples of the DTFT, but the amplitudes of those samples still follow the amplitude of the DTFT the finite duration signal.
3. Yes, the first six elements of x3 are the same as the elements of x1. As for scaling the results of the DFT, what is the definition of "amplitude of the physical quantity" in either the first or second case?
Rahul Kale
Rahul Kale 2023년 1월 30일
Paul,
As for scaling the results of the DFT, what is the definition of "amplitude of the physical quantity" in either the first or second case?
- It is the 0-pk magnitude that a signal producing mechanism would produce. e.g. a shaker vibrating a structure.
Thanks
Rahul

댓글을 달려면 로그인하십시오.

채택된 답변

Matt J
Matt J 2023년 1월 28일
편집: Matt J 2023년 1월 28일
@Rahul Kale It occurs to me now that you are not really trying to compute the Fourier transform of a finite signal. You are instead trying to compute the Fourier series coefficients of a periodic signal, where x(n), n=0,...,N-1 is a sampling of one period of the signal. Note that the Fourier transform is not the usual way to decompose an infinite periodic signal, because it has infinite energy, although there are certainly mathematical relationships between the continuous Fourier transform a continuous Fourier series.
To get the discrete Fourier series (DFS) coefficients of a signal, then the way to do that using Matlab's fft() command is
DFS=fft(x)/N
So, 1/N is probably the normalization factor you are looking for. Zero-padding x will indeed change the Fourier series coefficients. However, it is important to understand that padded zeros have a different interpretation when applying the FFT to non-periodic signal analysis than to periodic signal analysis.
  댓글 수: 36
Rahul Kale
Rahul Kale 2023년 2월 3일
How to do that? I would certainly want to give him credit. Rahul
Bruno Luong
Bruno Luong 2023년 2월 3일
@Rahul Kale please don't bother. Let's move on.

댓글을 달려면 로그인하십시오.

추가 답변 (3개)

Bruno Luong
Bruno Luong 2023년 1월 29일
편집: Bruno Luong 2023년 1월 29일
If your (windowed) signal dies out to 0 in a smooth way then the padding with 0s on the tail changes little to the FFT spectrum (amplitude), it is equivalent to densify the spectrum (interpolation) without scaling it.
N=70;
x=(0:N-1)-(N-1)/2;
E=exp(-(x./(N/6)).^2);
y = E.*randn(size(x));
f = (0:N-1)/N;
yhat = fft(y);
figure
plot(f,y,'b',f,E,'c');
NPad = 2^nextpow2(N);
ypad=y; ypad(NPad)=0;
fpad = (0:NPad-1)/NPad;
ypadhat = fft(ypad);
figure
plot(f, abs(yhat), 'b-+', fpad, abs(ypadhat), 'r-.')
  댓글 수: 5
Rahul Kale
Rahul Kale 2023년 1월 31일
Bruno,
Thanks much for this code..I am still trying to understand each step...will let you know if there are more questions.
I am not sure why in the scale factor there is N and not Npad.
Thanks
Rahul
Bruno Luong
Bruno Luong 2023년 1월 31일
편집: Bruno Luong 2023년 1월 31일
@Rahul Kale "I am not sure why in the scale factor there is N and not Npad."
This has been somewhat explained in my various posts if you follow closely.
What happens here is the Parseval equality, or equivalent in discrete is the so called Rayleigh Energy Theorem and DFT that defined with the factor 1/sqrt(N) so make the transformation symetric in time and Fourier domain on the energy point of view.
Combine with the second fact that 0-padding does not change the amplitude of the FFT as I show in my answer, but just make the frequency vector with finer step (denser) on the same frequency domain [0,Niquist], there is no reason to adjust the amplitude with the length of the zero-pad signal.

댓글을 달려면 로그인하십시오.


Paul
Paul 2023년 2월 5일
편집: Paul 2023년 2월 11일
This thread started with one question about zero-padding the Discrete Fourier Transform (DFT) and the answers and comments brought up many points on other related concepts including the Continuous Fourier Series (CFS), Continuous Time Fourier Transform (CTFT), windowing, and possibly many others. As might be expected for a thread with so many comments and additional questions in the comments, it might be hard to follow how all of the concepts relate and flow to each other. Consequently, I thought it might be good to show exactly that, at least for some of these of these concepts.
I'm posting this as its own answer because I wasn't sure it would flow nicely as a comment in any of the other answers.
Also, this answer does not add to any of the material already posted. Again, it's only trying to use Matlab to illustrate some of the points already made by (in no particluar order) @Mathieu NOE, @Star Strider, @Bjorn Gustavsson, @Matt J and @Bruno Luong and me.
Define a simple, continuous-time signal
syms t f real
Pi = sym(pi);
T1 = sym(2);
T2 = sym(3);
x(t) = 12*sin(2*Pi/T1*t) + 4*cos(2*Pi/T2*t);
x(t) is the sum of two periodic signals. Because the ratio T1/T2 is a rational number, we know that x(t) is periodic. By inspection, its fundamental period is
T = sym(6)
which we can verify by
simplify(x(t) - x(t+T))
ans = 
0
Because x(t) is periodic, its CTFT is a train of scaled Dirac delta functions at the fundamental frequency (1/6 Hz) and its harmonics. In this case, we only get non-zero harmonics at f = +-1/3 and f = +-1/2
X(f) = expand(fourier(x(t),t,2*Pi*f))
X(f) = 
Because x(t) is periodic, it has an exponential CFS represenation. The CFS coefficients as a function of n are computed by
syms n integer
C(n) = int(x(t)*exp(-1j*2*Pi*n/T*t),t,0,T)/T
C(n) = 
Or, we could have used fourier to compute them just as well.
% C(f) = fourier(x(t)*rectangularPulse(0,T,t),t,2*Pi*f)/T;
% C(n) = C(n/T);
Because the integrand in int has n as a parameter, it doesn't capture the special cases where abs(n) = 2 or abs(n) =3. Clean that up and C(n) is pretty simple
[~,denC] = numden(C(n));
npoles = solve(denC);
for npole = npoles(:).'
C(n) = piecewise(n == npole,limit(C(n),n,npole),C(n));
end
C(n) = simplify(C(n))
C(n) = 
The CFS coefficients are exactly the cofficients in the impulse train that defines X(f), a fact that will be taken advantage of below.
Plot the amplitudes of the four CFS coefficients.
nval = -5:5;
figure
stem(nval/T,abs(C(nval)))
ylim([0 7])
xlabel('Hz')
Because we are using the exponential CFS, the amplitudes of C(n) and C(-n) are each 1/2 of the amplitude of the corresponding harmonic. Therefore, the plot tells us that the amplitude of the harmonic at 1/2 Hz is 12 and at 1/3 Hz is 4, which is exactly what we expect based on the definition of x(t).
Suppose we only have available to us a portion of x(t). The next step shows how Fourier analysis can still extract the amplitude information from the signal.
Define a window function that covers the portion of the signal that we have for processing. The triangular window function (for illustration) is parameterized by the number of periods of x(t) that it covers.
syms numperiod
w(t,numperiod) = triangularPulse(0,numperiod*T,t);
Suppose we have 10 periods of x(t)
nperiod = 10;
The CTFT of the window signal is
W(f) = simplify(fourier(w(t,nperiod),t,2*Pi*f));
Because w(t) covers a significant number of periods, its bandwidth in the frequency domain looks small relative to the frequency spacing between harmonics of x(t)
figure
fplot(abs(W(f)),[-0.3 0.3]);
xline( double( 1/2/T),'r', '1/(2T)');
xline( double(-1/2/T),'r','-1/(2T)');
xlabel('Hz')
ylim([0 31])
The peak of W(f) at f = 0 is
W0 = limit(W(f),f,0)
W0 = 
30
which is also area under the window.
int(w(t,nperiod),t,0,nperiod*T)
ans = 
30
Apply the window to x(t) and compute the CTFT of the result
xw(t) = x(t)*w(t,nperiod);
XW(f) = fourier(xw(t),t,2*Pi*f);
Multiplication in the time domain corresponds to convolution in the frequency domain. Convolution of W(f) with the Dirac delta functions in X(f) correponds to shifting W(f) to the corresponding harmonics and multiplying by the corresponding Dirac coefficients. But those Dirac coefficients are just the CFS coefficients. Therefore, we can just as easily express the CTFT of xw(t) directly as
nval = [-2 2 -3 3];
XWalt(f) = sum(C(nval).*W(f-nval/T));
Show that XW(f) and XWalt(f) are equivalent.
figure
hold on
% fplot was taking too long.
%fplot(abs([XW(f) XWalt(f)]))
fplot(matlabFunction(abs(XW(f))),[-1 1]);
fplot(matlabFunction(abs(XWalt(f))),[-1 1]);
xlim([-1 1]*3/4)
xlabel('Hz')
Because the bandwidth of W(f) is small relative to the harmonic spacing, XW(f) is essentially composed of replicas of W(f) shifted to the harmonic frequencies and scaled by the CFS coefficients. This result suggests that we can recover the amplitudes of C(n) at the harmonic frequencies by simply dividing XW(f) by W0.
figure
fplot(abs(XW(f)/W0))
hold on
stem(nval/T,abs(C(nval)))
xlim([-1 1]*3/4)
ylim([0 7])
xlabel('Hz')
If all we have is the blue curve, we would surmise that the underlying periodic signal is composed of two sinusoids at 1/2 Hz and 1/3 Hz with amplitudes 12 and 4 respectively.
With that framework in continuous-time, now move to discrete-time.
Define functions to evaluate numerically (not stricly necessary, just runs a bit faster than evaluating in symbolic and converting)
xfunc = matlabFunction(x(t));
wfunc = matlabFunction(w(t,numperiod),'Vars',[t numperiod]);
Assume a sampling period of
Ts = 0.01;
The number of samples to cover that many periods is
N = nperiod*double(T)/Ts;
and the associated time samples are
tval = (0:(N-1))*Ts;
Evaluate the x(t) and w(t) at the sample times.
xval = xfunc(tval);
wval = wfunc(tval,nperiod);
Compute the DFT
XWdft= fftshift(fft(xval.*wval));
nfft = numel(XWdft);
fval = ceil(linspace(-nfft/2,nfft/2-1,nfft))/nfft/Ts;
Now, to plot the DFT with the appropriate scaling, we need to divide out the effect of the window function. In discrete-time, without any other scaling on the output of fft, that corresoponds to dividing by the sum of the window samples
figure
stem(fval,abs(XWdft)/sum(wval));
xlim([-0.75 0.75])
xlabel('Hz')
If we had used a rectangular window, then sum(wval) = N and we'd just divide by N.
Overlay with the scaled version of XW(f)
hold on
fplot(abs(XW(f)/W0),'g')
The scaled DFT samples are essentially frequency domain samples of the XW(f)/W0 (in general, this relationship is only approximate where the approximation gets better as Ts gets smaller and N gets larger).
Zoom in to see things a little better
copyobj(gca,figure);
xlim([0.2 0.6])
Because the DFT samples effectively recover XW(f)/W0, we again know the harmonic frequencies and the corresponding amplitudes.
Four DFT samples lay exactly on the harmonic frequencies because of how N was computed to cover exactly 10 periods for the given the value of Ts. In general, we cannot guarantee the data covers an integer number of periods; after all, the period is what we are trying to determine.
Also, the selection of Ts in this example is such that x[n] = x(n*Ts) is periodic in discrete-time.
Suppose we have bit more than 10 periods of samples.
nperiod = 10 + 0.25;
N = round(nperiod*double(T)/Ts);
Going through the same process
tval = (0:(N-1))*Ts;
xval = xfunc(tval);
wval = wfunc(tval,nperiod);
XWdft= fftshift(fft(xval.*wval));
nfft = numel(XWdft);
fval = ceil(linspace(-nfft/2,nfft/2-1,nfft))/nfft/Ts;
figure
stem(fval,abs(XWdft)/sum(wval));
xline(.5,'r')
xline(1/3,'r')
xlabel('Hz')
xlim([0.2 0.6])
ylim([0 7])
shows that the DFT samples are not so "nice," particularly around 1/3 Hz, which is basically halfway between two DFT samples.
Now, we finally get to the original question about zero padding, which we can use to fill in DFT samples in the frequency domain
N = round(nperiod*double(T)/Ts);
tval = (0:(N-1))*Ts;
xval = xfunc(tval);
wval = wfunc(tval,nperiod);
XWdft= fftshift(fft(xval.*wval,2^14)); % use a lot of padding
nfft = numel(XWdft);
fval = ceil(linspace(-nfft/2,nfft/2-1,nfft))/nfft/Ts;
figure
stem(fval,abs(XWdft)/sum(wval));
xline(.5,'r')
xline(1/3,'r')
xlabel('Hz')
xlim([0.2 0.6])
ylim([0 7])
As has been stated repeatedly, zero padding has no impact on the unscaled amplitude spectrum, and we recover the correct scaled amplitude spectrum as long as we divide by sum(wval), not the number of samples in DFT (2^14 in this instance). The zero padding recovers the information we seek, which is that the harmonics are somewhere very near 1/2 Hz and 1/3 Hz, with amplitudes nearly 12 and 4 respectively.
Although it looks like the DFT samples are samples of XW(f), that's not true because XW(f) is based on a 60 second window and the DFT sampes are based on a 60 + 6/4 = 61.5 second window, which has a slightly different CTFT. However, as long as the peak of CTFT of the window (whatever it may be) is at f = 0, the same process applies, i.e., using the peaks of XWdft(n)/sum(wval) to estimate the harmonic frequencies and amplitudes, as long as the DFT is normalized based on the actual window used, as done in this example.

Matt J
Matt J 2023년 1월 26일
편집: Matt J 2023년 1월 27일
Here is an example showing that the amplitude of an FFT does not change due to zero-padding. In all cases, the peak amplitude is 5.
x=[1,1,1,1 1];
N=(1:5)*5;
for n=N
plotFcn(x,n); %zero-pad x by 2*n and plot FFT
end
hold off; legend("n="+N)
function plotFcn(x,N)
z=zeros(1,N);
x=[z,x,z]; %zero pad
M=floor(numel(x)/2);
X=abs(fftshift(fft(x)));
plot((-M:M)/numel(X), X,'o'); hold on
end
  댓글 수: 29
Paul
Paul 2023년 1월 28일
Another way to say this is that if c[k] are the exponential CFS coefficients of a periodic signal x(t) with period T, then the CTFT of that periodic signal is
X(w) = 2*pi*sum( c[k]*delta(w - k*w0) ) (w here in rad/sec, non-unitary transform)
where w0 = 2*pi/T and the summation is taken over all integers k.
It's been a good discussion, IMO.
Bjorn Gustavsson
Bjorn Gustavsson 2023년 1월 30일
By taking the weekend off I missed a lot of discussion. Since this is already a rather sprawling stream of fourier-transform-information I will not be ably to catch up and will stand by and read, but first a last comment on @Paul's reply:
"From my perspective, we are not necessarily taking the Fourier transform of a constant-valued function. Matt J has a finite sequence of data: [1, 1, 1, 1, 1]. Before doing any analysis I think about what I'm assuming about the rest of the signal that is unknown. If these five points are a subset of a constant signal, that's one thing. But if these five points are a subset of a rectangular pulse, then that's a different thing. There may be other things that these five points represent. Whatever that is, it guides me to how I would do the analysis, e.g., zero-padding or not, window or not, etc."
That is a somewhat valid point-of-view, but from an Occhamian (Akaike/Bayesian model-selection) view it will be very hard to justifying additional fourier-components outside of the DC one, since that is enough to perfectly explain the available data, extending outside of that is very dodgy.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Spectral Measurements에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by