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.
