Python中4阶正弦曲线拟合的可行实现方案咨询

背景 #

Python目前没有类似MATLAB中直接调用的正弦和拟合函数（如sin4），需要手动实现4阶正弦曲线拟合，对应MATLAB中的实现逻辑如下：

ft = fittype('sin4');
opts = fitoptions('Method', 'NonlinearLeastSquares');
[fitresult, gof] = fit(xData, yData, ft, opts);

% values for [a1, b1, ..., c4] are calculated

% fitted result
x=1:length(yData);
y=a1*sin(b1*x+c1) + a2*sin(b2*x+c2) + a3*sin(b3*x+c3) + a4*sin(b4*x+c4);

已实现的1阶正弦拟合 #

已基于非线性最小二乘法完成1阶正弦拟合，代码如下：

import numpy as np
from scipy.optimize import curve_fit

def fit_sin(tt, yy):
    '''Fit sinewave to input time sequence'''
    tt = numpy.array(tt)
    yy = numpy.array(yy)
    ff = numpy.fft.fftfreq(len(tt), (tt[1]-tt[0]))
    Fyy = abs(numpy.fft.fft(yy))
    guess_freq = abs(ff[numpy.argmax(Fyy[1:])+1])
    guess_amp = numpy.std(yy) * 2.**0.5
    guess_offset = numpy.mean(yy)
    guess = numpy.array([guess_amp, 2.*numpy.pi*guess_freq, 0., guess_offset])
    def sinfunc(t, A, w, p, c):  return A * numpy.sin(w*t + p) + c
    popt, pcov = scipy.optimize.curve_fit(sinfunc, tt, yy, p0=guess)
    A, w, p, c = popt
    f = w/(2.*numpy.pi)
    fitfunc = lambda t: A * numpy.sin(w*t + p) + c
    return {"amp": A, "omega": w, "phase": p, "offset": c, "freq": f, "period": 1./f, "fitfunc": fitfunc, "maxcov": numpy.max(pcov), "rawres": (guess,popt,pcov)}

x = np.arange(0, len(yData))
res = fit_sin(x, yData)
y = res['fitfunc'](x)

当前问题 #

尝试通过迭代方式逐步提升拟合阶数（用1阶拟合结果作为2阶的初始参数，以此类推到4阶），但该方法效率极低，且拟合结果与MATLAB方案不符。

可行实现方案 #

1. 直接构造4阶正弦函数+FFT初始参数估计 #

直接定义4阶正弦拟合函数，通过FFT提取数据的主要频率分量作为初始参数，避免非线性优化陷入局部最优：

import numpy as np
from scipy.optimize import curve_fit

def fit_sin4(tt, yy):
    tt = np.array(tt)
    yy = np.array(yy)
    dt = tt[1] - tt[0]
    n = len(tt)
    
    # FFT提取前4个主要频率分量
    ff = np.fft.fftfreq(n, dt)
    fft_vals = np.abs(np.fft.fft(yy))
    # 跳过直流分量，取前4个幅值最大的频率
    sorted_indices = np.argsort(fft_vals[1:])[::-1][:4] + 1
    guess_freqs = np.abs(ff[sorted_indices])
    guess_amps = fft_vals[sorted_indices] * 2 / n  # FFT幅值转换为正弦振幅
    guess_phases = np.zeros(4)
    guess_offset = np.mean(yy)
    
    # 构造初始参数数组：[A1, w1, p1, A2, w2, p2, A3, w3, p3, A4, w4, p4, c]
    guess = []
    for amp, freq in zip(guess_amps, guess_freqs):
        guess.extend([amp, 2*np.pi*freq, 0.0])
    guess.append(guess_offset)
    
    # 定义4阶正弦拟合函数
    def sin4func(t, *params):
        c = params[-1]
        res = c
        for i in range(4):
            A = params[3*i]
            w = params[3*i+1]
            p = params[3*i+2]
            res += A * np.sin(w * t + p)
        return res
    
    # 执行非线性拟合
    popt, pcov = curve_fit(sin4func, tt, yy, p0=guess)
    
    # 构造拟合结果函数
    def fitfunc(t):
        c = popt[-1]
        res = c
        for i in range(4):
            A = popt[3*i]
            w = popt[3*i+1]
            p = popt[3*i+2]
            res += A * np.sin(w * t + p)
        return res
    
    # 整理返回结果
    result = {"offset": popt[-1], "fitfunc": fitfunc, "maxcov": np.max(pcov)}
    for i in range(4):
        idx = 3*i
        result[f"amp{i+1}"] = popt[idx]
        result[f"omega{i+1}"] = popt[idx+1]
        result[f"phase{i+1}"] = popt[idx+2]
        result[f"freq{i+1}"] = popt[idx+1]/(2*np.pi)
        result[f"period{i+1}"] = 1/result[f"freq{i+1}"]
    
    return result

# 使用示例
x = np.arange(0, len(yData))
res = fit_sin4(x, yData)
y_fit = res['fitfunc'](x)

2. 添加参数约束提升拟合稳定性 #

如果拟合结果仍不理想，可通过以下方式优化：

• 限制频率为正：通过curve_fit的bounds参数设置各频率参数的最小值为0
• 固定谐波关系：若已知频率为基频的整数倍，可直接约束对应参数
• 增加迭代次数：调整curve_fit的maxfev参数，确保优化过程收敛

3. 使用lmfit库简化复杂拟合 #

lmfit库提供更灵活的参数管理和拟合控制，适合多参数非线性拟合场景：

import numpy as np
from lmfit import Model, Parameters

def sin_term(t, A, w, p):
    return A * np.sin(w*t + p)

def fit_sin4_lmfit(tt, yy):
    tt = np.array(tt)
    yy = np.array(yy)
    dt = tt[1]-tt[0]
    n = len(tt)
    
    # FFT估计初始参数
    ff = np.fft.fftfreq(n, dt)
    fft_vals = np.abs(np.fft.fft(yy))
    sorted_indices = np.argsort(fft_vals[1:])[::-1][:4]+1
    guess_freqs = np.abs(ff[sorted_indices])
    guess_amps = fft_vals[sorted_indices]*2/n
    guess_offset = np.mean(yy)
    
    # 创建4阶正弦模型
    model = Model(lambda t, c, A1, w1, p1, A2, w2, p2, A3, w3, p3, A4, w4, p4: 
                  c + sin_term(t,A1,w1,p1) + sin_term(t,A2,w2,p2) + sin_term(t,A3,w3,p3) + sin_term(t,A4,w4,p4))
    
    # 设置参数及约束
    params = Parameters()
    params.add('c', value=guess_offset)
    for i in range(4):
        params.add(f'A{i+1}', value=guess_amps[i])
        params.add(f'w{i+1}', value=2*np.pi*guess_freqs[i], min=0)  # 强制频率为正
        params.add(f'p{i+1}', value=0.0, min=-np.pi, max=np.pi)  # 相位约束在[-π, π]
    
    # 执行拟合
    result = model.fit(yy, params, t=tt)
    
    # 构造拟合函数
    def fitfunc(t):
        return result.eval(t=t)
    
    return {"fitfunc": fitfunc, "params": result.params, "gof": result.goodness_of_fit}

# 使用示例
res_lmfit = fit_sin4_lmfit(x, yData)
y_fit_lmfit = res_lmfit['fitfunc'](x)

内容的提问来源于Stack Exchange，提问作者Mirka