我想获得自定义函数的傅立叶级数，该函数每 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
 
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):
    return -0.0000156*(x - 20)** 4 + 2.5
 

def fourierCoeffs(li, lf, n, f):
    l = 500#(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 = 40
    l = 1000#(lf-li)/2.0
 
    # Number of harmonic terms
    n = 1
    for n in range(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.05
 
        # 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]
        print(x)
        print(y1_fourier)
 
        x_plot =[]
        # Sawtooth
        y_plot1 = []
        y_plot1_fourier = []
     
        x_l_plot = x_l - 13
        x_u_plot = x_l_plot + 14
        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,c='darkkhaki',label='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()
 
         
        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 似乎更有趣，因为您可以快速看到高阶项如何越来越好地逼近信号。