Простой тест scipy curve_fit не дает ожидаемых результатов - PullRequest
Новая
       20

Простой тест scipy curve_fit не дает ожидаемых результатов

0 голосов
/ 21 сентября 2018

Я пытаюсь оценить амплитуду, частоту и фазу входящего сигнала около 50 Гц на основе измерения только нескольких циклов.Частота должна быть точной до 0,01 Гц.Поскольку сам сигнал будет довольно четкой синусоидальной волной, я пытаюсь подгонять параметры под SciPy's curve_fit.Я никогда не использовал его раньше, поэтому я написал функцию быстрого теста.

Я начинаю с генерации выборок из одного цикла фиктивной косинусной волны 

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

Затем я пытаюсь соответствоватьэти образцы к функции, адаптируя код из официального примера, предоставленного 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)

Это дает мне popt = [1., 1.57079633, 9.9] и сюжет

вывод графика

Мой вопрос: почему моя частота отключена?Я инициализировал функцию curve_fit точными параметрами волны косинуса, поэтому разве первая итерация алгоритма LM не должна понимать, что существует нулевой остаток и что он уже пришел к правильному ответу?Это похоже на случай амплитуды и фазы, но частота слишком низкая 0,1 Гц.

Я ожидаю, что это глупая ошибка кодирования, так как исходная волна и подгонка четко выстроены на графике.Я также подтвердил, что разница между ними была нулевой по всей выборке.Если бы они действительно были на частоте 0,1 Гц в противофазе, в моем окне 100 мс был бы сдвиг фазы на 3,6 градуса.

Будем очень благодарны за любые мысли!

1 Ответ

0 голосов
/ 21 сентября 2018

Проблема в том, что ваш массив t неверен.Последнее значение в вашем t равно 0,1, но с периодом выборки 1 / фс = 0,001, последнее значение в t должно быть 0,099.То есть время 100 выборок составляет [0, 0,001, 0,002, ..., 0,098, 0,099].

Вы можете правильно создать t с помощью 

t = np.linspace(0, T, int(fs*T), endpoint=False)

или 

t = np.arange(windowlength)/fs  # Use float(fs) if you are using Python 2
...