2
$\begingroup$

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 wavelengths of the tone. 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 corresponding 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 matlab 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);
$\endgroup$

1 Answer 1

3
$\begingroup$

The FFT of a strictly real input results in a conjugate mirrored result. For a sinewave, both a positive and a negative frequency peak will appear. If the input is not integer periodic over the full aperture or is zero-padded, the two response peaks will be Sinc shaped, with lots of extended ripples in the frequency domain. The positive and negative frequency Sinc images overlap and thus can constructively or destructively interfere, depending on the phase of the input sinewave and/or the window. For the lowest few FFT bins (and those very near N/2), this overlap can be quite substantial, and thus so is the potential peak shift of the summation.

$\endgroup$
3
  • $\begingroup$ Thanks! That could explain why certain windowing helped but didn't entirely alleviate the problem. Do you know of any source where I could read more about this problem/ phenomena? $\endgroup$ Aug 10, 2019 at 21:37
  • 1
    $\begingroup$ a complex sine would have a single component $\endgroup$
    – user28715
    Aug 11, 2019 at 2:30
  • 1
    $\begingroup$ Applying a Hilbert transform to filter negative frequencies solves the problem of finding the input frequency by peak picking of the fft M = 100; N = 100*M; m = ((0:M-1)/M)-1/2; f = 2; x = sin(2*pi*f*m); x2 = hilbert(x); FFT = fft(x2,N); figure,plot(abs(FFT)) [~,ind] = max(abs(FFT)); $\endgroup$ Aug 14, 2019 at 23:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.