Matlab傅里叶级数拟合良好，Python curve_fit实现效果不佳求助

我在Matlab 2016b的Curve Fitting Tool中用4项傅里叶级数可以完美拟合观测信号，但用Python代码实现时效果很差，代码如下：

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

def fourier(x,*a):
# maybe a = [ a0 , w , a1 , b1 , a2 , b2 ....]
    ret = a[0]
    w   = a[1]
    deg = int(len(a)/2-1)
    
    for i in range(1,deg+1):
        ret += a[2*i] * np.cos(i*w*x) + a[2*i+1] * np.sin(i*w*x)
    return ret

phase1_gtc202301 = np.loadtxt('XSM-20230111233742-20230209154234-0.1Hz.txt')[:,11]

dt = 1/0.1/60/60/24                         # DAY
N  = 247679
t  = np.linspace(dt,N*dt,N)

[A,B] = curve_fit(fourier,t,phase1_gtc202301,[0.2]*8,absolute_sigma=True)
phasefit = fourier(t,*A)

plt.plot(t,phase1_gtc202301,t,phasefit)

拟合效果对比：

问题分析与解决 #

1. 傅里叶级数系数索引错误

你的fourier函数中参数索引逻辑错误，当i从1到4时，2*i和2*i+1会超出8个初始参数的合理范围（初始参数应为a0, w, a1, b1, a2, b2, a3, b3），导致错误调用内存数据。

修正后的函数：

def fourier(x, *a):
    ret = a[0]  # 直流分量a0
    w = a[1]    # 基频w
    deg = int((len(a) - 2) / 2)  # 傅里叶项数
    
    for i in range(1, deg + 1):
        # 第i次谐波的余弦/正弦系数对应a[2*i-1]和a[2*i]
        ret += a[2*i - 1] * np.cos(i * w * x) + a[2*i] * np.sin(i * w * x)
    return ret

2. 初始参数设置不合理

[0.2]*8的初始值完全没有针对性，尤其是基频w偏离实际值太远，会导致curve_fit无法收敛到最优解。建议先观察信号周期，估算基频：

• 假设信号周期为T（天），则基频w=2*np.pi/T
• 其他系数初始值设为小数值或0

示例初始参数（假设周期为1天）：

initial_guess = [0.0, 2*np.pi, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]

3. 拟合迭代次数不足

默认迭代次数可能不够，需手动增加maxfev参数确保收敛：

[A,B] = curve_fit(fourier, t, phase1_gtc202301, initial_guess, absolute_sigma=True, maxfev=10000)

完整修正代码

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

def fourier(x, *a):
    ret = a[0]
    w = a[1]
    deg = int((len(a) - 2) / 2)
    
    for i in range(1, deg + 1):
        ret += a[2*i - 1] * np.cos(i * w * x) + a[2*i] * np.sin(i * w * x)
    return ret

# 加载数据
phase1_gtc202301 = np.loadtxt('XSM-20230111233742-20230209154234-0.1Hz.txt')[:,11]

# 时间轴计算（直接用数据长度避免手动输入错误）
dt = 1/(0.1*60*60*24)  # 转换为天
N = len(phase1_gtc202301)
t = np.linspace(dt, N*dt, N)

# 初始参数（根据信号周期调整基频）
initial_guess = [0.0, 2*np.pi, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
# 拟合，增加迭代次数
[A,B] = curve_fit(fourier, t, phase1_gtc202301, initial_guess, absolute_sigma=True, maxfev=10000)

phasefit = fourier(t, *A)

# 绘图对比
plt.figure(figsize=(12,6))
plt.plot(t, phase1_gtc202301, label='原始信号', alpha=0.7)
plt.plot(t, phasefit, label='拟合信号', linewidth=2)
plt.legend()
plt.xlabel('时间（天）')
plt.ylabel('相位')
plt.show()

补充说明 #

Matlab的Curve Fitting Tool会自动处理傅里叶拟合的参数索引和初始值估计，因此更容易收敛；而Python的curve_fit高度依赖合理的初始参数和正确的函数定义，否则容易陷入局部最优或不收敛。

内容的提问来源于stack exchange，提问作者EverChyi