zero padding shifting signal peak?

Question

I am trying to find the fundamental frequency of a low-frequency signal. I need an estimate that is precise to .01 Hz based on only a few cycles, so I'm trying to code up a fft in Python. The signal itself will be almost a pure tone, so I'm testing it out on a simple simulated wave.

I'm not seeing the results I expect, which may just be a dumb coding mistake. I'm new enough to DSP, though, that it's also possible I'm fundamentally misunderstanding what the transform is doing. Any help would be very much appreciated!

from math import *
import numpy as np
import matplotlib.pyplot as plt

I start by setting some system parameters

fs = 10000 # sampling rate (Hz)
T = 0.1 # length of collection (s)
assert fs*T % 1 == 0 # make sure we have an integer number of samples
windowlength = int(fs*T) # total number of samples

Then I create a sampled cosine wave for testing. Its period isn't an integer multiple of my window, so I'll need to zero-pad

f0 = 20.1 # fundamental frequency in Hz
checkwave = [0]*windowlength
for x in range(windowlength):
    checkwave[x] = cos(2*pi*f0*x/fs)

Plotting my wave in the time domain looks fine

t = np.linspace(0,T,int(fs*T)) # time domain x axis
plt.figure()
plt.plot(t,checkwave)

I then carry out the zero-padding and take the dft with numpy's rfft module

zp = 10 # factor to use in zero padding
dft = np.fft.rfft(checkwave,zp*len(checkwave))
fmaxindex = np.abs(dft).argmax() # find the peak frequency's index
fmax = fs*fmaxindex/(zp*len(checkwave)) # calculate the peak frequency
print(fmaxindex) # Higher than I expected
print(fmax)

But my plot in the frequency domain is off. I thought I should be looking at a sampled sinc wave centered on 20.1Hz, and since I'm sampling well above the Nyquist limit the centering should be perfect, right? But, if that were the case, 20Hz should be the largest bin (translates to bin 20). Instead, the largest bin is 21Hz. Plotting, we can confirm that's the case.

f = np.linspace(0,fs/2,int(fs*T*zp/2)+1)# frequency domain x axis
freqspan = 4 # number of samples to plot around the peak, so we get some detail
plt.figure()
plt.plot(f[fmaxindex-freqspan:fmaxindex+freqspan+1],abs(dft[fmaxindex-freqspan:fmaxindex+freqspan+1]))

My first thought was I made an off-by-one indexing error, but when I reproduce my last two blocks of code with an extra factor of ten on the zero padding, I see that things get even worse.

zp = 100 # factor to use in zero padding
dft = np.fft.rfft(checkwave,zp*len(checkwave))
fmaxindex = np.abs(dft).argmax() # find the peak frequency's index
fmax = fs*fmaxindex/(zp*len(checkwave)) # calculate the peak frequency
print(fmaxindex) # Higher than I expected
print(fmax)

Plotting in the frequency domain:

f = np.linspace(0,fs/2,int(fs*T*zp/2)+1)# frequency domain x axis
freqspan = 6 # number of samples to plot around the peak, so we get some detail
plt.figure()
plt.plot(f[fmaxindex-freqspan:fmaxindex+freqspan+1],abs(dft[fmaxindex-freqspan:fmaxindex+freqspan+1]))

Now the peak frequency is 20.8 (bin 208), about 6 bins away from where the sinc should be putting it.

Thanks in advance for any thoughts.

Answer 1

The FFT of a strictly real input signal is conjugate symmetric. What that means is that, for a single sinusoidal frequency, the FFT result includes not just one positive frequency peak, but also a another peak at the negative frequency. For signals that are not strictly integer periodic in the FFT length, the peak can be a very broad Sinc after zero-padding. For extremely low numbers of signal periods (and for frequencies that are very near Fs/2), the negative frequency peak is close to and thus can interfere with the positive frequency peak and distort its exact location.

So, for very small numbers of periods that are not known to be strictly integer periodic in aperture, it may be better to use a parametric estimator of some sort, rather than an FFT, even windowed or zero-padded.