傅里叶级数拟合 Python
Posted
技术标签:
【中文标题】傅里叶级数拟合 Python【英文标题】：Fourier Series Fit in Python 【发布时间】：2019-03-02 15:44:55 【问题描述】：我有一些数据想要使用 2 阶、3 阶或 4 阶的傅立叶级数来拟合。
虽然this 堆栈溢出问题和答案接近我想要使用 scipy 做的事情，但他们已经将系数预先定义为 tau = 0.045。我希望我能找到具有 95% 置信区间的可能系数（a0、w1、w2、w3 等），就像傅立叶级数的 MATLAB curve fit 等效项一样。我看到的另一个选项是使用fourier_series from sympy，但是这个函数只适用于适合定义函数的符号参数，而不是原始数据。
1) sympyfourier_series 有没有办法接收原始数据，而不是使用这个库的函数或其他解决方法？
2) 考虑到有多个未知数（系数），或者对数据进行 scipy 曲线拟合
【问题讨论】：
"我有一些数据" => 你不需要 SymPy。这不是为了处理数字。 是什么阻止您拟合基频和系数？curve_fit 不关心函数的形式，它只找到一组输入，使与一组预期输出的差异最小化。
我的开源Python在线曲线拟合器有傅里叶级数，你可以在线尝试，如果他们做了你需要的，每页底部都有源代码链接。 3 期链接：zunzun.com/Equation/2/FourierSeries/3%20Term%20Standard 和 4 期链接：zunzun.com/Equation/2/FourierSeries/4%20Term%20Standard
【参考方案1】：
这是我最近使用的一个简单代码。我知道频率，所以稍微修改一下以适应它。不过，最初的猜测必须非常好（我认为）。 另一方面，有几种获得基频的好方法。
import numpy as np
from scipy.optimize import leastsq
def make_sine_graph( params, xData):
"""
take amplitudes A and phases P in form [ A0, A1, A2, ..., An, P0, P1,..., Pn ]
and construct function f = A0 sin( w t + P0) + A1 sin( 2 w t + Pn ) + ... + An sin( n w t + Pn )
and return f( x )
"""
fr = params[0]
npara = params[1:]
lp =len( npara )
amps = npara[ : lp // 2 ]
phases = npara[ lp // 2 : ]
fact = range(1, lp // 2 + 1 )
return [ sum( [ a * np.sin( 2 * np.pi * x * f * fr + p ) for a, p, f in zip( amps, phases, fact ) ] ) for x in xData ]
def sine_residuals( params , xData, yData):
yTh = make_sine_graph( params, xData )
diff = [ y - yt for y, yt in zip( yData, yTh ) ]
return diff
def sine_fit_graph( xData, yData, freqGuess=100., sineorder = 3 ):
aStart = sineorder * [ 0 ]
aStart[0] = max( yData )
pStart = sineorder * [ 0 ]
result, _ = leastsq( sine_residuals, [ freqGuess ] + aStart + pStart, args=( xData, yData ) )
return result
if __name__ == '__main__':
import matplotlib.pyplot as plt
timeList = np.linspace( 0, .1, 777 )
signalList = make_sine_graph( [ 113.7 ] + [1,.5,-.3,0,.1, 0,.01,.02,.03,.04], timeList )
result = sine_fit_graph( timeList, signalList, freqGuess=110., sineorder = 3 )
print result
fitList = make_sine_graph( result, timeList )
fig = plt.figure()
ax = fig.add_subplot( 1, 1 ,1 )
ax.plot( timeList, signalList )
ax.plot( timeList, fitList, '--' )
plt.show()
提供
<< [ 1.13699742e+02 9.99722859e-01 -5.00511764e-01 3.00772260e-01
1.04248878e-03 -3.13050074e+00 -3.12358208e+00 ]
【讨论】：
【参考方案2】：我认为使用this API 在 python 脚本上调用 MATLAB 函数来完成所有繁重的工作可能会更容易，而不是处理 sympy 和 scipy 的所有细节。
我的解决方案如下：
-
安装matlab及拟合函数工具箱
通过在 extern/engines/python 下的 matlab 根文件夹中运行
python setup.py install 来安装适用于 python 的 matlab 引擎。注意：它仅适用于 python 2.7、3.5 和 3.6
使用以下代码：
import matlab.engine
import numpy as np
eng = matlab.engine.start_matlab()
eng.evalc("idx = transpose(1:)".format(len(diff)))
eng.workspace['y'] = eng.transpose(eng.cell2mat(diff.tolist()))
eng.evalc("f = fit(idx, y, 'fourier3')")
y_f = eng.evalc("f(idx)").replace('ans =', '')
y_f = np.fromstring(y_f, dtype=float, sep='\n')
几点说明：
eng.workspace['myVariable'] 是使用 python 结果声明 matlab 变量，以后可以通过 evalc 调用
eng.evalc 以 'ans = ...' 形式返回一个字符串
这段代码中的diff 只是一些数据与其最小二乘线之间的差异，它属于系列类型。
Python 列表等效于 MATLAB 中的单元格类型。
【讨论】：
【参考方案3】：如果愿意，您可以非常接近 sympy 代码进行数据拟合，使用我为此目的编写的名为 symfit 的包。它基本上使用sympy 接口包装scipy。使用symfit，您可以执行以下操作：
from symfit import parameters, variables, sin, cos, Fit
import numpy as np
import matplotlib.pyplot as plt
def fourier_series(x, f, n=0):
"""
Returns a symbolic fourier series of order `n`.
:param n: Order of the fourier series.
:param x: Independent variable
:param f: Frequency of the fourier series
"""
# Make the parameter objects for all the terms
a0, *cos_a = parameters(','.join(['a'.format(i) for i in range(0, n + 1)]))
sin_b = parameters(','.join(['b'.format(i) for i in range(1, n + 1)]))
# Construct the series
series = a0 + sum(ai * cos(i * f * x) + bi * sin(i * f * x)
for i, (ai, bi) in enumerate(zip(cos_a, sin_b), start=1))
return series
x, y = variables('x, y')
w, = parameters('w')
model_dict = y: fourier_series(x, f=w, n=3)
print(model_dict)
这将打印我们想要的符号模型：
y: a0 + a1*cos(w*x) + a2*cos(2*w*x) + a3*cos(3*w*x) + b1*sin(w*x) + b2*sin(2*w*x) + b3*sin(3*w*x)
接下来，我将把它拟合成一个简单的步进函数，向你展示它是如何工作的：
# Make step function data
xdata = np.linspace(-np.pi, np.pi)
ydata = np.zeros_like(xdata)
ydata[xdata > 0] = 1
# Define a Fit object for this model and data
fit = Fit(model_dict, x=xdata, y=ydata)
fit_result = fit.execute()
print(fit_result)
# Plot the result
plt.plot(xdata, ydata)
plt.plot(xdata, fit.model(x=xdata, **fit_result.params).y, color='green', ls=':')
这将打印：
Parameter Value Standard Deviation
a0 5.000000e-01 2.075395e-02
a1 -4.903805e-12 3.277426e-02
a2 5.325068e-12 3.197889e-02
a3 -4.857033e-12 3.080979e-02
b1 6.267589e-01 2.546980e-02
b2 1.986491e-02 2.637273e-02
b3 1.846406e-01 2.725019e-02
w 8.671471e-01 3.132108e-02
Fitting status message: Optimization terminated successfully.
Number of iterations: 44
Regression Coefficient: 0.9401712713086535
并产生以下情节：
就这么简单！我把剩下的留给你想象。欲了解更多信息，您可以找到documentation here。
【讨论】：
我喜欢接近 OP 的简单方法，但想提一下，使用更多的'non-linear' version，即使用阶段而不是sin-cos 组合有时具有更好的收敛性。
【参考方案4】：
傅立叶曲线拟合具有封闭形式的解。这个函数可以帮你计算出来。
def fourier_curve_fit(ser, no_fourier=3, display_latex=True, series=False):
"""
Apply fourier curve fitting to series.
ser: pandas.Series
Contains data stored in Series.
no_fourier: int
degree of fourier series to be used.
Returns
pandas.Series
"""
from IPython.display import display, Math
import numpy as np
import pandas as pd
ser_ = ser.reset_index(drop=True)
x = ser_.index
y = ser_.values
A = (2/len(y))*np.matmul(y.reshape(1,-1),
np.array([np.cos(2*np.pi/len(y)*n*np.arange(len(y))) for n in np.arange(no_fourier+1)]).T
).flatten()
B = (2/len(y))*np.matmul(y.reshape(1,-1),
np.array([np.sin(2*np.pi/len(y)*n*np.arange(len(y))) for n in np.arange(no_fourier+1)]).T
).flatten()
L = max(x)-min(x)
if display_latex:
Omega = '\\frac2*\\pi'+str(L)+''
fourier_equation = ''
for i,(a,b) in enumerate(zip(A,B)):
if i==0:
fourier_equation += f'f(x)=a:.4f \\\\'
else:
fourier_equation += f' a:+.4f*\\cos(Omega*i*x) b:+.4f*\\sin(Omega*i*x) \\\\'
display(Math(fourier_equation))
coeff = dict([(f'ai',j) for i,j in enumerate(A)]+[(f'bi',j) for i,j in enumerate(B)])
if series:
As = np.matmul(A.reshape(1,-1), np.array([np.cos(2*np.pi/len(y)*n*np.arange(len(y))) for n in np.arange(no_fourier+1)])).flatten()
Bs = np.matmul(B.reshape(1,-1), np.array([np.sin(2*np.pi/len(y)*n*np.arange(len(y))) for n in np.arange(no_fourier+1)])).flatten()
return pd.Series(data=(As+Bs), index=ser.index, name='Fourier Fitted')-coeff['a0']/2
else:
return pd.Series(data=coeff, name='Fourier Coefficients')
【讨论】：
【参考方案5】：我可以通过将我的猜测作为 *args 传递并在函数中检索所需的参数来使用 scipy 做到这一点：
def fourier(x, *params):
params = np.array(params).reshape(-1,3)
a = params[:, 0]
b = params[:, 1]
c = params[:, 2]
ret = a[0] * np.sin(np.pi / b[0] * x) + c[0]
for deg in range(1, len(a)):
ret += a[deg] * np.sin((deg+1) * np.pi / b[deg] * x) + c[0]
return ret
# num fourier terms
terms = 20
params, covariance = curve_fit(
fourier, x_data, y_data, maxfev=10000,
# Initial guesses for a, b, c
p0=[0.3, 1.0, 0.5]*terms
)
【讨论】：
以上是关于傅里叶级数拟合 Python的主要内容，如果未能解决你的问题，请参考以下文章