如何获取自定义函数的傅里叶级数，该函数每 1000 毫秒生成“心电图”波的周期性“QRS”段？

Question

背景

我一直在阅读一篇文章，它解释了人类“ECG”波（特别是“QRS”波）的数学重建。用于重建“QRS”波的函数是以下形式的四次多项式： f(t) = a(t - b)**4 + h, 其中“a”是系数，“ b" 是半周期，"h" 是幅度。

现在，我想要获得傅立叶级数的函数是 f(t) = -0.0000156(t − 20)⁴ + 2.5。并且该函数应该每 1000 毫秒重复一次。

曲线应从 (0,0) 开始，由于脉冲长 40 ms，因此它将穿过 (40,0)，并以 t = 20 为中心。

现在，为了计算傅里叶系数并绘制傅里叶级数，在遵循详细的文章

之后，我用Python编写了以下程序

程序运行成功，但没有每 1000 毫秒生成一次“QRS”信号。看来我没能找出错误在哪里。我会请求 StackOverflow 社区查看这个问题并帮助我解决问题。

这是我编写的用于生成函数的傅立叶近似的代码：

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

 
# Function that will convert any given function 'f' defined in a given range '[li,lf]' to a periodic function of period 'lf-li' 
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)
 

def ecgFunc(x):
    if x < 40:
        return -0.0000156 * (x - 20) ** 4 + 2.5
    else:
        return 0
 

def fourierCoeffs(li, lf, n, f):
    l = 40#(lf-li)/2
    # Constant term
    a0=1/l*integrate.quad(lambda x: f(x), li, lf)[0]
    # Cosine coefficents
    A = np.zeros((n))
    # Sine coefficents
    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]
 
# This functions returns the value of the Fourier series for a given value of x given the already calculated Fourier coefficients
def fourierSeries(coeffs,x,l,n):
    
    value = coeffs[0]
    for i in range(1,n+1):
        value = 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('dark_background')
    plt.style.use('seaborn')
     
 
    # Limits for the functions
    li = 0
    lf = 1000
    l = 40#(lf-li)/2.0
 
    # Number of harmonic terms
    
    n = 100
 
    plt.title('Fourier Series Approximation\nECG Wave\n n = '+str(n))
        
    # Fourier coeffficients for various functions
    coeffsEcgFunc = fourierCoeffs(li,lf,n,ecgFunc)


    # Step size for plotting
    step_size = 0.5

    # Limits for plotting
    x_l = 0
    x_u = 4000

    # Sample values of x for plotting
    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 = []
    # Sawtooth
    y_plot1 = []
    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])
        # Actual function values
        y_plot1.append(y1[i])

        # Values from fourier series
        y_plot1_fourier.append(y1_fourier[i])


    #     #ecg wave

        plt.plot(x_plot,y_plot1_fourier,c='forestgreen',label='Fourier Approximation')

        
            
        x_l_plot = x_l_plot + step_size
        x_u_plot = x_u_plot + step_size
        plt.xlim(x_l_plot,x_u_plot)
        plt.pause(0.001)
        if i==0:
            plt.legend()
    print(x_plot)
    print(y_plot1_fourier)
    plt.clf()


plt.show()

Answer 1

为了生成您想要的心电图信号，需要更改两个部分，并且这两个部分都围绕这样一个事实：您不希望只重复 ecgFunc 抛物线，而是希望在 1000 毫秒的剩余时间内在基线顶部出现该尖峰。

因此需要修改功能：

def ecgFunc(x):
    if x < 40:
        return -0.0000156 * (x - 20) ** 4 + 2.5
    else:
        return 0

稍后您可以设置重复该功能的限制。当然，您想重复 0 到 1000 之间的所有内容：

if __name__ == "__main__":
    (...)
    # Limits for the functions
    li = 0
    lf = 1000 # <-- changed from 40
    l = 500 #(lf-li)/2.0 <-- 500 seemed to be correct in my tests.

对你的循环的一些评论： 循环谐波项的数量和循环数据点将变得非常冗长。我建议仅循环这些参数中的一个，对我来说，循环 n 似乎更有趣，因为您可以快速看到高阶项如何越来越好地逼近信号。