如何为周期1000ms的自定义QRS心电函数求傅里叶级数及排错

问题背景 #

用四次多项式f(t) = -0.0000156(t − 20)⁴ + 2.5重构QRS波，要求该函数每1000ms重复一次：脉冲时长40ms（0-40ms有波形，40-1000ms为0），以t=20为中心，起点和终点均为(0,0)、(40,0)。编写Python代码计算傅里叶级数并绘制，但无法按1000ms周期生成波峰。

核心错误分析 #

1. 周期参数l的硬编码错误
   傅里叶级数计算中，l应为周期的一半（周期T=1000ms，故l=T/2=500），但代码中错误地将l固定为40，导致基波频率完全偏离，无法匹配1000ms的周期。
2. 参数不统一
   fourierCoeffs和fourierSeries函数中的l必须保持一致，且严格对应周期的一半，否则傅里叶级数的合成结果会完全错误。

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt 
import scipy.integrate as integrate
 
fig = plt.figure(figsize=(7, 7), dpi=120)

# 生成周期函数的工具函数
def periodicf(li,lf,f,x):
    if x>=li and x<=lf :
        return f(x)
    elif x>lf:
        x_new=x-(lf-li)
        return periodicf(li,lf,f,x_new)
    elif x<(li):
        x_new=x+(lf-li)
        return periodicf(li,lf,f,x_new)

def ecgFuncP(li,lf,x):
    return periodicf(li,lf,ecgFunc,x)

# QRS波原函数：0-40ms有波形，其余为0
def ecgFunc(x):
    if 0 <= x < 40:
        return -0.0000156 * (x - 20) ** 4 + 2.5
    else:
        return 0

# 计算傅里叶系数
def fourierCoeffs(li, lf, n, f):
    T = lf - li  # 周期
    l = T / 2    # 周期的一半
    # 常数项
    a0 = 1/l * integrate.quad(lambda x: f(x), li, lf)[0]
    A = np.zeros((n))  # 余弦系数
    B = np.zeros((n))  # 正弦系数
     
    for i in range(1, n+1):
        A[i-1] = 1/l * integrate.quad(lambda x: f(x)*np.cos(i*np.pi*x/l), li, lf)[0]
        B[i-1] = 1/l * integrate.quad(lambda x: f(x)*np.sin(i*np.pi*x/l), li, lf)[0]
 
    return [a0/2.0, A, B]

# 傅里叶级数合成
def fourierSeries(coeffs,x,l,n):
    value = coeffs[0]
    for i in range(1,n+1):
        value += coeffs[1][i-1] * np.cos(i*np.pi*x/l) + coeffs[2][i-1] * np.sin(i*np.pi*x/l)
    return value
      
if __name__ == "__main__":
    plt.style.use('seaborn')
     
    # 周期区间：0-1000ms
    li = 0
    lf = 1000
    T = lf - li
    l = T / 2.0
 
    # 谐波项数量
    n = 100
 
    plt.title(f'Fourier Series Approximation\nECG Wave\n n = {n}')
        
    # 计算傅里叶系数
    coeffsEcgFunc = fourierCoeffs(li, lf, n, ecgFunc)

    # 绘图参数
    step_size = 0.5
    x_l = 0
    x_u = 4000  # 绘制4个周期
    x = np.arange(x_l, x_u, step_size)
    
    # 生成原周期函数和傅里叶近似值
    y1 = [ecgFuncP(li, lf, xi) for xi in x]
    y1_fourier = [fourierSeries(coeffsEcgFunc, xi, l, n) for xi in x]

    # 动画绘图设置
    x_plot = []
    y_plot1_fourier = []
    x_l_plot = x_l - 13
    x_u_plot = x_l_plot + 300
    plt.xlim(x_l_plot, x_u_plot)
    plt.ylim(-3, 3)
    
    for i in range(x.size):
        x_plot.append(x[i])
        y_plot1_fourier.append(y1_fourier[i])

        plt.plot(x_plot, y_plot1_fourier, c='forestgreen', label='Fourier Approximation')
        
        x_l_plot += step_size
        x_u_plot += step_size
        plt.xlim(x_l_plot, x_u_plot)
        plt.pause(0.001)
        if i == 0:
            plt.legend()

    plt.show()

修正说明 #

• 将l的赋值改为l = (lf - li)/2，对应1000ms周期的一半（500），确保傅里叶级数的基波频率正确。
• 优化了原函数ecgFunc的条件判断，明确0<=x<40的范围，避免边界歧义。
• 统一了所有涉及l的计算，保证系数计算与级数合成使用相同的周期参数。

内容的提问来源于stack exchange，提问作者Vijith Kumar V