Python - SciPy - ECG信号的谱分析及数字滤波

本节的阅读需要傅里叶级数及傅里叶变换的相关数学知识。

示范代码目录下有一个ecgsignal.dat文件，这里存储了作者采集的一段人体心电信号-ECG。这个文件以4字节浮点数存储样本，单位为μV，采样总数 = 文件大小 / 4，采样频率 = 2000样本/秒。需要说明的是，这个心电信号不是标准的医用心电信号，作者在一台其它用途的医用电生理设备上，用左手拿着正电极，右手拿着负电极，简单记录了上述信号。而且，作者故意没有涂用于皮肤电极的导电膏，以便引入“工频干扰”。

版权声明

本文可以在互联网上自由转载，但必须：注明出处(作者：海洋饼干叔叔)并包含指向本页面的链接。

本文不可以以纸质出版为目的进行改编、摘抄。

本文节选自作者的《Python编程基础及应用》视频教程。

1. 分析信号的频谱 - spectrum

心电信号可以视为定义在时间t上的函数，把这个函数进行傅里叶级数展开，可以将其表达成不同频率/周期的正弦函数的和。而这些正弦项的系数，则表明了该种频率的正弦项在信号构成中的重要程度。所谓的频谱分析，就是把时间t上的函数，转换成频率f上的函数，即把信号从时域转换到频域。这种转换，可以通过傅里叶变换实现。在离散的计算机世界里，对应的算法工具称为快速傅里叶变换 - Fast Fourier Transform，简称FFT。

#Spectrum.py
import numpy as np
from matplotlib import pyplot as plt

iSampleRate = 2000					#采样频率,每秒2000个样本
x = np.fromfile("ecgsignal.dat",dtype=np.float32)
iSampleCount = x.shape[0]			#采样数
t = np.linspace(0,iSampleCount/iSampleRate,iSampleCount)

xFFT = np.abs(np.fft.rfft(x)/iSampleCount)  #快速傅里叶变换
xFreqs = np.linspace(0, iSampleRate/2, int(iSampleCount/2)+1)

plt.figure(figsize=(10,6))
ax0 = plt.subplot(211)             #画时域信号
ax0.set_xlabel("Time(s)")
ax0.set_ylabel("Amp(μV)")
ax0.plot(t,x)

ax1 = plt.subplot(212)             #画频域信号-频谱
ax1.set_xlabel("Freq(Hz)")
ax1.set_ylabel("Power")
ax1.plot(xFreqs, xFFT)
plt.show()

执行结果：

先点下方工具栏中的放大镜，然后按住鼠标左键框选放大，可以帮助我们帮察信号及频谱的细节：

np.fromfile(“ecgsignal.dat”,dtype=np.float32)从文件读取信号数据，numpy会以二进制格式打开文件，读取数据，并以4个字节为单位，逐一转换成np.float32类型，然后返回一个一维数组。该数组包含全部采样。 np.float32类型表示一个采样浮点数由32个bit，即4个字节构成。

np.linspace(0,iSampleCount/iSampleRate,iSampleCount)则生成了与信号x相同维度的一维数组t，其中的数据对应每个样本的采样时间。其中，iSampleCount为采样总数，iSampleRate为采样频率。如果采样总数=6000，则信号的总时间长度为6000/2000(采样频率)=3秒。后续代码中的ax0.plot(t,x)则以时间t为横轴，信号振幅x为纵轴，在ax0子图上画出时域信号。

np.fft模块支持快速傅里叶变换，其中的rfft()函数对实数信号进行FFT运算。根据计算公式，还需要将返回结果除以iSampleCount以便正确显示波形能量。

rfft()函数的返回值是iSampleCount/2+1个复数，分别表示从0(Hz)到iSampleRate/2(Hz)的频率能量。np.abs()对rfft()返回数组中的复数进行求模-abs。受限于香农采样定理，采样频率为2000Hz的信号，有效的信号频率最高为1000Hz。

xFreqs = np.linspace(0, iSampleRate/2, int(iSampleCount/2)+1)则负责生成与xFFT中的能量相对应的频率。此时,(xFreqs,xFFT)可视为定义在频率域xFreqs上的“信号”，即原时域信号(t,x)的频谱。后续代码中的ax1.plot(xFreqs, xFFT)则以频率为横轴，能量为纵轴，在ax1子图上画出频谱。

顺便再次指出，plt.subplot(211)将当前图按2行1列划分，在1号位置（即上方位置）创建一张子图ax0，同理，plt.subplot(212)则在2号位置（即下方位置）创建一张子图ax1。

在下方频谱图中，我们看到在50Hz处存在一个巨大的峰，这些能量对应着那些叠加在原始信号上的有规律的周期小波，这就是所谓的工频干扰。读者可以数一数下方上图中从3.0秒至3.2秒的周期小波的周期(峰到峰）个数，应该是10个，0.2秒10个周期对应1秒种50个周期。我国的市电是220V/50Hz交流电，交变的电流流过市电导线，向空间辐射能量，这些辐射能量借助于人体，电极线等与市电导线的耦合电容，进入了信号放大器，形成“干扰”。此外，我们在100Hz处也可以看到频谱能量的峰，100Hz正好是50Hz 的倍频，也是所谓工频干扰的一部分。可以想象，如果同样的干扰发生在美国，其频率应为60Hz，因为当地的市电是110V/60Hz交流电。这里，50Hz的工频干扰我们将使用带阻滤波器，也叫陷波器滤除。

在频谱图中，我们还看到在0Hz附近也有较多能量，这些低频成分对应原始信号里缓慢变化的基线。这些低频成分也不具备诊断意义，需要滤除。我们选择使用一个带通滤波器，滤除3Hz以下的低频信号，同时滤除70Hz以上的高频信号。

2. 滤波器设计

滤波是个宏大的课题，这里我们只能描述一种简便的应用方法，不描述背后的数学原理。

我们需要两个滤波器，其中一个是3-70Hz的带通滤波器，它保留信号中3-70Hz的频率成分，去除低于3Hz的低频部分以及高于70Hz的高频部分。另外一个是48-52Hz的带阻滤波器，别名50Hz陷波器，它去除信号中48-52Hz的成分。

下图展现了我们设计的带通滤波器(左)及带阻滤波器(右)的频率响应。横轴为频率，纵轴为滤波器对该频率信号的增益-gain，当增益为1.0时，说明该频率信号无障碍通过滤波器，当增益小于1.0时，说明通过滤波器时，该频率信号被衰减。可以看到，当滤波器的阶-order越高，则滤波器性能越好(频响曲线陡峭)，但计算量也会增加。

相关滤波器设计代码如下:

#FilterResponse.py
def butterBandPassFilter(lowcut, highcut, samplerate, order):