fourier3系列的数据拟合总是产生一条直线。

Question

我有我想要拟合Fourier3系列的数据，我寻找了这个答案：这里，并尝试了来自不同包的不同算法(比如symfit，和scipy)。但是当我绘制数据时，不同的包给我的结果是：在这里输入图像描述

目前，我使用的是curve_fit包，下面是我的代码：

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

def fourier(x, *as_bs):
    sum_a = 0
    sum_b = 0
    j = 1
    w = as_bs[0]
    a0 = as_bs[1]
    for i in range(2, len(as_bs)-1, 2):
        sum_a += as_bs[i] * np.cos(j * w * x)
        sum_b += as_bs[i+1] * np.sin(j * w * x)
        j = j + 1

    return a0 + sum_a + sum_b

T = pd.read_excel('FS_data.xlsx')
A = pd.DataFrame(T)

xdata = np.array(A.iloc[:, 0])
ydata = np.array(A.iloc[:, 1])

# fits
popt, pcov = curve_fit(fourier, xdata, ydata, [np.random.rand(1)] * 8)
print(popt)

data_fit = fourier(ydata, *popt)
print(data_fit)

plt.plot(ydata)
plt.plot(data_fit, label='after fitting')
plt.legend()
plt.show()

因此，我的代码基本上会读取随机8个数字，并将它们分别赋值为(f，a0，a1，b1，a2，b2，a3，b3)的初始猜测。

我试着在Matlab上对数据进行拟合，以检查数据是否可以与fourier3相匹配，结果非常好：在这里输入图像描述。

我在Python和Matlab上都打印了输出来进行比较，下面是两者的结果: Python：

w = 5.66709943e-01
a0 = 3.80499132e+01 
a1 = 5.56883486e-04
b1 = -3.88408379e-04
a2 = -3.88408379e-04
b2 = 3.32951592e-04
a3 = 3.15641900e-04
b3 = 1.96414168e-04

Matlab：

   a0 =       38.07  (38.07, 38.08)
   a1 =      0.5352  (0.4951, 0.5753)
   b1 =     -0.5788  (-0.5863, -0.5714)
   a2 =     -0.3728  (-0.413, -0.3326)
   b2 =      0.5411  (0.492, 0.5901)
   a3 =      0.2357  (0.2226, 0.2488)
   b3 =     0.05895  (0.02773, 0.09018)
   w =   0.0003088  

因此，正如所指出的，只有a0的值是正确的，但其他值与Matlab相去甚远。那么为什么我要在Python中得到这个结果呢？我做错什么了？

以下是那些喜欢测试它的人的数据：

https://docs.google.com/spreadsheets/d/18lL1iMZ3kdaqUUtRDLNRK4A3uCPzOrXt/edit?usp=sharing&ouid=112684448221465330517&rtpof=true&sd=true

Answer 1

我不喜欢Matlab，所以我不知道，在估计非线性拟合的起始值时，Matlab fit会做哪些额外的工作。不过，我可以说，curve_fit一点也不做，即所有的值都假定为1的顺序。最简单的方法，应该是重新确定x轴到范围[0, 2 pi]。因此，OP的问题再次是错误的起始值。然而，重新标度需要知道要拟合的主波大约是数据集的宽度。此外，我们还需要假设所有其他的fit参数也都是1级的。幸运的是，情况就是这样，所以这样做是可行的：

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

xdat, ydat = np.loadtxt( "data.tsv", unpack=True, skiprows=1 )


def fourier(x, *as_bs):
    sum_a = 0
    sum_b = 0
    j = 1
    w = as_bs[0]
    a0 = as_bs[1]
    for i in range(2, len( as_bs ) - 1, 2 ):
        sum_a += as_bs[i] * np.cos( j * w * x )
        sum_b += as_bs[i+1] * np.sin( j * w * x )
        j = j + 1
    return a0 + sum_a + sum_b

"""
lets rescale the data to get the base frequency in the range of one
"""

xmin = min( xdat )
xmax = max( xdat )
xdat = ( xdat - xmin ) / (xmax - xmin ) * 2 * np.pi

popt, pcov = curve_fit(
    fourier,
    xdat, ydat,
    p0 = np.ones(8)
)
### here I assume that higher order are similar to lower orders
### but slightly smaller. ... hoping that the fit correts errors in
### this assumption

print(popt)
### scale back w noting that it scales inverse to x
print( popt[0] * 2 * np.pi  / (xmax - xmin ) )
data_fit = fourier( xdat, *popt )

如果我们不能进行上述假设，我们可能只假设有一个基频对信号的主要贡献(注意，这并不总是正确的)。在这种情况下，我们可以用一种非迭代的方式预先计算开始猜测。

解决方案看起来有点复杂：

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

xdat, ydat = np.loadtxt( "data.tsv", unpack=True, skiprows=1 )


def fourier(x, *as_bs):
    sum_a = 0
    sum_b = 0
    j = 1
    w = as_bs[0]
    a0 = as_bs[1]
    for i in range(2, len( as_bs ) - 1, 2 ):
        sum_a += as_bs[i] * np.cos( j * w * x )
        sum_b += as_bs[i+1] * np.sin( j * w * x )
        j = j + 1
    return a0 + sum_a + sum_b

#### initial guess
"""
This uses the fact that if y = a sin w t + b cos w t + c we have
int int y = -y/w^2 + c/2 t^2 + d t + e
i.e. we can get 1/w^2 as linear fit parameter without the danger of
a non-linear fit iterative process running into a local minimum
for details see:
https://scikit-guess.readthedocs.io/en/sine/_downloads/4b4ed1e691ff195be3ca73879a674234/Regressions-et-equations-integrales.pdf
"""

Sy = cumtrapz( ydat, xdat, initial=0 )
SSy = cumtrapz( Sy, xdat, initial=0 )

ST = np.array( [
    ydat, xdat**2, xdat, np.ones( len( xdat ) )
] )
S = np.transpose( ST )
eta = np.dot( ST, SSy )
A = np.dot( ST, S )
sol = np.linalg.solve( A, eta )

wFit = np.sqrt( -1 / sol[0] )
### linear parameters
"""
Once we have a good guess for w we can get starting guesses for
a, b and c from a standard linear fit
"""
ST = np.array( [
    np.sin( wFit * xdat ), np.cos( wFit * xdat ), np.ones( len( xdat ) )
])
S = np.transpose( ST )
eta = np.dot( ST, ydat )
A = np.dot( ST, S )
sol = np.linalg.solve( A, eta )

a1 = sol[0]
b1 = sol[1]
a0 = sol[2]

### final non-linear fit
"""
Now we can use the guesses from above as input for the final
non-linear fit. Hopefully, we are now close enough to the global minimum
and have the algorithm converge reasonably
"""
popt, pcov = curve_fit(
    fourier,
    xdat, ydat,
        p0=[
            wFit, a0, a1, b1,
            a1 / 2, b1 / 2,
            a1 / 4, b1 / 4
        ]
)
### here I assume that higher order are similar to lower orders
### but slightly smaller. ... hoping that the fit correts errors in
### this assumption

print(popt)
data_fit = fourier( xdat, *popt )

plt.plot( xdat, ydat, ls="", marker="o", ms=0.5, label="data" )
plt.plot( xdat, data_fit, label='fitting')
plt.legend()
plt.show()

两者都提供了基本相同的解决方案，后者的代码适用于更多的情况，且假设较少。