如何从python中的离散数据中获得周期性曲线拟合？

Question

我想在像 ODE 这样的动态系统上进行集成

x_ddot + d*x_dot + k*x = a*sin(omega*t)

有一个修改：外力 a * sin(omega * t) 必须由另一个基于某些测量数据的周期信号代替。这意味着：我需要对我的离散测量数据进行曲线拟合，并且结果函数必须是周期性的。我有两个想法如何解决这个问题：

1) 使用傅里叶变换 (numpy.fft)。但是使用离散傅里叶变换很难从离散数据中生成连续函数。

2) 使用曲线拟合，函数如 a1+a2*sin(omega*t)+a3*sin(2*omega*t)+a4*sin(.... 其中 omega 为 2*pi/ （测量数据的长度）。不幸的是，这也不是很成功。我已经尝试结合 sin 和 cos 项并达到非常高的顺序 (...sin(10*omega*t)) 没有不要让它变得更好。

Plot

我真的很感激一些提示，请在下面查看我的示例数据代码：

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import math as math
import scipy.special as sp

# create data
f =    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1.1,1.3,1.7,1.9,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2.1,2.3,2.9,4.1,4.3,4.4,4.5,4.4,4,3.6,3.2,2.8,2.4,2.0,1.6,1.3,1.2,1.1,1,1,1,1,1]
alpha   = np.linspace(0,len(f),len(f))

omega = 2*np.pi/len(f)

def func(alpha, a1, ac2, ac3, ac4, ac5, ac6, ac7):
    fit     = a1 + ac2*np.sin(omega*alpha) + ac3*np.sin(omega*2*alpha) + ac4*np.sin(omega*3*alpha) + ac5*np.sin(omega*4*alpha) +      ac6*np.sin(omega*5*alpha) + ac7*np.sin(omega*6*alpha)
    return fit

popt, pcov = curve_fit(func, alpha, f) 

print popt

y_fit= func(alpha,popt[0],popt[1],popt[2],popt[3],popt[4],popt[5],popt[6])

plt.plot(alpha,f,'bo',label='discrete data')
plt.plot(alpha,y_fit,'r',label='periodic fit')
plt.legend()
plt.show()

Answer 1

如果您对数值函数感到满意，基于插值的简单方法如下：

第一步：根据数据集定义一个插值函数，确定函数的一个周期。 （我使用了来自 Scipy 的简单 interp1d，Scipy 中也提供了更复杂的替代方案。）

from  scipy.interpolate import interp1d    
f_interp_one_period = interp1d(alpha, f, 'cubic')

第2步：根据上述定义一个周期函数（我根据数据使用值55作为周期）

def f_periodic(alpha):
        return f_interp_one_period(numpy.mod(alpha,55))

下图显示了 -100 到 200 范围内的 f_periodic 以及原始数据