通过样本恢复正弦函数的和

Question

我们有一个函数的 1000 个随机点 (x, y)，还有 x 的可能范围 a_n 和 b_n。给定此示例，以编程方式近似 a_n 和 b_n 变量的原始 10 长度数组。怎么样？

Answer 1

您可以将值重新采样到统一网格（恒定 x 步长），然后进行常规 FFT。

在理想情况下，您将得到一组峰值，其中第 i 个峰值的位置对应于频率 b(i)，幅度与 a(i) 相关（遗憾的是，在现实中不一定有效）

另一种方法 - 使用总和计算 non-uniform Fourier transform - 对于您的样本数 1000 和有限的频率范围来说，这是非常可靠的。请注意，您必须为每个需要的频率 b(k) 计算结果 F（与 a(k) 相关的幅度）。

F(b(k)) = Sum[j=0..N-1](Y[i]*exp(-2*i*Pi*X[j]*b(k)))

忘了说 - 也许有些频率太接近了，你可能需要一些优化方法，比如 Levenberg-Marquardt one 来分开它们。

Answer 2

函数scipy.optimize.curve_fit可以完成这项工作。如果 x 坐标均匀分布，可以对信号加窗并应用 DFT，或者可以通过 DFT 计算信号的自相关以确定其频率。

我们来 curev_fit() 试试看：

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

def func(x, a1, b1, a2, b2, a3, b3):
    # x is an array of x value
    return a1*np.sin(b1*b1*x)+a2*np.sin(b2*b2*x)+a3*np.sin(b3*b3*x)

n=1000
#np.random.seed(1729)
xdata = np.random.uniform(0,10,size=n)
y = func(xdata, 1.,1., 10.,3., -4.,1.7)

y_noise = 0.8 * np.random.normal(size=xdata.size)
ydata = y + y_noise

popt, pcov = curve_fit(func, xdata, ydata, p0=[1.,1.,2.,2.,3.,3.],method='lm')
print 'optimal parameters are ', popt

perr = np.sqrt(np.diag(pcov))
print 'standard deviation of parameters ',perr

#sorting xdata and ydata
for i in range(xdata.size):
    for j in range(i):
        if xdata[i]<xdata[j]:
             temp=xdata[i]
             xdata[i]=xdata[j]
             xdata[j]=temp
             temp=ydata[i]
             ydata[i]=ydata[j]
             ydata[j]=temp

plt.plot(xdata, ydata, 'b-', label='data')
plt.plot(xdata, func(xdata, *popt), 'r-',label='fit')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

提供的示例显示了一些技巧：

由于定义了函数，可以通过平方来限制参数为正。 np.sin(b1*b1*x)。约束ba也可以直接在curve_fit()中引入，但是它阻止了Levenberg-Marquardt方法的使用。
如 curve_fit() 的文档中所述，可以使用 pcov 中返回的值来估计已识别参数的误差。 perr = np.sqrt(np.diag(pcov))

如果为了拟合2*10个参数，这个方法在1000个点上运行，很可能一下子就不行了。您可以尝试减少正弦波的数量或将您可以获得的任何知识倒入拟合函数中。

通过样本恢复正弦函数的和

Restore sum of sinusoids function by its sample

python

math

trigonometry

fft