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

3vpjnl9f 于 7个月前发布在其他

关注(0)|答案(2)|浏览(96)

背景

我一直在经历一个article，它解释了人类“心电图”波的数学重建，具体来说，“QRS”波。用于重建“QRS”波的函数是以下形式的四次多项式：*f（t）= a（t - B）*4 + h， 其中“a”是系数，“B”是半周期，“h”是幅度。
现在，我想得到的傅里叶级数是 f（t）= -0.0000156（t − 20）<$+ 2.5。并且该函数应每1000 ms重复一次。
曲线应该从（0，0）开始，由于脉冲是40毫秒长，它将通过（40，0），并且以 t = 20为中心。
现在，为了计算傅立叶系数，并绘制傅立叶级数，我在Python中编写了以下程序，在详细的article之后，
程序成功运行，但不会每1000 ms产生一次“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()

scipy

来源：https://stackoverflow.com/questions/76900684/how-to-obtain-a-fourier-series-for-a-custom-function-that-generates-the-periodic

2条答案

按热度按时间

w9apscun1#

需要改变两个部分，以生成所需的ECG信号，并且两者都围绕以下事实展开：您不希望仅重复ecgFunc抛物线，而是在1000 ms的剩余时间内重复基线上的光点。
因此，需要修改函数：

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上循环似乎更有趣，因为你可以很快看到高阶项如何越来越好地逼近信号。
关于您在评论中的后续问题：如果你想保持轴固定，你不应该不断地改变它的极限

x_l_plot = x_l_plot + step_size
        x_u_plot = x_u_plot + step_size

此外，更新数字变得越来越慢。我猜在每一次迭代中，它会在现有的那些上面画一个图，把所有的都保存在内存中。因此，我建议将plt.clf()引入循环：

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()
        plt.clf()

或者，您可以按照其他人的建议查看动画API。
可视化不同数量的n我基于你的原始代码实现如下：

if __name__ == "__main__":

    # plt.style.use('dark_background')
    plt.style.use('seaborn')

    # Limits for the functions
    li = 0
    lf = 1000
    l = 500  # (lf-li)/2.0

    for n in range(2, 30):

        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 = 1

        # 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 = []
        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(-1, 3)

        plt.plot(x, y1, c='darkkhaki', label='ecg Wave')
        plt.plot(x, y1_fourier, c='forestgreen', label='Fourier Approximation')
        plt.legend()
        plt.xlim([0, 2000])
        plt.draw()

        plt.pause(0.1)
        plt.clf()

这里，如所建议的，ecg信号的长度保持恒定，而n的数量在每次迭代中增加。
希望能帮上忙。

赞(0）回复(0）举报 8个月前

30byixjq2#

感谢@Flow的帮助，也感谢@D.L建议使用matplotlib.animation来处理傅立叶级数。
我浏览了matplotlib.animation的文档，并在此提供了心电图傅立叶级数近似的工作动画。

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=96)

ax = fig.add_subplot(1, 1, 1)
ax.set_facecolor((0.945, 0.831, 0.886))

 
 
 
# 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 = 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 = 1000
    l = 15#(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 = []
    
    y_plot1 = []
    y_plot1_fourier = []
    
    x_l_plot = x_l-10
    x_u_plot = x_l_plot + 300
    
    ax.set_xlim([x_l_plot, x_u_plot])
    ax.set_ylim([-3, 3])
    
    #ecg wave
    plt.grid(visible=True, which='major', color='k', linestyle='-', alpha=0.3)
    plt.grid(visible=True, which='minor', color='r', linestyle='-', alpha=0.2)
    plt.minorticks_on()
    
    #plotting the graph
    im, = ax.plot([], [], c="black", linewidth=0.5)
    
    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])
    print(x.size)
    
    x_plot = np.array(x_plot)
    y_plot1_fourier=np.array(y_plot1_fourier)
    

    
    
    def func(n):
        
        im.set_xdata(x[0:n])
        print(x[n])
        
        im.set_ydata(y_plot1_fourier[0:n])
        #print(y_plot1_fourier[n])
        
        if x[n] > x_u_plot:
            
            lim = ax.set_xlim(x[n]-x_u_plot, x[n])
            
        else:
            
            lim = ax.set_xlim(x_l_plot, x_u_plot)
            
        return im
    
    
    ani = animation.FuncAnimation(fig, func, frames = x.size, interval = 100, blit=False)
    
    plt.show()
    # 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.grid(visible=True, which='major', color='k', linestyle='-', alpha=0.3)
    #     plt.grid(visible=True, which='minor', color='r', linestyle='-', alpha=0.2)
    #     plt.minorticks_on()
    #     plt.plot(x_plot,y_plot1_fourier,c='black',label='Fourier Approximation', linewidth=0.7)

        
            
    #     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()