傅里叶三角级数拟合

柿子柿子柿子

原本在想着怎么通过傅里叶变换找有序数列的周期、振幅和相位的，找了一圈没找到网上现成的程序，但是通过一堆零零碎碎的教程给我指向了《高数》……于是跟着《高数》的内容调试出了这个段拟合代码
【吐槽】想往这个发帖的框框里填代码是真的麻烦……


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)