将波函数拟合到时间序列数据

Question

我有随时间推移采样的速度数据，我想找到一些可用于描述它的方程/函数。我的数据看起来像这个图

正如您从图像中看到的，数据不是简单的正弦或余弦波，我正在努力寻找一种方法来查找/拟合数据的某些方程。我的信号处理知识非常有限。

这里以下采样数据为例：

xdata = np.array([11.79, 11.87, 11.95, 12.03, 12.11, 12.19, 12.27, 12.35, 12.43,
       12.51, 12.59, 12.67, 12.75, 12.83, 12.91, 12.99, 13.07, 13.15,
       13.23, 13.31, 13.39, 13.47, 13.55, 13.63, 13.71, 13.79, 13.87,
       13.95, 14.03, 14.11, 14.19, 14.27, 14.35, 14.43, 14.51, 14.59,
       14.67, 14.75, 14.83, 14.91, 14.99, 15.07, 15.15, 15.23, 15.31,
       15.39, 15.47, 15.55, 15.63, 15.71, 15.79, 15.87, 15.95, 16.03,
       16.11, 16.19, 16.27, 16.35, 16.43, 16.51, 16.59, 16.67, 16.75,
       16.83, 16.91, 16.99, 17.07, 17.15, 17.23, 17.31, 17.39, 17.47,
       17.55, 17.63, 17.71, 17.79, 17.87, 17.95, 18.03, 18.11, 18.19,
       18.27, 18.35, 18.43, 18.51, 18.59, 18.67, 18.75, 18.83, 18.91,
       18.99, 19.07, 19.15, 19.23, 19.31, 19.39, 19.47, 19.55, 19.63,
       19.71, 19.79, 19.87, 19.95])

ydata = np.array([ 1.86470801e-05, -3.05185889e-03, -4.53502752e-03, -5.01501449e-03,
       -7.61753339e-03, -7.89916120e-03, -8.45710261e-03, -7.64640792e-03,
       -7.28613761e-03, -6.07402134e-03, -4.21708665e-03, -2.53126644e-03,
       -1.38970318e-03,  1.59394526e-05,  1.91879565e-03,  3.10854836e-03,
        3.37327421e-03,  4.56715556e-03,  5.59283055e-03,  7.38842610e-03,
        7.62706763e-03,  9.14228858e-03,  1.24410442e-02,  1.47384372e-02,
        1.50136837e-02,  1.14957746e-02,  9.03580024e-03,  7.04710182e-03,
        6.94429062e-03,  6.57961389e-03,  5.25393124e-03,  4.75389627e-03,
        2.49195903e-03,  2.41027520e-03,  1.67260849e-03,  5.02479585e-04,
        2.35275116e-04, -1.91204237e-03, -1.38367516e-03, -1.02639516e-03,
        2.78570931e-04, -4.42114657e-04,  2.86560704e-04,  5.69435589e-04,
       -2.94260316e-04,  6.68917718e-05, -1.61045579e-03, -2.32730345e-03,
       -2.55154534e-03, -4.01547893e-03, -4.72977248e-03, -5.12847064e-03,
       -7.10974182e-03, -6.99822300e-03, -8.07665704e-03, -8.84851229e-03,
       -9.09903075e-03, -1.03180810e-02, -1.14544281e-02, -1.21251714e-02,
       -1.30148026e-02, -1.29700705e-02, -1.25293732e-02, -1.12340153e-02,
       -1.05980279e-02, -8.51107110e-03, -6.05923880e-03, -4.50314785e-03,
       -3.11505449e-03, -1.61344075e-03, -1.90893378e-04,  8.53961182e-04,
        2.53364048e-03,  3.00061370e-03,  4.13971717e-03,  5.83572401e-03,
        7.97330030e-03,  9.11719022e-03,  1.07586221e-02,  1.20721136e-02,
        1.24900520e-02,  1.07001078e-02,  9.77437191e-03,  8.30650537e-03,
        6.17321981e-03,  4.27379777e-03,  3.53587465e-03,  2.37557043e-03,
        1.66114222e-03,  6.83006538e-04, -6.38602576e-04, -1.54135169e-03,
       -9.86915409e-04, -1.58287464e-03, -2.02820728e-03, -1.53065658e-03,
       -8.52157455e-04,  1.62949595e-03,  8.56897304e-04,  1.20745000e-03,
       -1.06239388e-03, -2.39230381e-03, -2.39669751e-03])

我的其他数据非常规则，我通常可以使用

scipy.optimize

来拟合 sin 或 cos 函数，如下例所示。

def cos_func(x, D, E):
    y = D*np.cos(E*x)
    return y

guess = [0.01, 4]
parameters, covariance = curve_fit(cos_func, xdata, ydata, p0=guess)
print(parameters)

fit_D = parameters[0]
fit_E = parameters[1]

fit_cosine = cos_func(xdata, fit_D, fit_E)

plt.plot(xdata, ydata, 'o', label='data')
plt.plot(xdata, fit_cosine, '-', label='fit')

但是在这个数据上却惨遭失败：

我应该用什么函数或类型的波来拟合这里的数据？我认为它需要一些具有不同幅度的正弦/余弦复合波才能创建这样的波，但我不知道如何在使用

scipy

优化参数之前自动找到正确的方程类型。

Answer 1

我强烈建议您不要将

curve_fit

用于周期性信号。例如，检查这篇文章，还有许多其他具有相同答案的示例：

curve_fit

无法准确估计频率，这意味着随着时间的推移，拟合和实际偏离。

相反，这似乎是应用高斯回归的一个非常好的数据集。这是一个例子：

# Do Gaussian Regression
kernel = 1 * RBF(length_scale=1.0, length_scale_bounds=(1e-2, 1e2)) # define the kernel
GP = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9) # define the GP model
x_train = x[::5] # get training x as every 5th point
y_train = y[::5] # get training y as every 5th point
GP.fit(x_train, y_train) # fit the GP
muPredictions, sigmaPredictions = GP.predict(x, return_std=True) # predict
# plotting
plt.scatter(x, y, label = "Original") # scatter original data
plt.plot(x, muPredictions, label="Mean predictions", color = "r") # plot the mean of predictions
plt.fill_between( # fill between confidance interval
    x.ravel(),
    muPredictions - 1.96*sigmaPredictions,
    muPredictions + 1.96*sigmaPredictions,
    alpha=0.5,
    label=r"95% confidence interval",
)
plt.legend() # add legend
plt.grid() # and grid

结果将如下所示：

您可以根据选择的训练点的稀疏程度来使其更细或更粗，目前是每 5 个元素。

以下是进口：

import numpy as np
%matplotlib notebook
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF

看来StackOverflow真的很喜欢动画，所以这里展示一下采样的效果：

我认为最佳值是 3：