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.