How do I calculate peak amplitude of the signal components after zero padding and FFT?

Question

I am learning about DFT and trying to apply it to some audio processing. I am new to DSP but experienced in programming and have some background in math and physics. The FFT algorithm I use (lomontFFT) doesn't have 1/N scaling in front of the DFT sum (ie. it's defined just like on wikipedia DFT page). I will use simple example to describe my issue.

Let's say I have have 1024 samples of the pure 100Hz sine signal with peak amplitude of 1, and let the sample rate be 10240Hz. So that means that frequency resolution is 10240Hz/1024 = 10Hz. So when I perform FFT I get peak in 10th frequency bin which is correct. (since the signal is real I am only looking at the first half of the spectrum). Now, the magnitude of this peak in 10th bin is 512. How should I interpret this value? Is this magnitude somehow related to the energy of the 100Hz sine wave component in the signal?

Now, how can I get a spectrum to show amplitudes? I'd like to see what components are contained in the signal and what are their amplitudes. So in this example I should have one component at the 10th frequency bin with the value of 1 which is equal to the amplitude. I know I can divide the current value I got by 512 which seems to be the N/2. Is this correct way to calculate it?

If so, what about zero padding? If I zero-pad the same signal with additional 1024 zeros now I have total of 2048 samples. This means that frequency resolution is now 5Hz and my sine component is in the 20th bin. The magnitude is again 512, however now if I divide it by N/2=2048/2=1024 I get 0.5 instead of 1 for the peak amplitude. It seems like I should not count zero padded samples in this scaling? But what about if I just get the signal from somewhere and I don't know if signal is zero padded?

Also it seems that zero padding changes magnitudes in other bins. So in the first case without zero padding I get magnitude of 512 in the 10th bin and all other magnitudes are nearly zero (eg. values of the order 10^-10). But in the second case with zero padding, beside the magnitude value of 512 in the 20th bin I also get large magnitudes in bins around it (eg. 19th bin has magnitude of 334 and 21th bin has magnitude of 318 etc...). I guess these are side-lobes of spectral leaking but why are they showing up if I only have one sine frequency of 100Hz and frequency resolution is 5Hz so it falls exactly in the 20th bin?

I wish I could show some plots but I am programming this in C# console application so I can't provide visuals.

Note: I am using "frequency resolution" term here to denote "frequency bin spacing", I am aware that zero-padding doesn't increases frequency resolution since it doesn't add any new information to the signal.

Answer 1

What you have observed is why I prefer a $1/N$ normalization factor. If your signal is a pure tone, it is very well behaved in a DFT. For a frequency which is a whole number of cycles within your sample frame, the magnitude of the corresponding bin will be $1/2$. This is due entirely to the fact that a sinusoid is the average of two complex signals. The most straightforward way to see it is in the exponential definition of the cosine function:

$$ \cos( \theta ) = \frac{ e^{ i\theta } + e^{ -i\theta } }{2} $$

This also explains why the upper half of the DFT is the complex conjugate mirror of the lower half for real valued signals.

When your pure tone is not a whole integer number of cycles in the sample frame a phenomenon known as "leakage" occurs. This is not a flaw, it is how the DFT works. You can find the equation for leakage values in my blog article: DFT Bin Value Formulas for Pure Real Tones. As far as I know, this is the only place you will find these exact formulas.

Zero padding distorts things.

You will find a bunch of articles concerning the DFT in my blog. I recommend you start at the beginning and read them all.

Answer 2

Please have a look at the following MATLAB / OCTAVE code that performs a DFT analysis of the windowed sine wave including adjustable zero padding. As you can see (without any theoretical analysis) the scaling factor is not related to the zero padded final length but the original nonzero samples length. Furthermore you can see that MATLAB's built in FFT also uses no $1/N$ weighting, as apparent from the explicitly computed result. If your DFT/FFT routine uses any weight $K$ then you should only consider that $K$ in addition. Note that in order to make the plots clear, I've chosen the signal frequency and sampling frequency in rational proportion so that the peak of the DFT coincides with the actual sine wave frequeny. Otherwise there will be further changes in the DFT peak amplitude due to the shift in location of the samples of DFT according to theoretical DTFT.

clc; clear all; close all;

% Analog signal of duration Td, sampled at Fs Hz, producing L samples:
Td = 0.1;               % observaton interval in seconds
f0 = 50;                % input sine frequency in Hz.
Fs = 500;               % sampling frequency in Hz.
t = 0:1/Fs:Td;          % sampling instants time
L = length(t);          % length of actually sampled signal

A = 1.5;                % amplitude of the analog sine wave
x = A * sin(2*pi*f0*t); % samples of the analog signal...

M = 2*L;                % M: zero padding (on those L samples of x[n])
N = L + M;              % N: final DFT / FFT length, which is zero padded to N samples
xw = [x , zeros(1,M)];  % zero padded signal generated explicitly


% Compute (inefficently) N-point DFT without (1/N) weighting
sum = 0;
Xk = zeros(1,N);
for k = 1:N
    sum = 0;
    for n = 1:N
        sum = sum + xw(n) * exp(-j*2*pi*(k-1)*(n-1) / N) ;
    end    
    Xk(k) = sum;
end

% Display the results:
figure,subplot(3,1,1)
stem(xw); 
title(['L = ', num2str(L), ' point sequence x[n], padded with M = ',...
       num2str(M), ' samples into length of N = ',num2str(N),' samples']);
subplot(3,1,2)
stem(linspace(-Fs/2,Fs/2,N) , abs(fftshift(fft( xw , N ))) / (L/2) );
title('Magnitude, N-point DFT (matlab) of zero padded xw[n]'); 
subplot(3,1,3)
stem(linspace(-Fs/2,Fs/2,N) , abs(fftshift(Xk) / (L/2) )  );
xlabel('Magnitude, N-point DFT (explicit) of zero padded xw[n]'); 

Answer 3

Not yet mentioned clearly so I will add that zero-padding is an excellent, fast and robust approach for making very accurate estimates of the peak amplitudes of the signal components.

In order to do this properly you will likely need to use a more advanced window to avoid the high spectral side-lobes, which can also contribute to errors in the stronger signals due to sidelobes from other signals (or in very low number of bins even with a single real tone as I will demonstrate below, due to the image of a real signal).

Further realize that for a general sinusoid of the form $Acos(\omega t)$, the result when scaled by the typical $1/M$ will result in a magnitude of $A/2$ since each bin in the DFT represents an $e^{j\omega t}$ not a cosine or sine. (Euler's identity relates the two where we see the half value). Without zero=padding, we can only get this result accurately if the signals are an even number of cycles within the DFT frame. With zero=padding this is not the case (see other posts on this site relating the DTFT to the DFT and how zero padding gives us more samples of the DTFT).

BEST APPROACH FOR QUICKLY GETTING AN ESTIMATE OF ALL AMPLITUDES FROM DFT:

Window AND Zero Pad!

The following summarizes my go-to approach for getting the correct amplitudes for each tone in a stationary $M$ sample sequence using an $N$ sample DFT with $N>>M$ (zero-padding):

First choose a window that will ensure low sidelobes. This is NOT the rectangular window! I like the Kaiser window since you can trade the main lobe width and sidelobe rejection directly through it's $\beta$ knob: using $ \tt kaiser(M, \beta\tt)$ for an $M$ sample window. The sidelobes are made to be low enough so that bins in closer proximity don't interfere - but know we lose frequency resolution in the trade (be careful about the cyclical nature of the DFT, bin $N-1$ is adjacent to bin $0$!

The window length $M$ is the number of actual samples, NOT samples with zero padding!!

The scaling that would have been $1/M$ is changed to be $1/S$ where $S$ is the summation of the window sammples. This properly compensates for the signal we remove in time from the window. Such using a windowed FFT, the DFT is computed using:

$$X[k] = \frac{1}{S}\sum_{k=0}^{N-1}x[n]w[n]e^{-j2\pi k n/N}$$

Where

$x[n]$ is the input sequence

$w[n]$ is the window

$S = \sum_{k=0}^{M-1}w[n]$: is the proper scaling factor based on window used

$M$ is the number of samples in the data and window

$N$ is the total number of samples with zero-padding ($M-N$ zeros appended to the time data)

---

DEMONSTRATION

Below is a MATLAB/Octave demonstration of two non-integer cycle tones at different amplitudes with only $21$ bins. The first sinusoid has an amplitude of $10$ and a frequency of $0.205 \text{ cycles/sample}$ (or $k = 4.305$ in bin counts), and the second sinusoid has an amplitude of $2$ and a frequency of $0.352 \text{ cycles/sample}$. ($k = 7.392$ in bin counts). As real signals, we would expect each tone to produce two impulses in frequency each with a magnitude of half of the amplitude, or $5$ and $1$ in this example.

Here is the MATLAB/Octave code:

M= 21;
N =200;
n=0:M-1;
f1 = .205;
f2 = .352;
sig1 = 10*cos(2*pi*f1*n);
sig2 = 2*cos(2*pi*f2*n);
out = sig1+sig2;

win = kaiser(M, 6);
fout = 1/sum(win)*fft(win'.*out, N);
faxis = (0:N-1)/N;

First for comparison, below is the result if we did neither zero-padding nor a more advanced window, This is specifically:

1/M*abs(fft(out))

Result with neither zero-padding or further windowing:

Next only the zero padding is added where we see a dramatic improvement in both the ability to determine the frequency and the amplitude, yet there is still a noticeable error, and our lower level signal is indistinguishable from the other side-lobes of the larger signal:

1/M*abs(fft(out, N))

Finally the code as above is run showing the significant improvement in providing a reasonable estimate of both the larger and smaller signal amplitudes. The actual results from the max of the two signals in the plot below came to:

Larger signal: $4.999$

Smaller Signal: $0.993$

Here even in an example like this with very small number of bins or similarly what would be relatively closely spaced tones in a larger data set we are still able to estimate the amplitude very accurately, the lower level signal had an error of only -0.06 dB which is more than suitable for most practical applications.

Answer 4

It looks to me like your frequency resolution is sufficient and you are getting the correct value of the peak magnitude in your FFT but your scaling factor is incorrect.

You are right in dividing the magnitude of your FFT peak by N (1024 in your case). This number however, divides the power across the positive and negative sides of the frequency spectrum.

If you only want to look at one side of the FFT, you should multiply your peak magnitude by 2/N to get the correct amplitude that you are expecting. See this answer for another example.

e.g.

512 * 2 / 1024 = 1 (the answer that you were expecting)