Effect of zero padding on FFT amplitude

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

[Mathieu NOE]

hello

1 / i am also a noise & vibration engineer and I very rarely need to zero pad my data (I always take much longer records compared to what I really need in terms of resolution. I use the longer records to do linear averaging of fft to improve my SNR rather than the resolution )

2/ no there is no loss of amplitude by zero padding. yes you increase your fft resolution but you don't add any new information to the fft output.

some readings that may still interest you :

roboticist LEE, TAEWOO :: MATLAB - FFT and Zero Padding (ru.nl)

Does zero-padding a FFT matter? — Matt Gaidica, Ph.D.

There is still one case where I find a practical use of zero padding, is when you do modal analysis on very lightly damped structures. the best test signal is random burst with enough zero padding to let the structure response decay to zero. You can do multiple averages but always let the response decay before next record (average).

[Rahul Kale]

Mathieu,

Thanks for your reply. However, by adding zeros, we are actually changing the original signal. I think the Matlab example that you showed has scaled(normalized) it by the length and thereby the amplitude remained same.

Thinking from the viewpoint of a physical entity, say excitation force , if FFT amplitude remained the same for arbitrary zero padding then the duration between the two square pulses (for example) would not matter. However, structural response for the two excitations with different duration between pulses (and not the pulse duration itself) would differ and would be dependent on the time duration between the pulses.

How to address this apparent contradiction?

Thanks

Rahul

[Star Strider]

Zero-padding using:

NFFT = 2^nextpow2(L)

where ‘L’ is the length (or in a matrix, row size) of the signal array has the added advantage (in addition to increasing the frequency resolution) of making the fft calculation much more efficient. This is due to the way the Fast Fourier Transform is calculated. I also usually window the fft. This corrects for the fft being finite rather than infinite, and (in my opinion) improves the result.

.

[Bjorn Gustavsson]

@Star Strider: If I've understood things right the efficiency-gain with padding to the next power of 2 is less important with fftw - it is rather efficient as long as the size of the input has "few large prime-factors"?

Windowing is important, whether zero-padding or not. And to my understanding it is a bit imprecise to claim that no information is added by zero-padding - the discontinuity (in some derivative) at the transition between the and the zeros, does affect the transform.

[Paul]

@Rahul Kale

Can you clarify the concern about a "contradiction?" An actual example would be better.

The DFT output of the fft is always a set of frequency domain samples of the DTFT of the signal. The effect of zero-padding can only change where in the frequency domain the DTFT is sampled. So zero-padding can change how the DFT samples appear on a plot, but they will always be frequency-domain samples of the same DTFT.

[Rahul Kale]

Well, by contradiction, I meant that FFT of zero padded signal is the same as the original one except adding values at points on the X axis where they were not present( in the FFT of original signal, thus essentially improving the resolution). However, when we think of its physical interpretation, it means changing the signal in time domain by zero padding would not affect its frequency domain values but would cause different response.

To illustrate this, if we consider a bell (old fashioned one having very less damping and not electronic one) making sound for a strike( it is hit with a hammer) say every 2 sec and compare its sound with a bell that is struck one time and then again struck one time after say every 10-15 min. There is definitely different sound response that we would feel. However, we can consider the two types of strikes the same because from time 2sec to 10 min is equivalent to zero padding the first type of strike to get the second type of strike

I am not sure if I have conveyed properly what I wanted to ask.

Thanks

Rahul

[Paul]

Signals go on forever, both into the past, and into the future. However, the DFT is only defined for signals x[n] with finite duration of length N, where the signal is 0 for n < 0 and the signal is 0 for n >= N. It's important to understand that this condition does not preclude a run of zeros for n <= N-1.

One way to interpret the DFT is that the results are always the Discrete Fourier Series coefficients of the periodic extension of x[n] with period Nfft >= N. If Nfft = N, then we are computing the DFT w/o zero-padding and if Nfft > N then we are using zero-padding.

Assuming for this discussion the we are limited to deterministic signals, we can use the DFT in one of two ways.

First, x[n] is one period of a periodic signal with period N, and it is this periodic signal we wish to characterize. In this case, we set Nfft = N and compute the DFT. We now have the DFS coefficients of the periodic signal that we actually care about. Suppose we strike the bell three for three beats, starting from n = 0, and then rest for three beats, and we do this for all time into the past and future. So one period of our striking sequence would look like x1 = [1 1 1 0 0 0], with N=6. Those last three zeros are not zero-padding. They are very important in defining one period of the signal. We can take the DFT of that sequence, without zero-padding, to obtain its DFS coefficients, from which we could compute x[n] for any value of n. These six DFS coefficients are all we need.

Second, x[n] is actually a finite duration signal. In this example, we don't strike the bell until n = 0, then srike it twice more, then stop. So x[n] looks something like this x[n] = [... 0 0 0 0 0 0 1 1 1 0 0 0 ...] where the ellipses indicate never ending and the first 1-entry correponds to n = 0. Now, this signal has frequency content at all frequencies, as could be seen by computing its DTFT. However, we want to use the DFT to compute samples of one period of the DTFT. The way we do that is to start with the sequence x2 =[1 1 1] but take Nfft large, i.e., lots of zero padding, so the DFT of x3 = [x2 zeros(1,Nfft-3)]. Now, when we take the DFT of x3, we are getting the DFS coefficients of the periodic extension of x3, with period Nfft, and these are still frequency domain samples of one period of the DTFT of x[n]. As Nfft >> N, then the periodic extension of x3 becomes closer and closer to actually being x[n] (imagine what happens as Nfft -> inf), and the DFT samples become more and more dense giving us a better picture of one period of the DTFT of x[n].

[Rahul Kale]

Paul,

Thanks much for such a detailed explanation. However, I have some followup questions. Hope you would not mind. I apologize for dragging this Post longer.

1) In your first example, the statement,"Those last three zeros are not zero-padding" is very important. So the question is, isn't this(x1) the subset of the seriesx[n] in your second example(x3)? Or, are you trying to say that, since the Nfft>>N in the second case, the zeros(that are not due to zero padding in the first case) do not matter?

2) My confusion is exactly this that how to determinethe period of the signal for FFT scaling when we zero-pad a signal. Whether we need to take the period of the whole signal (original + padded zeros) or just that of the original signal. If its former then wouldn't the FFT become subjective to who does it with how many zeros because the scaling factor and thus the amplitude of FFT depend on the period of the original signal?

3)"Now, when we take the DFT of x3, we are getting the DFS coefficients of the periodic extension of x3, with period Nfft, and these are still frequency domain samples of one period of the DTFT of x[n]"

- Here, x3 is comprized of x1 as the original signal. So how would you scale FFT of x3 if we are interested in absolute amplitude of the physical quantity(say force) of x3 as a funciton of frequency? How much would be the signal length that you would use for dividing(6 which is of x1 or 3 which is of x2)?

Thanks

Rahul

[Paul]

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]

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

Answer 1

[Matt J]

@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.

[Matt J]

Because different signal processing texts define transforms slightly differently, I will add here that the above comes from a discretization of the formula for the k-th exponential Fourier series harmonic as defined here.

where T is the period of x(t). We discretize this integral into a sum by writing dt=T/N, t=n*T/N, and =x(n*T/N) leading to the approximation,

[Rahul Kale]

Matt,

Yes, I am trying to do an FFT of a periodic infinite wave by taking its 1 period length of that signal. Your statement,

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

is very important. Can you please elaborate on this?

What I ultimately want to do is given below,

1) I have two signals(that can be considered periodic) with different length in time, different sampling interval and probably different frequency content.

2) I want to compare the two signals w.r.t their amplitudes at certain frequencies, so, doing FFT is imperative.

Now that I have two signals that have diffrent length, I do not want to have their lengths affecting their FFT amplitudes. So I am going to zero-pad those so that they are of the same length and then do the FFT operation on each one. Also, I guess scaling(1/N) would be required for comparison.

Thanks

Rahul

[Bruno Luong]

When you pad 0s to a signal you change the amplitude of ALL frequencies. So it is a bad idead if your goal is to compare the amplitude of certain frequency.

[Rahul Kale]

A small correction. I do not take 1 period of two signals for FFT. Because, I do not know how many frequencies are in each of them. I take some time length of each signal and append zeros. I take x1[N1] + <zeros_1> and x2[N2] + [zeros_2] such that N1 + length(<zeros_1>) = N2 + length<zeros_2>

Thanks

Rahul

[Bruno Luong]

I would rather take the FFT of the signals (without padding 0s, which IMO is horrible modification for your goal), probably windowing them to avoid artefact of truncation, then interpolate the spectra in order to compare them.

[Rahul Kale]

Bruno,

Yes, I am using window but after padding (maybe that is effective for the beginning of the signal only, in this case, because padding makes it already zero at one end).

The reason for padding is I want to take both the signals to the length NFFT = 2^nextpow2(L); So I am aim to find a common NFFT for both the signals by padding.

If what I am doing makes sense then I would like to know what should be the normalizing factor after padding, whether the original signal length or the padded signal length. I am also awaiting Matt's answer to my question in my previous comment.

Thanks

Rahul

[Matt J]

"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." is very important. Can you please elaborate on this?

Well, imagine for a moment that we have a continuous time signal and that we are going to sample it, but the time sampling rate is fixed and given at (the exact value won't matter) and we have no control over this.

If the continuous signal is periodic, then we are going to take its Fourier Series. To do that, we sample one period of the signal at intervals to obtain a sequence and use the discretized formula I derived above for the Fourier series coefficients,

In this case, because the signal and its time sampling rate are fixed and given, the length N of the sequence is not under our control. The derivation of c_k required that we integrate over a complete period, so with fixed, the number of samples N over the period is fixed as well.

On the other hand, suppose instead that the signal is of finite duration. So, we are not going to compute the Fourier series, but instead the Fourier transform which is this integral:

I am now going to discretize this integral, again under the assumption that second is fixed and given. However, in this case, because the signal is finite duration, I can replace the infinite integral limits with a finite limit D/2 so long as D is large enough to cover the duration of the signal,

I discretize the integral with t=n , dt=, f=k/D, and N=D/ leading to

There are two differences to notice between the periodic case and the finite duration case. First, the normalization of the FFT is different when computing as compared to . Second, in the finite duration case, the parameter D and therefore N=D/ is under our control even though the sampling rate is fixed. Increasing or decreasing D results in more or less padded zeros in . It also has the affect of changing the frequency sampling interval =1/D, but not the underlying continuous spectrum, as has been demonstrated elsewhere.

Now let's go back and relax the assumption that is fixed. In the case of the periodic sginal, since we now control , we now also control N in the calculation of the Fourier series coefficients . However, the proportion of samples which are pulse samples and the proportion which are padded zeros remain the same, assuming the period of the given continuous signal has not changed. Conversely, in the case of the finite duration signal, both the truncation interval D and are freely variable, so the proportion of pulse samples and padded zeros can be controlled arbitrarily.

[Paul]

@Rahul Kale

"1) I have two signals(that can be considered periodic) with different length in time, different sampling interval and probably different frequency content."

For clarification, the "two signals" are continuous time signals, correct?

What does "considered" mean in this context?

Regarding "different length in time", this means the window of time over which the samples have been collected is different for each signal?

"2) I want to compare the two signals w.r.t their amplitudes at certain frequencies, so, doing FFT is imperative."

If two continuous-time signals x1(t) and x2(t) are periodic with periods T1 and T2, the only frequency content in each are at the CFS frequencies k/T1 and k/T2 (Hz) respectively. Obviously, if T1 = T2, then we can compare at all the CFS frequencies. If T1 ~= T2, then we can only compare, at best, at the subset of k1,k2 where k1/T1 == k2/T2. In light of this, what are the "certain frequencies" for the comparison?

[Rahul Kale]

Paul,

1)The two signals have been obtained from a sensor during testing and they are digitally sampled. By "considered" I mean, I am going to process them assuming that they are periodic because the sources generating them are supposed to produce periodic signals in theory.

2) "different length in time", this means the window of time over which the samples have been collected is different for each signal? - Yes

3)"If two continuous-time signals x1(t) and x2(t) are periodic with periods T1 and T2, the only frequency content in each are ...." - There is not just 1 frequency present in the signals. They would have many frequencies.

4)"certain frequencies" - These are frequencies that are expected to have significant amplitudes

Please let me know if more info is needed.

Thanks

Rahul

[Matt J]

By "physically correct" I mean, if, let's say, a person weighing 50kg is squatting sinusoidally over a board, then the FFT of the force measured by a load cell attached between his feet and the board should show a peak corresponding to his mass of 50kg.

That is more indication that you are trying to compute the discrete Fourier series. The DFS coefficients will give the amplitudes of sinusoidal components of the input,

N=100;
t=linspace(0,1,N+1); t(end)=[]; 

force=50*sin(2*pi*t);

X=fft(force,N)/N; %Discrete Fourier Series 
fAxis=(0:N-1);

stem(fAxis,2*abs(X)); xlim([0,10]); xlabel ('Frequency (Hz)')

The factor of 2 in 2*abs(X) is due to the fact that a real sinusoid of frequency f0 will be decomposed into 2 complex sinusoids of frequencies +f0 and -f0.

[Bruno Luong]

"There is not just 1 frequency present in the signals. They would have many frequencies."

Any period signal can have infinity frequencies (the Fourier series has ininity term), all are harmonics of the base frequency f1=1/T.

Paul speaks about mathematical periodic signal, if you have data from sensor, they are NOT periodic signal due to transient regime, incertitude of measurement, non perfect excitation, noise ... SO forget about consideration as (mathematical) periodic. It just confuse people who try to answer you.

As I see, the underline of the FFT is the signal is cycling of period N, where N is the sample length, or your signal is periodic by repeating the original sequence. If you do 0 padding, then actually you don't care if your signal periodic or not. So just leave it out the discussion.

[Rahul Kale]

Bruno,

I have not yet gone through your code. Will go through later.

My question to Matt (on his explanation on scaling) is, in the DFS example you have given above, what should be the scaling factor N in the case when 1) Original signal of length N is used 2) original + zero padded signal is used. Should N be always the length of original signal irrespective of whether I zero-pad or not? or, should I add zero padded length to N and then divide by the total length?

Bruno, you have advised me to integrate FFT to get the correct amplitude. However, it seems Matt has gotten amplitude right without integration. Am I missing something?

THanks

Rahul

[Matt J]

@Rahul Kale

Since your signal is periodic (you've been saying it is), then there shouldn't be any distinction between the "original signal" and the "zero-padded" signal. If your continuous-time signal consists of a train of pulses repeating every 10 seconds, then your discrete vector should be constructed from taking samples at over exactly a 10 second interval and the length of should be exactly N=10/. Any zeros present in will therefore be there because they were there in the continuous signal as well, and not because you added them after the fact.

[Bruno Luong]

"Bruno, you have advised me to integrate FFT to get the correct amplitude. However, it seems Matt has gotten amplitude right without integration. Am I missing something?"

Sorry I don't know in detail what is Matt logic.

[Rahul Kale]

Matt,

As said before I want to compare amplitudes of two signals with different sampling rate and total time, at each frequency. Although the signals contain many frequencies, I do not exactly know beforehand what the frequencies are in the two signals so I do not exactly know what total time should be truncating the two signals at and the deltaT I should be using. Those are decided based on prior theoretical knowledge of the excitation producing mechanism. So whatever time I truncate those at, I would introduce zero-padding if I need to select NFFT = 2^nextpow2(L);

In this case, the scaling factor would inherently include the length of padded zeros as well!

This means, based on your comment above, amplitude at each frequency would be scaled differently due to the zero padding length that each frequency component would see depending on how many integer periods fit in the selected total time (keeping aside effect of windowing for the sake of this discussion with zero padding). Is my understanding correct?

Thanks

Rahul

[Matt J]

@Rahul Kale It sounds like you haven't committed to a decision on whether you are doing periodic signal analysis or non-periodic, finite duration signal analysis. This decision must be made first before anything else, and it must be made in continuous time. Once you know the kind of continuous-time analysis you are trying to approximate, decisions can be made on how best to discretize the analysis using FFTs

Although the signals contain many frequencies, I do not exactly know beforehand what the frequencies are in the two signals

If the signals are periodic, as you have been saying, then you should know. A continuous signal of period T contains a discrete set of harmonic frequencies 1/T, 2/T, 3/T,.... If the two signals have two different periods T1 and T2, then their harmonic frequencies will not coincide and there is no way you can compare them.

This means, based on your comment above, amplitude at each frequency would be scaled differently due to the zero padding length that each frequency component would see depending on how many integer periods fit in the selected total time (keeping aside effect of windowing for the sake of this discussion with zero padding).

I made many "comments above" at this point, so I don't know which one you are referencing. If you are doing periodic signal analysis, as you have been saying, then there is no reason to be windowing the signal down to finite length. The signals should be viewed as infinitely repeating.

If you are windowing your signal because you do in fact want to do finite-duration signal analysis, then in continuous-time you need to be using the Fourier transform to analyse the signals, not the Fourier series. And, if you are to be using Fourier transforms, all of the expectations you've expressed about the amplitude of the spectrum are invalid. In particular, if you take the Fourier transform of 50*sin(2*pi*t), it should not be expected to have a peak of 50.

[Rahul Kale]

Matt, Sorry if I have confused you. Firstly, the two signals will have same frequency content. Secondly, signals are periodic but samples ought to be finite in time because they are the outcome of a sensor used in the test. My total time can be decided based on lowest frequency (1 period of it). But I can't guarantee that all frequency components will fit integer number of periods in the time of the sample collected because there are other non harmonic frequencies present. Hope this clarifies. Please let me know if more explanation is needed Thanks Rahul

[Rahul Kale]

I repeat the following for clarity. Each signal has many frequencies but the two signals contain same set of frequencies. The two data samples are obtained from a test using sensors so the samples are of finite length containing periodic quantities and my interest is to compare amplitudes at each frequency contained in the signals. The two samples have different lengths and sampling rate.

[Matt J]

I think you are saying that signal #1 has period T1 and signal #2 has period T2<T1, but there are no integers m1 and m2 such that m1*T1=m2*T2?

If so, I think you will need to choose a window that truncates one of the signals mid-cycle. That will have an impact on the amplitude of its spectral components, but one which has nothing to do with zero-padding selection.

[Rahul Kale]

Yes correct..

Instead of truncating, I select NFFT = 2^nextpow2(L); This results in zero padding.

One question to your earlier comment, "If you are windowing your signal because you do in fact want to do finite-duration signal analysis, then in continuous-time you need to be using the Fourier transform to analyse the signals, not the Fourier series. And, if you are to be using Fourier transforms, all of the expectations you've expressed about the amplitude of the spectrum are invalid"

Any data processing, even of an infinite periodic wave, will have to be done on a finite data only, correct?

In the sine wave with 50 amplitude example, you have used finite data only.

What difference in the processing would it make if the the same two signals containing periodic components were to be processed 1) using infinite data length(hypothetically speaking) 2) using finite data length. Ultimately the end result should be the same i.e. correct amplitude...isn't it? Would fourier transform not give correct amplitude and fourier series give that instead? Or, are you suggesting that if I aim for getting correct amplitude, I should do data processing in such a way that I end up using fourier series and not transform?

Thanks

Rahul

[Bruno Luong]

1. The Fourier series is defined for periodic signal. It can be approximated by DFT analysic on one period of the signal. The more sample you have the better approximation will be.
2. The Continuous Fourier transform (CFT) is defined for ANY function that is square integrable, might or might have infinity length.
3. If the function is periodic of period T, the Fourier transform is NOT defined (the function is not square integrable)
4. However if we window this infinite periodic function, you can take CFT. Now if you take the window larger and larger such that it covers more and more periods, the CFT will converge to a Dirac-comb weighted by coefficients that are proportional to the Fourier series coefficient described in 1. The frequency of the first dirac is at the fundamental frequency f=1/T.
5. The shape of the CFT around the fundamental frequency and its harmonic (k*f with k integer) are the CFT of the envelop window function. When the window get wider, the amplitude of those spectrum peaks with go up and the width will shrinking down so that the integral remains (more or less) constant. Therefore the amplitude is hard to interpret. You select the window function by making a compromise between two things: (A) nice decay in frequency domain, but (B) do not cut down too much the original signal in time domain. NOTE: If you truncate abrutely your time signal, meaning using a rectangular window function, you get all kind of strange effect because the CTF of a rectangular is a sinc function that decay slowly to 0, which is not desired.
6. You can do zero-pad in the FFT during analysis of the windowed signal that will have as effect of increasing the resolution in the Fourier domain. The amplitude of the spectrum does NOT affected by zero padding.

[Matt J]

Any data processing, even of an infinite periodic wave, will have to be done on a finite data only, correct?

Yes, but as I showed you in earlier posts, the Fourier Series integral (for periodic signals) and the Fourier Transform integral (for non-periodic signals) are different, and lead to different discretizations.

What difference in the processing would it make if the the same two signals containing periodic components were to be processed 1) using infinite data length(hypothetically speaking) 2) using finite data length.

It depends how you are truncating the periodic signal to finite length. For the example signal 50*sin(2*pi*t), I was able to get the peak Fourier Series coefficient of 50 because I truncated to exactly one period. I would also have gotten the same peak value if I'd truncated to an integer multiple of the period. In any other case, though, I would not have gotten the same peak coefficient.

This is relevant to you because you have two signals of different periods and might have attempted truncating one of them to a fraction of a period in order to get both of them in a common time window.

Instead of truncating, I select NFFT = 2^nextpow2(L); This results in zero padding.

But do you start padding the shorter-period signal after one period, or several? One might suppose that you are letting the shorter signal run for as many periods as it can complete in t<L before zero-padding starts. An argument can be made for both padding schemes.

[Rahul Kale]

Matt and Bruno,

I really appreciate the patience with which both of you have continued answering my questions although some of them might be too trivial for you. But for me, your answers are helping a lot.

Matt, I make sure that a signal has sufficient number of periods included before zero padding starts. Time length of at least 2 periods of the lowest possible frequency is what I used. So, that would include many periods of higher frequencies. Zero padding is done to the signal as measured (which is the combination of all frequency components) and not to individual periods. What argument would you make for both padding schemes w.r.t getting amplitudes right?

Thanks Rahul

[Matt J]

Time length of at least 2 periods of the lowest possible frequency is what I used. So, that would include many periods of higher frequencies.

Your terminology is very hard to follow. I think you meant to say "Both signals were truncated to two periods of the slower-oscillating signal. So, that would include many periods of the faster-oscillating signal."

[Matt J]

What argument would you make for both padding schemes w.r.t getting amplitudes right?

We can see the differences in the example below. Signals x1 and x2 both have period T=1. However, T=1 is not the fundamental period of the faster-oscillating signal x2, but rather T2=1/2. If I take the Discrete Fourier Series truncating both signals to a horizon of T=1, I get the result below. The harmonic frquencies of both signals have the expected Fourier Series coefficient amplitudes.

N=100;
t=linspace(0,1,N+1); t(end)=[]; 
fAxis=(0:N-1);

close all

x1=100*(sin(2*pi*t) + sin(8*pi*t));  %Slower signal
x2= 60*(sin(4*pi*t) + sin(12*pi*t)); %Faster signal

X1=fft(x1,N)/N; %Discrete Fourier Series 
X2=fft(x2,N)/N; %Discrete Fourier Series 

tiledlayout(1,2)
nexttile; 
plot(t,x1,t,x2); 
xlabel 'Time'; legend('Slower','Faster');axis square

nexttile
stem(fAxis,abs(X1)); hold on
stem(fAxis,abs(X2)); hold on
xlabel ('Frequency (Hz)'); legend('Slower','Faster');axis square
xlim([0,15])

sgtitle('Multi-Period Truncation')

However, if I do the anaysis of x2 by letting it run only for 1 period (0.5 sec.), and zero-padding the rest, each signal's coefficients have the same amplitudes as before at their respective harmonics. However, at 1 Hz and 4 Hz (which are harmonics of x1 but not x2), the coefficients of x2 get sinc-extrapolated, resulting in a much more continuous spectrum. If continuity is important to the comparison you are doing, this might be the better sampling scheme.

figure;
x2(end/2:end)=0;
X2=fft(x2,N)/(N/2); %Discrete Fourier Series -interpolated

[Paul]

Referring to @Bruno Luong's comment

Considering items 3 and 4 together, it's true that the CTFT integral does not converge for a periodic signal, but we still do define a CTFT for a periodic signal as an impulse train spaced by the harmonic frequencies and scaled by exponential CFS coefficients of the signal. This definition can be arrived at by the limiting process as Bruno described, which I'd imagine is the rigorous approach. But it can be intuitively understood by taking that scaled impulse train back through the inverse CTFT integral, which yields the CFS representation of the periodic signal. Off topic, but worth noting, is that a CTFT can be defined for other, aperiodic signals for which the CTFT integral doesn't converge, such as

x(t) = heaviside(t) -> X(f) = dirac(f)/2 - 1i/(2*f*pi)

Considering item 5, it seems to me that as long as the bandwidth of the CTFT of the window is enough less than 1/T, and if the peak of the CTFT of the window is at f = 0, then the spectrum peaks of the CTFT of the product of the periodic signal and the window will be at the harmonic frequencies, and the amplitude of those peaks will the product of the amplitude of the CFS coefficients and the peak of the CTFT of the window, which is known. Basically, I'm trying to understand why "the amplitude [of those spectrum peaks] is hard to interpret," assuming that the window bandwidth is narrow relative to 1/T. (Note: I would have shown an example with code, but the RUN feature isn't working for me.)

[Bruno Luong]

Matt's example is a perfect illustration of the "note" of my point 5. Truncating the signal at half (effectively one) period, the two frequencies are spread out by a sinc function then blended together (convolution between 2 diracs and sinc).

Taking more periods (>=10), using better window function than the rectangle one, you'll get a nice spcectrum curves with two clearly separate peaks.

N=1000;
t=linspace(0,10,N+1); t(end)=[]; 

x2= 60*(sin(4*pi*t) + sin(12*pi*t)); %Faster signal
x=(0:N-1)-(N-1)/2;
E=exp(-(x./(N/6)).^2);
XW2=fft(E.*x2);

plot(abs(XW2(1:ceil((end+1)/2))))

Then you can do zero padding on top of it.

[Rahul Kale]

Matt,

Yes, you have understood my "terminology" correctly.

Now, here is one important and subtle observation that I would like to point out.

To your comment, "However, if I do the anaysis of x2 by letting it run only for 1 period (0.5 sec.), and zero-padding the rest, each signal's coefficients have the same amplitudes as before at their respective harmonics"- This is true only when you used scale factor as N/2 and not N, for fft.

In your example where you truncated x2 and zero padded it, you divided its fft by N/2 and not N as in earlier case when there was no truncation. This was exactly my initial question that to get right amplitudes, whether I should consider signal length including the length of zero paddings for scaling the spectrum. If you had scaled using N, in second case, the signal amplitude would be half of the actual.

In my case,

N1=(length of sensor output data)

N2= N1 + length of padded zeros due to making NFFT = 2^nextpow2(L);

So, whether I should be scaling by N1 or N2?

According to your example it appears to me that I should be using the actual length of signal even for FFT of zero padded signal and that is N1 in my case.

Do you agree?

Thanks

Rahul

[Matt J]

This was exactly my initial question that to get right amplitudes, whether I should consider signal length including the length of zero paddings for scaling the spectrum.

Well, then it appears I've answered your question. I hope you will Accept-click the answer. :-)

According to your example it appears to me that I should be using the actual length of signal even for FFT of zero padded signal and that is N1 in my case.

It really depends on your attitude toward the padded zeros and whether or not they're supposed to be considered part of the period of the wave. We don't always just ignore zeros. Suppose my signal was the double sawtooth below. The period of this wave is 6, and you would have to normalize by 1/6 in a proper Fourier series calculation. If you were tempted to say the period is 3 because the length of the non-zero part of the oscillation is 3, you would be wrong. If you were tempted to say the period is 4 because only the length between the start of the front-facing tooth and the back of the rear-facing tooth matters, you would also be wrong.

The reason the interval of zeros doesn't matter to you when you consider the earlier sinusoidal example is because you apparently don't consider them part of the period of the wave. The zeros there serve a different purpose. If you did consider them part of the wave, then you would be wrong to consider the "correct amplitudes" of the harmonic coefficients of x2 to be 30.

  f=@(t) max((abs(t)-0.5).*(abs(t)<=2),0);
  x=@(t) f(mod(t,6)-2);
  
  fplot(x,[-0.1,18]);  ylim([0,3])

[Rahul Kale]

Paul, Matt and Bruno

Thanks much for your help. I am accepting the answer. However, I do not see the "Accept" button.

Rahul

[Matt J]

@Rahul Kale You have accepted Bruno's answer, so it looks like you found the button. Was that the answer you meant to accept?

[Rahul Kale]

Actually, I do not remember when I pressed that button. Bruno's help and his code has added a lot of insight for me. However, I feel your last code showing multi-period fft is the direct answer to my question. So I would accept that as the answer.

Thanks

Rahul

[Rahul Kale]

Bruno,

I appreciate help from both of you, you and Matt. I would like to accept answers from both of you if that was possible. However, I think there is provision to accept only one answer and as said in my reply to Matt, I feel Matt's answer directly addresses my question. Hence accepting his answer.

I am sorry...and thanks a lot

Rahul

[Matt J]

@Rahul Kale You can upvote Bruno's answer, as I have, if you want to award credit to it. Also, he will not lose the points that were awarded to him from your Acceptance, even though you've now withdrawn it.

[Rahul Kale]

How to do that? I would certainly want to give him credit. Rahul

[Bruno Luong]

@Rahul Kale please don't bother. Let's move on.

Answer 2

[Bruno Luong]

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-.')

[Rahul Kale]

Bruno,

I did try to do resampling by interpolating to get NFFT = 2^nextpow2(L) within the signal length(original and no zero padding) but that is tedious and probably(I think) would introduce unknown frequencies in the FFT so I would rather zero-pad if that is acceptable by FFT procedure to produce meaningful amplitudes. Please note, I am also interested in getting physically correct amplitudes and not just comparative amplitudes with some scale factors. By "physically correct" I mean, if, let's say, a person weighing 50kg is squatting sinusoidally over a board, then the FFT of the force measured by a load cell attached between his feet and the board should show a peak corresponding to his mass of 50kg.

Thanks

Rahul

[Bruno Luong]

"I did try to do resampling by interpolating to get NFFT = 2^nextpow2(L)"

Why you do that? I show zero-padding IS equivalent to interpolating in frequency doman IF you signal decays to 0.

Interpolation frequency does NOT introduce new frequency (how on earth linear interpolation can?), whereas zeros padding with non decay signal will do all sort of effect of the spectrum.

If you are windowing you signal correctly (as you claim) then your signal will decay to 0, then just 0-pad as in my code show, it is more or less interpolation in spectral domain: the amplitude is unchanged, the frequency has to be adjusted with the new array indices.

Afterall that is the main purpose of windowing when study signal in Fourier domain: avoid artefact.

The question of how to scale the FFT to get "physically correct" is a separate question, regardless if you zero pad or not the signal. It seems Matt has address the normalization for your goal (though I haven't read his answer carefully).

The one sided spectrum of y

% yhat = fft(y)*2*sqrt(dt/(N-1))

has unit force/sqrt(Hz) if you original signal is Force, and associated with the corresponding frequency vector f. If you want to get back to force unit you must integrate the yhat^2 on a certain non-zero interval (the width is ~ bandwidth) and then take a square-root of the integration. There are probaby a tone of discussions on this topic. The bandwidth is non-zeros since your signal is finite in length, even for perfect input sine wave, since it is truncated.

% sample frequency and period
dt=0.1;
fs=1/dt;

% signal length and time vector
N=4000;
t=dt*(0:N-1);

% generate tested signal mono frequency
T=10;
A=2023;
y=A*sin(2*pi/T*t);

% Fourier transform and one-sided spectrum
df=fs/(N-1);
yhat=fft(y)*2*sqrt(dt/(N-1));
yhat=abs(yhat(1:ceil((end+1)/2)));
nf = length(yhat);
iseven = mod(N+1,2);
if iseven
    duplicatedside = union(1,nf);
else
    duplicatedside = 1;
end
yhat(duplicatedside) = yhat(duplicatedside)/sqrt(2);
% associate frequency vector
f=df*(0:nf-1);
plot(f, yhat);
xlabel('f [Hz]')
ylabel('spectrum [V/sqrt(Hz)]')
xline(f(end),'-','Nyquist','LabelHorizontalAlignment','left')

% Retrieve amplitude
fres=1/T;
bw=0.2;
fint = fres+bw/2*[-1 1];
keep = f>=fint(1) & f<=fint(2);
Amplitude = sqrt(trapz(df,yhat(keep).^2))

Amplitude = 2.0235e+03

xlim(fint)

The direct interpretation of the amplitude of the FFT is always problematic, only the energy and RMS value via L^2 integration in fourier domain via and parseval is mathematically correct.

[Bruno Luong]

Here is the full code that does windowing and then regroup both facts

• zeros padding => amplitude FFT does not change on windowing signal
• Correct scaling to get back amplittude after integration

%% 
%clear
%close all

% sample frequency and sample period
dt=0.3;
fs=1/dt;

% signal length and time vector
N = 600;
t = dt*(-(N-1)/2:(N-1)/2);

% generate tested signal mono frequency of period T
T=10;
A=2023;
y=A*sin(2*pi/T*t);

% Windows function
Ewidth = t(end)/2;
E = exp(-(t/Ewidth).^2);

% Windowed signal to make them decay to 0 smoothly
ywin = y.*E;

figure
h = plot(t, [y(:) ywin(:) A*E(:)]);
xlabel('t [s]')
ylabel('signal [V]')
legend(h, 'original','windowed', 'window');

% Fourier transform of 0 padded signal and one-sided spectrum
Npad = 8*2^nextpow2(N); % 8 factor for fine resolution in Fourier domain
df = fs/Npad;
yhat = fft(ywin,Npad);
scalefft = 2*sqrt(dt/(N-1)); % N, not NPad
yhat = scalefft*abs(yhat(1:ceil((end+1)/2)));
nf = length(yhat);
iseven = mod(Npad+1,2); % EDIT fix a small glitch here
% The DC and Nyquist are common for both sides
if iseven
    duplicatedside = union(1,nf);
else
    duplicatedside = 1;
end
yhat(duplicatedside) = yhat(duplicatedside)/sqrt(2);
% associate frequency vector
f = df*(0:nf-1);

figure
plot(f, yhat);
xlabel('f [Hz]')
ylabel('spectrum [V/sqrt(Hz)]')
xline(f(end),'-','Nyquist','LabelHorizontalAlignment','left')

% Retrieve amplitude from Parseval
fres=1/T; % we can also detect the frequency of the peak of yhat, but we know it's 1/T
bw=sqrt(2)/Ewidth*sqrt(2);
fint = fres+bw/2*[-1 1];
keep = f>=fint(1) & f<=fint(2);
RMSE = sqrt(trapz(E.^2)/(N-1));
Amplitude = sqrt(trapz(df,yhat(keep).^2)) / RMSE % Estimate of A

Amplitude = 2.0230e+03

% Zoom in the peak of the spectrum
xlim(fint)

[Rahul Kale]

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]

@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.

Answer 3

[Paul]

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.

Answer 4

[Matt J]

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

[Bjorn Gustavsson]

This is a wonderful, great and horibly uggly example of what happens with zero-padding!

The uggliness stands out ever more starkly if one labels the frequencies on the x-axis - DC-component has the right magnitude, and then there is a first harmonic component with close to the same amplitude...

[Paul]

@Bjorn Gustavsson

I can't tell if your comment is implying that zero-padding is good or bad (or perhaps something else?). The zero-padding is filling in the gaps in the frequency domain, as it should.

[Bjorn Gustavsson]

@Paul: This is true. But Matt J's example also shows that without windowing the fft of a constant (all 1s in x!) signal suddendly also has almost as high amplitudes at the first 2 non-zero-frequency components, and then at all higher frequencies you have some power - corresponding to a sinc - due to the two steps (0 to 1 and 1 to 0) in x - no windowing corresponds to windowing with a top-hat filter with the length of the signal. Zero-padding can be used, but one should know what it means, and what it leads to, and it ought to be combined with windowing of the signal.

[Rahul Kale]

Please see my comment in response to Paul's question.

In fact the example given by Matt J essentially shows that zero padded spectrum is different. So if X axis were frequency then it is different than the original signal. And, if this were a force acting on a structure, then the structure's frequency response would be different for the original and zero padded signals

Thanks

Rahul

[Bjorn Gustavsson]

@Rahul Kale, yes the spectrum is different. But you also need to think about why and how it is different and what that means. In your "force on structure" scenario the x-only signal would correspond to a static force of 1. The zeropadded signal would correspond to a force that suddendly appears with a magnitude of 1 at "time" 5, then dissapears at time 9. Whe you look at the fouriertransform you have the spectrum of a signal that is periodic with period-time of 13 (or 14, I will count all the samples) and alternates between zero and one, and that is obviously different from a static force.

[Paul]

@Bjorn Gustavsson

" Zero-padding can be used, but one should know what it means" - Agree 100%

"and what it leads to," - Agree 100%

"and it ought to be combined with windowing of the signal." Depends on the problem to be solved.

Suppose the problem statment is: Use the DFT to get frequency domain samples of the DTFT of a finite-duration sequence. In this case, using anything other than a rectangular window would be inappropriate, IMO. Of course, there may very well be, or actually are, other problem statement for which windowing is appropriate. Depends on the problem to be solved and the assumptions of the underlying signal.

[Matt J]

In fact the example given by Matt J essentially shows that zero padded spectrum is different. So if X axis were frequency then it is different than the original signal.

I have refined my example to better show that the underlying spectrum and its amplitude is not changing with zero padding. You're just seeing more samples of it.

[Paul]

Excellent update. Except now there may be several comments under this answer that were referring to the original example and might not be as meaningful relative to the update.

[Bjorn Gustavsson]

@Matt J: Unfortunately you also removed the fft of the orignal signal (x all ones) which only have a non-zero amplitude at DC - which it corresponds to a Dirac-spike at 0 as it should. When zero-padding such a signal one fundamentally change the signal from a constant-valued signal (that in the fft-sense also is periodic over the length of the signal) to a square-wave with the period of the zero-padded signal-length (which obviously have a sinc as its fft.)

@Paul, 1: yeah, my rant now definitely comes raging out of nowhere. 2: Yes given that problem-description a rectangular window is natural, but "on an exam" I'd hope the next question is something along the lines of "explain how the implicit periodic assumption of the fft leads to aliasing of high-frequency components due to end-effects, and how to control that problem with windowing".

...and if I think we might kind of argue not so much about how many angels should be on the head of this pin, but what colour their wings should be?

This might be a script that illustrates the windowing issues a bit:

t = linspace(0,5*pi,512);
f1 = fftfreq(mean(diff(t(1:end-1))),numel(t(1:end-1)));
f0 = fftfreq(mean(diff(t)),numel(t));

figure
subplot(4,2,1)
% Signal that's periodic and its spectral amplitudes
plot(t(1:end-1),sin((2)*t(1:end-1)))  
subplot(4,2,2)
plot(f1,log10(abs(fftshift(fft(sin((2)*t(1:end-1)))))),'.-')  
axis([-5 5 -16 1.2*log10(numel(t))])
subplot(4,2,3)
% Signal that's one sample off being periodic
plot(t(1:end),sin((2)*t(1:end)))                            
subplot(4,2,4)
plot(f0,log10(abs(fftshift(fft(sin(2*t(1:end)))))),'.-')      
axis([-5 5 -5 1.2*log10(numel(t))])
subplot(4,2,5)
% signal a quarter period from being periodic
plot(t(1:end-1),sin((2.5)*t(1:end-1)))                      
subplot(4,2,6)
plot(f1,log10(abs(fftshift(fft(sin((2.5)*t(1:end-1)))))),'.-')
axis([-5 5 -5 1.2*log10(numel(t))])
w1 = tukeywin(numel(t(1:end-1)),0.75);
w2 = chebwin(numel(t(1:end-1)),50);
subplot(4,2,7)
% Windowed versions and their spectral amplitudes
plot(t(1:end-1),sin((2.5)*t(1:end-1)).*w1')
hold on
plot(t(1:end-1),sin((2.5)*t(1:end-1)).*w2')
subplot(4,2,8)
plot(f1,log10(abs(fftshift(fft(w1.'.*sin((2.5)*t(1:end-1)))))),'-') 
hold on
plot(f1,log10(abs(fftshift(fft(w2.'.*sin((2.5)*t(1:end-1)))))),'-')
axis([-5 5 -5 1.2*log10(numel(t))])


function f = fftfreq(dt,n)
% FFTFREQ - fft-frequency-array
%   
% Calling:
%   f = fftfreq(dt,n)
% Input:
%   dt - sampling duration (s), double scalar
%   n  - number of frequencies,  scalar int
% Output:
%   f - frequency array (Hz), dobule array

if rem(n,2) == 0
  f = [0:n/2-1,-n/2:-1] / (dt*n); %   if n is even
else
  f = [0:(n-1)/2,-(n-1)/2:-1] / (dt*n); %   if n is odd
end
end

Should produce a figure something like this:

[Matt J]

When zero-padding such a signal one fundamentally change the signal from a constant-valued signal to a square-wave

@Bjorn Gustavsson I see your point, but I think it's a distinction without a difference. The Dirac spectrum that you get from a constant-valued signal can also be seen as a poorly sampled version of the square wave's spectrum. In fact, that's probably a better way to view it, because a proper Dirac spectrum, such as you would get from a true, constant-valued continuous signal would be infinite at DC whereas the spectrum of a discrete sequence of 1's is not..

[Bjorn Gustavsson]

@Matt J, yeah, that is true, but in my defense seeing a fourier-transform of a constant-valued function with amplitudes for components at any frequency except at dc just make my brain hurt - and I know about some of the definitions of the Dirac, so that's not a problem.

[Matt J]

@Paul Yes, I've done it.

[Paul]

@Bjorn Gustavsson

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.

[Rahul Kale]

Hello Paul,

I would really appreciate your response to my questions, especially on the scale factor to be used for a zero padded signal (3 or 6?)

Thanks in advance

Rahul

[Rahul Kale]

Saying in other words, I think if the function is actually zero valued or not in the zero paddded region is determined by the period we define. So, in your example if x1 is [111000], then period N=6. If I just reduce N to 4 then rest two zeros are due to padding!

Is my understanding correct? If yes, then the normalizing factor becomes subjective.

Thanks

Rahul

[Matt J]

@Rahul Kale Paul has already explained to you that the FFT does not view x1 as a single period of a longer periodic signal. The FFT views x1 as a finite sequence. If your signal consists of multiple pulses, the pulses must be explicitly present in the input x vector if you want the FFT to be aware of them. In other words, if you want to find the spectrum of a 3 pulse signal with 50% duty cycle, you would do,

x=[1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0];

X=fft(x)

The spectrum X and its amplitude will indeed depend on the duty cycle, i.e., if I add or remove zeros between the pulses, the amplitude will change. However, increasing the number of zeros between the pulses is not the same thing as zero-padding. Zero-padding refers to adding initial and/or trailing zeros exclusively, and has no effect on the amplitude of the spectrum, as I demonstrated to you in my original answer above.

[Rahul Kale]

Matt,

Thanks for your explanation. I do not have any confusion if the zeros are added in between. I am interested to know if the signals actually has zeros at the beginning or end then how to distinguish those from zeros due to padding.

Thanks

Rahul

[Rahul Kale]

Matt...Just forgot to add, how does FFT know that in your signal x(with 50% duty cycle) the last three zeros are actually present in the signal and is not zero padding?

Thanks

Rahul

[Matt J]

how does FFT know that in your signal x(with 50% duty cycle) the last three zeros are actually present in the signal and is not zero padding?

It doesn't need to know. You are free to view those zeros either as part of the signal or as padded zeros. The result is the same.

[Matt J]

Below I have modified my original example to show the effect of varying the duty cycle of a series of dirac pulses on the amplitude spectrum. I was incorrect to say that the amplitude changes, even though as we see, the spectrum does change overall.

x=[1,1,1,1 1];
T=1:3;
for t=T
    plotFcn(x,t); %insert t-1 zeros between the x(i) and plot spectrum
end
hold off; legend("duty cycle="+100./T+"%")

function plotFcn(x0,t)

 x=zeros(1,t*numel(x0));  
 x(1:t:end)=x0;

 N=numel(x)+1e4;
 freqAxis=((0:N-1)-ceil((N-1)/2))/N;
 
 X=fft(x,N);

 
 plot(freqAxis, fftshift(abs(X)),'-'); hold on
 xlabel 'Freq (Hz)'
 ylabel 'Spectrum Amplitude'
end

[Paul]

@Matt J

To the contrary, I said quite the opposite: ".. the results [of the DFT] are always the Discrete Fourier Series coefficients of the periodic extension of x[n] ..." To be clearer, I should have said "DFS coefficients (to within a scale factor)"

We can see this in this example.

Let the signal x1[n] = [1 1 1 0 0 0] for n = 0:(N1-1), N1 = 6, and zero otherwise

Let the signal x2[n] = [1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0] for n = 0:(N2-1), N2 = 18, and zero otherwise.

Let x1p[n] be the periodic extension of x1[n] with period N1.

Let x2p[n] be the periodic extension of x2[n] with period N2.

It should be clear that x1p[n] == x2p[n]

Consequently, the DFTs of one period x1p[n] and x2p[n] should contain the same information

x1 = [1 1 1 0 0 0]; N1 = 6; % one period of x1
x2 = [x1 x1 x1];    N2 = 18;% one period of x2
X1 = fft(x1); f1 = (0:(N1-1))/N1;
X2 = fft(x2); f2 = (0:(N2-1))/N2;
figure
subplot(211);plot(f1,abs(X1),'-o',f2,abs(X2),'-x');
subplot(212);plot(f1,angle(X1),'-o',f2,angle(X2),'-x');

The ampltidues are identical to within a factor of N2/N1 = 3, and the phases are identical. The DFT of one period of x2[n] is the identical (to within that scale factor of 3) as DFT of one period of x1[n], because the periodic extensions of x1[n] and x2[n] are the same. By "identical" I mean wrt the four non-zero points that are encoding all of the information about the underlying periodic signals.

[Matt J]

@Paul Very well. Still, another way to view the DFT is the Fourier transform of the continuous signal formed from the sinc interpolation of the finite sequence ,

and then sampling at frequencies f=k/N, k=-N/2,1,...,N/2-1.

The continuous signal so defined is indeed not precisely finite, but it does decay to zero and is certainly not periodic. So, if the spectrum you are trying to compute is the continuous Fourier transform of a multi-pulse signal, the DFT will not be aware of the presence of multiple pulses unless they are explicitly sampled in .

[Rahul Kale]

Paul,

Thanks for the example.

Would it be possible for you to superimpose FFT of [x2 0 0 0 0 0 0] on the plots that you have given above? I would like to see what happens to the spectrum after zero padding in MATLAB when the signal has trailing zeros. I am not quite sure I agree with what Matt has said about the same FFT output irrespective of whether we consider the trailing zeros as a part the signal or a part of the zero padding.

Thanks

Rahul

[Matt J]

Would it be possible for you to superimpose FFT of [x2 0 0 0 0 0 0] on the plots that you have given above?

I have already done that for you in my last set of plots, except that instead of adding 6 zeros, I have added 10000. These were the lines of code.

 N=numel(x)+1e4;
 freqAxis=((0:N-1)-ceil((N-1)/2))/N;
 
 X=fft(x,N); %zero pads x to N=numel(x)+1e4

[Paul]

@Rahul Kale

Sure, we can do that.

x1 = [1 1 1 0 0 0]; N1 = 6; % one period of x1
x2 = [x1 x1 x1];    N2 = 18;% one period of x2
x3 = [x2 zeros(1,6)]; N3 = N2 + 6; % requested by Rahul

X1 = fft(x1); f1 = (0:(N1-1))/N1;
X2 = fft(x2); f2 = (0:(N2-1))/N2;
X3 = fft(x3); f3 = (0:(N3-1))/N3;

Let's also compute the DTFT of x1[n] and x2[n]. The DTFT of x3[n] is exactly the same as that of x2[n], so we don't need to compute it.

[X1DTFT,f1dtft] = freqz(x1,1,2048,'whole',1);
[X2DTFT,f2dtft] = freqz(x2,1,2048,'whole',1);

First. plot the DFT's, and I'll muliplty the DFT of X1 by 3 to get everything on equal footing.

figure
hold on
plot(f1,abs(3*X1),'-o');
plot(f2,abs(X2),'-x');
plot(f3,abs(X3),'-*');
legend('X1*3','X2','X3')

We see now that the zero-padding on x3 has added more frequency domain samples, which is expected because N3 > N2. Why are those X3 frequency domain samples where they are? Well, they are the samples of the DTFT of x3, which as we stated above is identical to the DTFT of x2. Add the DTFTs of x1 and x2 to the plot.

figure
hold on
plot(f1,abs(3*X1),'-o');
plot(f2,abs(X2),'-x');
plot(f3,abs(X3),'-*');
plot(f1dtft,abs(X1DTFT)*3,f2dtft,abs(X2DTFT));
legend('X1*3','X2','X3','X1DTFT*3','X2DTFT')

Zoom in a bit for a better view

copyobj(gca,figure)

xlim([0 0.4])
legend('X1*3','X2','X3','X1DTFT*3','X2DTFT')

As must be the case, the samples of X1 lie on X1DTFT. Also, as must be the case, the samples of X2 and X3 lie on X2DTFT. As we can see, X2DTFT is sinc-y. The samples of X3 when viewed alone, as in the previous plot, might look kind of funky. But that's only because of where those X3 samples lie on X2DTFT (which is the same as X3DTFT if we were to compute it). If we add add more zero-padding to x3, the resulting X3 samples would look more and more sinc-y, just like the green curve.

Also, as is clear, 3*X1DTFT ~= X2DTFT. However, they are equal to each other at four frequencies, which are, not coincidentally, the frequencies of the non-zero samples of their respective DFTs.

[Paul]

@Matt J

Responding to this comment ....

Focusing on all but the last sentence, I agree. We can illustrate with a simple example.

syms t w f real
xn = 1:5; n = 0:4;
x(t) = sum(xn.*sinc(t-n));
figure
fplot(x(t),[-3 10])
hold on
stem(n,xn)

x[n] is a five-sample window of the infinite duration signal x(t),

Going to frequency domain

First the CTFT of x(t)

X(w) = simplify(fourier(x(t),t,w)); % CTFT
X(f) = X(2*sym(pi)*f);

Now the DTFT of x[n]

[Xdtft,fdtft] = freqz(xn,1,linspace(-1.5,1.5,500),1);

Now the DFT of x[n]

Xdft = fftshift(fft(xn));
fdft = (-2:2)/5;

Now plot. I'll lift up Xdtft a bit just so we can see what's going on underneath it

figure
hold on
fplot(abs(X(f)),[-1.5 1.5])
plot(fdtft,1.05*abs(Xdtft))
stem(fdft,abs(Xdft))
legend('CTFT','DTFT','DFT')

As expected, because of the specific definition of x(t), the DTFT of x[n] is an exact periodic extension of X(f).

So far, so good. But I don't see how the above relates to anything about inherit periodicity built into the DFT of a finite duration signal.

"So, if the spectrum you are trying to compute is the continuous Fourier transform of a multi-pulse signal, the DFT will not be aware of the presence of multiple pulses unless they are explicitly sampled in x[n]"

Please correct me if I'm wrong, but I think this statement is in regards to using sampling and the DFT to estimate the CTFT of a finite duration, multi-pulse signal (as opposed to an infinite-duration, periodic signal)

So something like this

x(t) = rectangularPulse(-0.5,2.5,t);
x(t) = x(t) + x(t-6) + x(t-12); % finite duration signal of three pulses
figure
hold on
fplot(x(t),[-1, 20]),grid
axis padded
% samples 
xn = repmat([1 1 1 0 0 0],1,3); n = 0:17;
stem((0:17),xn);

CTFT of x(t)

X(f) = subs(simplify(fourier(x(t),t,w)),w,2*sym(pi)*f);
figure
fplot(abs(X(f)),[-0.75 0.75])
ylim([0 9])

If we want to get the best estimate of X(f), we'd use xn, which contains samples of all three pulses, and compute the DTFT

figure
hold on
fplot(abs(X(f)),[-0.75 0.75])
[Xdtft,fdtft] = freqz(xn,1,linspace(-0.5,0.5,1001),1);
plot(fdtft,abs(Xdtft))

The DTFT is a pretty good approximation to the central part of the CTFT, i.e., -0.5 <= f <= 0.5. The DTFT is actually the same as what we'd get with fft(xn,1001), i.e., a massively zero-padded DFT of xn, so we wont plot that here.

Now, we know from above that the DFT of xn are frequency domain samples of the DTFT. Add those to the plot

Xdft = fft(xn);
stem((-9:8)/18,abs(fftshift(Xdft)))

As shown, the DFT of xn are frequency domain samples of the DTFT.

But we can get the exact same (to within the factor of 3) DFT information using only the first pulse

Xdft1 = fft(xn(1:6));
plot((-3:2)/6,3*abs(fftshift(Xdft1)),'-x')
legend('CTFT','DTFT','DFT of three pulses','3*DFT of one pulse')

If the objective is to estimate the CTFT of x(t) via a DFT w/o zero-padding, the only additional information we get using all three pulses is the amplitude scaling of the samples at the four frequencies with non-zero amplitude. Does this result align with what you expected from "the DFT [not being aware] of the presence of multiple pulses"? Or should not including those two additional pulses have a different or an additional effect on how the DFT of a single pulse compares to the CTFT of x(t)?

[Matt J]

@Paul I didn't follow all of that, but I can comment on this part:

Please correct me if I'm wrong, but I think this statement is in regards to using sampling and the DFT to estimate the CTFT of a finite duration, multi-pulse signal (as opposed to an infinite-duration, periodic signal)

Yes, because the CTFT of an infinite duration periodic signal is not physical. It will be the CTFT of a single pulse multiplied by a Dirac comb, which is basically zero or infinity at all frequencies. The OP's proposal here to analyze repeated bell pulses with the CTFT therefore only makes sense if there are a finite number of pulses.

[Paul]

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]

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.