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 pltiSampleRate = 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):"生成巴特沃斯带通滤波器"semiSampleRate = samplerate*0.5low = lowcut / semiSampleRatehigh = highcut / semiSampleRateb,a = signal.butter(order,[low,high],btype='bandpass')print("bandpass:","b.shape:",b.shape,"a.shape:",a.shape,"order=",order)print("b=",b)print("a=",a)return b,adef butterBandStopFilter(lowcut, highcut, samplerate, order):"生成巴特沃斯带阻滤波器"semiSampleRate = samplerate*0.5low = lowcut / semiSampleRatehigh = highcut / semiSampleRateb,a = signal.butter(order,[low,high],btype='bandstop')print("bandstop:","b.shape:",b.shape,"a.shape:",a.shape,"order=",order)print("b=",b)print("a=",a)return b,aiSampleRate = 2000 #采样频率plt.figure(figsize=(12,5))
ax0 = plt.subplot(121)
for k in [2, 3, 4]:b, a = butterBandPassFilter(3,70,samplerate=iSampleRate,order=k)w, h = signal.freqz(b, a, worN=2000)ax0.plot((iSampleRate*0.5/np.pi)*w,np.abs(h),label="order = %d" % k)ax1 = plt.subplot(122)
for k in [2, 3, 4]:b, a = butterBandStopFilter(48, 52, samplerate=iSampleRate, order=k)w, h = signal.freqz(b, a, worN=2000)ax1.plot((iSampleRate*0.5/np.pi)*w,np.abs(h),label="order = %d" % k)

控制台输出：

bandpass: b.shape: (5,) a.shape: (5,) order= 2
b= [ 0.00961483  0.         -0.01922967  0.          0.00961483]
a= [ 1.         -3.70024346  5.14358108 -3.18588885  0.74255496]
bandpass: b.shape: (7,) a.shape: (7,) order= 3
...
bandpass: b.shape: (9,) a.shape: (9,) order= 4
...

下述讨论中，f₀等于采样频率的一半，即代码中的iSampleRate*0.5，或者samplerate*0.5。

从控制台输出可以看到，signal.butter()函数生成滤波器的结果为两个一维数组b和a。数组里存储一些实数，滤波器的阶-order越大，b,a包含的实数个数越多。b,a是IIR滤波器的系数。btype参数指明了滤波器的类型，bandpass意为带通，bandstop意为带阻。[low,high]则指明了滤波器的截止频率。上述代码说明，这里的low和high应等于截止频率/f₀。当采样频率为2000时，f₀=2000*0.5=1000，则3Hz的截止频率对应low = 3/1000 = 0.003，70Hz的截止频率对应high = 70 / 1000 = 0.07。

freqz()函数用于计算滤波器的频率响应，返回w,h两个数组。其中，w是圆频率数组，通过w*f₀/π可以计算出与其对应的频率，h是w对应频率点的响应，它是一组复数，其幅值表示滤波器的增益特性，相角表示滤波器的相位特性。参数worN表明了要计算的频率的项数，该值越大，计算越精细。

上述代码中，我们生成了阶为2,3,4的3-70Hz带通滤波器共3个，阶为2,3,4的48-52Hz带阻滤波器共3个，然后分别生成并显示了其频率响应曲线，以供读者观察。上述代码中，np.abs()函数用于求复数数组的模。

上图中，我们还画了一条sqrt(0.5)的虚线，这里对应着所谓的-3dB增益。此处，对应频率的信号的功率缩减为最高功率的一半。

3. 滤波

我们应用上述滤波器对信号进行了滤波，结果如下图（局部放大）：

可以看到，在3-70Hz带通滤波器的作用下，0Hz附近的极低频成分消失了，70Hz以后的高频成分也得到有效抑制。同时，我们注意到100Hz的干扰成分仍有残留。

在48-52Hz带阻滤波器的作用下，50Hz附近工频干扰几乎完全消失。上图中，我们看到了基线不飘移，50Hz工频周期波完全去除后的“干净”的ECG信号。这个信号来自于心脏，常在医用心电监护仪上看到。每个“尖波”对应着一次心跳，读者可以计算一下作者记录这段信号时的心率。