FFT interpolation using zero-padding and the Chirp Z-Transform for a single tone sinusoid

Why would FFT interpolation by zero-padding or using the Chirp Z-Transform produce a maximum at a bin that corresponds to a frequency less than the input frequency of a single tone sinusoid?

I am attempting to precisely determine the frequency of a single tone sinusoid from a data set which captures just 1-2 oscilations of the wave. To do this I have tried to use zero-padding interpolation and the Chirp Z-Transform. When I perform my interpolation the maximum in the FFT or CZT falls into a bin that does not correspond to the frequency of the input sinusoid and instead always undershoots. Further, a convolution between the input sinusoid and the complex exponential of the same frequency is smaller in magnitude than the convolution between the input and the complex exponential of the frequency coresponding to the maximum FFT bin. This obviously goes against the properties of sinusoids.

I have tried to diagnose this issue by windowing, using different interpolation factors, using different length and frequency input signals, and using other fft codes, to no success.

Below is the code that performs the interpolation of the fft of an input sinusoid, plots the fft and czt, plots the input sinusoid along with the sinusoid corresponding to the max bin, and finally performs the convolution between the input and the sinusoids of the FFT max bin and the input frequency.

%create input sinusoid
M = 100;
m = linspace(-pi,pi,M);
f = 2;
x = sin(f*m);

%set interpolation factor and run zeropadded fft
N = M*100;
FFT = fft(x,N);

%define czt parameters
r1 = 0;
r2 = .1;
a = exp(2j*pi*(r1));
w = exp(-2j*pi*(r2-r1)/N);

CZT = czt(x,N,w,a);

%Plot interpolated FFT and CZT
figure, plot(abs(FFT))
hold on
plot(abs((CZT)))
legend('fft','czt')

% Find index of maximum bin for FFT and CZT
[~,indexFFT] = max(abs(FFT(1:N/2)));
[~,indexCZT] = max(abs(CZT(1:N/2)));


% Create sinusoids from the extracted bins to compare against input
% sinusoid
a = exp(-2j*pi*(indexFFT-1)*(linspace(1,M,N))/(N));
b = real(a);

a2 = exp(-2j*pi*(indexCZT-1)*(linspace(1,M,N))/(N/(r2-r1)));
b2 = real(a2);

a3 = exp(-2j*pi*(f)*(linspace(1,M,M))/(M));
b3 = real(a3);

figure,plot(linspace(1,M,N),b,'-.')
hold on
plot(linspace(1,M,N),b2,'--')
plot(linspace(1,M,M),b3,':k')
hold off

legend('FFT','CZT','input sinusoid')
title('the sinusoid associated with the maximum FFT bin') 

% Perform convolution between zeropadded input and the sinusoid from the
% fft max bin, and compare against the convolution with the true sinusoid with true frequency

x_Padded = zeros(1,N);
x_Padded(1:M) = x;

X = 0;X2 = 0;

for n=1:N
    X =  X + x_Padded(n)*exp(-2j*pi*(f)*(n-1)/(M));
    X2 =  X2 + x_Padded(n)*exp(-2j*pi*(indexFFT-1)*(n-1)/(N));
end

X2_mag = abs(X2);
X_mag = abs(X);

Answer 1

[David Goodmanson]

Hi Matt,

I made one minor change at the beginning of your code, replacing the first four lines with

M = 100;
m = linspace(0,1,M)-1/2;
f = 2;
x = sin(2*pi*f*m);

This does the same thing, but putting the 2 pi into the sine argument is more consistent with common practice and with the rest of the code. And it's useful later.

Since you have zerofilled by a factor of 100, frequency f = 2 appears at point 200 in the FFT array. Plus one more point since Matlab is one-based and f=0 is at point 1. So, point 201 is the goal.

The problem is that the oscillations are for sine, but you are looking at the abs amplitudes of the positive frequency components and picking the largest one of those. Writing the continuous-m integral expression for simplicity [ and using f0 in place of f ] , the fourier amplitude for frequency fn is

an = Integral (1/(2*pi)) * ...
         [ exp(2*pi*i*f0*m) -  exp(-*pi*i*f0*m) ] * exp(-2*pi*i*fn*m) dm  

The first term does have max abs value at fn = f0, just like you expect. However, the second term also makes a contribution, and you are finding max abs of the sum of the two terms. That peaks out at the wrong spot, point 195.

Replacing the signal with

x = exp(2*pi*i*f*m)

i.e. positive frequency only, gives the correct answer. Except not quite. Now it peaks out at point 203. The problem is that the linspace array does not represent the periodic wave correctly, because it includes both end points. But if you go to

M = 100;
m = ((0:M-1)/M)-1/2;
f = 2;
x = exp(2*pi*i*f*m);

which includes only the first end point, then it all works out. The max is at point 201. The max of the chirp transform is at point 2001, which sounds right. The waves in the second plot all align perfectly.

[Matt Southwick]

Thanks! Really helpful awnser. Is there a book you would reccomened looking at to learn more about this, and the DFT in general?

[David Goodmanson]

Hi Matt,

I mostly figured it out on my own, having used continuous Fourier transforms a lot previously. You can get hosed trying to take every idea over to the discrete case, but I think it's a good starting point.

I don't have any DFT books, although I do have one pretty old FFT book that mentions not duplicating the end points, but only in passing.

Schaum's outlines have a couple of books on the subject, which I don't have. But they have excellent books on, for example, complex variables and topology, there are lots of worked problems, and the price is right. If I were going to buy one today sight unseen that would be it.

Maybe some other people will weigh in with DFT recommendations.

When I saw things peaking at frequencies slighty low, it immediately reminded me of the forced damped harmonic oscillator. There is a resonant frequency, but if you look at the amplitude as a function of forcing frequency, that function peaks below the resonant frequency, for I believe basically the same reason. Again that's a clue that came from the continuous case.

[Matt Southwick]

Thanks again David! I'll look into Shaum's.

As a further proof to your awnser and work-around for the discrete problem, I found that applying a hilbert transform to the signal to filter out the negative frequencies allows the input frequency to be determined by peak picking of the fft magnitude.

M = 100;
m = ((0:M-1)/M)-1/2;
f = 2;
x = sin(2*pi*f*m);
x2 = hilbert(x);

N = 100*M;
FFT = fft(x2,N);
figure,plot(abs(FFT))
[~,ind] = max(abs(FFT));

[David Goodmanson]

Hi Matt,

That's an interesting use of the hilbert transform. Of course what Mathworks calls the hilbert transform and what the rest of the world usually calls the hilbert transform are not the same (I wish Mathworks would call it something different), but their function gets the job done here.