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

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

以下是下采样数据的示例：

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])

我的其他数据非常规则，我通常可以使用 sin 或 cos 函数，scipy.optimize如下例所示。


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。

14

  • 1
    首先对其进行 FFT 以获得主频率可能是个好主意。然后对其进行低通滤波以去除噪声。


    – 


  • 编辑标题并删除“复杂”，这可能会产生误导


    – 

  • 您是否拥有比此处显示的更多的可用数据？任何拟合重要波形的方法都将从使用更多波长中受益匪浅。 lastchance 建议使用 FFT 是一个很好的建议，但它不会告诉您太多您不知道的信息，除非您有更多的数据。


    – 

  • 1
    该波形的基本过程是什么以及什么“模型”有意义？我们讨厌猜测。需要详细信息。


    – 


  • 1
    您能否提供更多的数据周期，@jpmorr。这两个周期之间存在足够的差异，因此很难选择最合适的模型。老实说，我仍然会选择对任何周期信号进行 FFT。只需保留几个主要频率的组合（而不仅仅是单个正弦曲线），因为它会自然地扩展您的更简单的数据。有陡峭水波的理论（艾里波 -> 椭圆形波），但你的波浪看起来就像正在破裂。


    – 



2 个回答
2

确实如此。听取lastchance 的评论并执行简单的FFT 缩减。此处显示了几种不同的缩减长度：

import numpy as np
from matplotlib import pyplot as plt

xdelta = 0.08
xdata = np.arange(11.79, 20, xdelta)
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,
))

freq = np.fft.rfftfreq(n=xdata.size, d=xdelta)
fft = np.fft.rfft(a=ydata)
plt.plot(freq, np.abs(fft))

plt.figure()
plt.plot(xdata, ydata, label='full')
for length in range(14, 4, -3):
    fft[length:] = 0
    inverted = np.fft.irfft(fft)
    plt.plot(xdata[:-1], inverted, label=length)
plt.legend()
plt.show()

光谱

缩减时间序列

方程

当outer()和@都表示矩阵乘法时，实数 FFT 的逆表达式为

cut = 5
n = xdata.size

inverse = 2/n * (
    np.exp(
        2j*np.pi/n * np.outer(
            np.arange(n),
            np.arange(cut),
        )
    ) @ fft[:cut]
).real

3

  • 好的，这非常适合，但是拟合曲线的方程是什么 – 该方程对我来说是最重要的部分。


    – 

  • 正弦+余弦之和从根本上定义了傅里叶级数：维基百科中以欧拉形式显示的逆变换：


    – 


  • @jpmorr 添加了反方程。 Stack Overflow 不支持 MathJax，因此我将其显示为直接 numpy 表达式。


    – 

我强烈建议您不要使用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：

7

  • 这看起来是一个非常好的解决方案，但我实际上无法让代码工作。除了必须在各处重塑数组之外，将训练数据设置为子集还会给预测器带来问题ValueError: X has 804 features, but GaussianProcessRegressor is expecting 161 features as input.– 所以我不明白如何训练一部分数据？


    – 

  • 您可以在某处发布完整数据的链接吗？我可以稍后看一下。如果我记得的话，我必须使用 x.reshape(-1,1) 来重塑 x 。 y 保持原样


    – 

  • 1
    基本上，这些函数需要一列值，而当前您的数据存储为一行。因此，xdata.reshape(-1,1)从一开始就去做。


    – 


  • 有没有办法得到拟合曲线的方程？这就是我需要包含在模型中来表示的内容


    – 

  • 我不再靠近我的笔记本电脑，但高斯回归通常如下所示： P(y|f,X)~N(y|Hβ+f,σ^2.I)


    –