I am trying to estimate the amplitude, frequency, and phase of an incoming signal of about 50Hz based on measurement of only a few cycles. The frequency needs to be precise to .01 Hz. Since the signal itself is going to be a pretty clear sine wave, I am trying parameter fitting with SciPy's curve_fit. I've never used it before, so I wrote a quick test function.
I start by generating samples of a single cycle of a dummy cosine wave
from math import *
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
fs = 1000 # Sampling rate (Hz)
T = .1 # Length of collection (s)
windowlength = int(fs*T) # Number of samples
f0 = 10 # Fundamental frequency of our wave (Hz)
wave = [0]*windowlength
for x in range(windowlength):
wave[x] = cos(2*pi*f0*x/fs)
t = np.linspace(0,T,int(fs*T)) # This will be our x-axis for plotting
Then I try to fit those samples to a function, adapting the code from the official example provided by scipy: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html
# Define function to fit
def sinefit(x, A, ph, f):
return A * np.sin(2*pi*f * x + ph)
# Call curve_fit
popt,cov = curve_fit(sinefit, t, wave, p0=[1,np.pi/2,10])
# Plot the result
plt.plot(t, wave, 'b-', label='data')
plt.plot(t, sinefit(t, *popt), 'r-', label='fit')
print("[Amplitude,phase,frequency]")
print(popt)
This gives me popt = [1., 1.57079633, 9.9] and the plot
My question is: why is my frequency off? I initialized the curve_fit function with the exact parameters of the cosine wave, so shouldn't the first iteration of the LM algorithm realize that there is zero residual and that it has already arrived at the correct answer? That seems to be the case for amplitude and phase, but frequency is 0.1Hz too low.
I expect this is a dumb coding mistake, since the original wave and the fit are clearly lined up in the plot. I also confirmed that the difference between them was zero across the entire sample. If they really were .1 Hz out of phase, there would be a phase shift of 3.6 degrees over my 100ms window.
Any thoughts would be very much appreciated!
The problem is that your array t is not correct. The last value in your t is 0.1, but with a sampling period of 1/fs = 0.001, the last value in t should be 0.099. That is, the times of the 100 samples are [0, 0.001, 0.002, ..., 0.098, 0.099].
You can create t correctly with either
t = np.linspace(0, T, int(fs*T), endpoint=False)
or
t = np.arange(windowlength)/fs # Use float(fs) if you are using Python 2