- 积分
- 1965
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2014-4-13
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 柿子柿子柿子 于 2024-1-7 09:08 编辑
原本在想着怎么通过傅里叶变换找有序数列的周期、振幅和相位的,找了一圈没找到网上现成的程序,但是通过一堆零零碎碎的教程给我指向了《高数》……于是跟着《高数》的内容调试出了这个段拟合代码
【吐槽】想往这个发帖的框框里填代码是真的麻烦……
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
def fourier(x, *a):
aa = a[1:]
m = len(aa) // 2
ret = a[0]
for i,nn in enumerate(aa[:m],1):
ret += nn * np.cos(i*x)
for i,nn in enumerate(aa[m:],1):
ret += nn * np.sin(i*x)
return ret
def fit_fourier(x,y,n):
popt, pcov = curve_fit(fourier, x, y, [1.]*(2*n+1))
p = fourier(x,*popt)
# print(popt)
return popt
def test_func(x):
y = 1.5 + 2.5 * np.sin(x - np.pi/4)
# 添加噪声
y += 0.2 * np.random.normal()
return y
if __name__ == '__main__':
op = Path("pic")
op.mkdir(511,True,True)
x = np.linspace(0,2*np.pi,20)
y = test_func(x)
a = fit_fourier(x,y,2)
x1 = np.linspace(0,4*np.pi,50)
y1 = fourier(x1, *a)
plt.plot(x,y,'ko',label='data')
plt.plot(x1,y1,'r',label='fit')
plt.legend()
plt.savefig(op.joinpath("fitting.png"), dpi=300)
1.5+2.5sin(x-pi/4)
|
|