用傅立叶分量拟合时间序列:估计傅立叶级数系数

Question

问题:我有一组显示出周期性变化的测量值(时间、测量值、误差)，我想用以下形式的傅立叶级数来拟合它们

​

​

其中A0是我测量的平均值，t是时间，t0是(已知)参考时间，P是(已知)周期。我想要拟合系数A_k和phi_k。

这是我目前掌握的情况：

# Find Fourier components
# nfourier is the number of fourier components 

    def fourier(coeffs, time_data, epoch, period, nfourier, A0):
       import numpy as np
       omega = 2.0*np.pi/period
       fseries = np.zeros(len(time_data))
       fseries.fill(A0)
       for k in range(nfourier):
          ak = coeffs[k]
          phik = coeffs[k+1]
          time_diff = time_data - epoch
          fseries = fseries + ak * np.cos(k * omega * time_diff + phik)

       return fseries

我估计残差如下：

def residuals(coeffs, measurement_data, time_data, error_data, epoch, period, nfourier, A0):
   model = fourier(coeffs, time_data, epoch, period, nfourier, A0)
   result = measurement_data - model
   return result

然后我将它与以下内容相匹配：

def fit_it(coeffs, measurement_data, time_data, error_data, epoch, period, nfourier, A0):
   from scipy.optimize import leastsq
   opt_coeff = leastsq(residuals, coeffs, args=(measurement_data, time_data, error_data, epoch, period, nfourier, A0))
   return opt_coeff

程序成功完成，但拟合似乎失败，如下图所示：

​

我不确定我在这里做错了什么，但也许专家可以提供一些建议。如果有人愿意提供帮助，我很乐意提供测试数据集。

Answer 1

我不明白傅立叶拟合背后的原因。我想你想知道你的数据的调制频率分量。我建议你每次取数据的平均值并进行fft，它会让你对你的数据的频谱有更多的了解。

就您的代码而言，我想提出两点意见。首先，第k个元素的相位是k+1个元素的振幅。其次，error_data没有从函数残差中获取任何东西。你可以检查一下这些要点。

这更像是一个评论，但我没有足够的名气来发表评论。我只是想帮你。

问候

Answer 2

如果你知道数据的周期，你应该对x轴进行相位折叠。看起来x轴是在朱利安日。您应该在测量过程中计算相位。

阶段=((时间-参考时间)%周期)/周期

您将希望将傅立叶级数拟合到测量与相位图中，它看起来更像是一个周期信号。