1

我试图找到一个适合这个情节的傅立叶级数： 在此处输入图像描述

首先我使用了这个网站，我发现了以下傅立叶系数

a = [0.983, 0.292, 0.110, 0.006, -0.022, -0.018, 0.006]
b = [0.000, 0.246, 0.002, -0.021, 0.017, 0.046, 0.019]
omega = 0.045

产生了以下结果 在此处输入图像描述

这很好，但后来我决定尝试在 python 中编写我自己的曲线拟合器，并发现库 symfit 似乎是最好的选择。我什至重新利用了文档中的一个示例，但几乎没有改变。

果然，当我运行代码时，直接使用 fit 对象制作的图表中的 fit 似乎非常好：

plt.plot(xdata, fit.model(x=xdata, **fit_result.params).y, color='green', ls=':')

在此处输入图像描述

但后来我决定直接使用傅立叶系数，而不是依赖拟合对象来尝试重现函数。这就是代码打印这些系数的方式：

Parameter Value        Standard Deviation
a0        -1.206104e+03 nan
a1        1.799223e+03 1.087713e+05
a10       2.393892e+00 8.023119e+02
a2        -8.465307e+02 2.133670e+05
a3        6.798590e+02 nan
a4        -9.287399e+02 nan
a5        7.388156e+02 nan
a6        -2.099387e+02 1.376980e+05
a7        -9.838123e+01 9.799211e+04
a8        9.821148e+01 3.739932e+04
a9        -2.827002e+01 8.152094e+03
b1        -6.084419e+02 2.480938e+05
b10       4.586827e+01 3.396681e+02
b2        9.311030e+02 2.178771e+05
b3        -8.877734e+02 2.406878e+05
b4        8.810145e+02 4.462669e+05
b5        -1.288586e+03 4.376938e+05
b6        1.734395e+03 2.730319e+05
b7        -1.583775e+03 1.137506e+05
b8        9.115344e+02 3.084382e+04
b9        -3.045882e+02 4.905459e+03
w         -1.370387e-02 nan

标准差是一个巨大的危险信号，它们大得离谱，果然，这就是结果

在此处输入图像描述

这似乎是一个糟糕的笑话。

这是使用 10 个系数完成的，如果我只使用 6 个，我得到的东西几乎不能接受。

在此处输入图像描述

最糟糕的是，直接用 fit 对象制作的图表必须意味着代码确实产生了很好的拟合，但由于某种原因，它给我的傅立叶系数只是谎言。

这是我的代码。

首先，这就是我适合功能的方式

from symfit import parameters, variables, sin, cos, Fit
import numpy as np
import matplotlib.pyplot as plt
import math as m

def fourier_series(x, f, n=0):
    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)]))
    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=10)}
print(model_dict)

#This is the graph I'm trying to fit turned into data points
xlist=[381.854,385.274,387.55,388.973,390.656,392.671,394.281,396.37,398.014,402.047,408.288,411.885,414.452,418.099,\
       423.493,426.901,428.014,430.068,433.114,438.699,445.023,448.562,453.307,458.014,463.663,467.877,474.536,\
       479.384,485.409,489.247,496.283,501.164,508.191,514.315,520.1]
ylist=[0.02816,0.09827,0.18838,0.2948,0.36316,0.46243,0.53795,0.60116,0.64162,0.68361,0.71676,0.79285,0.87861,0.94579,\
       0.99422,0.95307,0.89017,0.85549,0.79285,0.76301,0.7273,0.65896,0.59621,0.53757,0.47969,0.43353,0.38501,\
       0.33526,0.28305,0.23121,0.18838,0.14451,0.11555,0.08092,0.07185]

xdata = np.array(xlist)
ydata = np.array(ylist)

# 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=':')

这就是我采用这些系数并尝试重现图表的方式。我有不同的 ab 和 omega 数据集来尝试重现其他傅立叶级数。阶梯楔（文档中的示例）是唯一有效的。

from symfit import parameters, variables, sin, cos, Fit
import numpy as np
import matplotlib.pyplot as plt
import math as m
from random import random

def f(x,a,b,omega):
    output = a[0]/2
    for k in range(1,len(a)-1):
        output = output + a[k]*m.cos(k*omega*x) + b[k]*m.sin(k*omega*x)
    return output

def monteCarlo(min,max,a,b,omega):
    x = min + random()*(max-min)
    return (x,f(x,a,b,omega))

xdata=[]
ydata=[]

#for step function between -3 and 3
#a=[5.000000e-01,1.106075e-10,2.364473e-11,1.697959e-10,2.822410e-11,3.860925e-11,-1.522386e-10]
#b=[0,6.267589e-01,2.020945e-02,1.846406e-01,3.623074e-02,8.726419e-02,4.518721e-02]
#omega=8.615115e-01

#for our spectrum with n=6
a = [-1.755069e+02,-1.097051e+02,1.571135e+02,6.491389e+01,1.220463e+02,1.778165e+02,4.564577e+01]
b = [0,-5.523502e+01,1.517397e+02,1.325198e+02,-1.015198e+02,-3.255781e+01,2.495038e+01]
omega = -5.353751e-03

#for our spectrum with n=10
#a = [-1.206104e+03,1.799223e+03,-8.465307e+02,6.798590e+02,-9.287399e+02,7.388156e+02,-2.099387e+02,-9.838123e+01,\
#    9.821148e+01,-2.827002e+01,2.393892e+00]
#b = [0,-6.084419e+02,9.311030e+02,-8.877734e+02,8.810145e+02,-1.288586e+03,1.734395e+03,-1.583775e+03,9.115344e+02,\
#     -3.045882e+02,4.586827e+01]
#omega = -1.370387e-02

#a = [0.983, 0.292, 0.110, 0.006, -0.022, -0.018, 0.006]
#b = [0.000, 0.246, 0.002, -0.021, 0.017, 0.046, 0.019]
#omega = 0.045

for n in range(0,1000):
    coor = monteCarlo(375,525,a,b,omega)
    xdata.append(coor[0])
    ydata.append(coor[1])

plt.plot(xdata, ydata, 'bo')
4

0 回答 0