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

 

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

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

版权声明

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

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

1. 分析信号的频谱 - spectrum

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#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()

执行结果:

1545882659338

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

1545882642607

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越高,则滤波器性能越好(频响曲线陡峭),但计算量也会增加。

1545886223319

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#FilterResponse.py
def butterBandPassFilter(lowcut, highcut, samplerate, order):
"生成巴特沃斯带通滤波器"
semiSampleRate = samplerate*0.5
low = lowcut / semiSampleRate
high = highcut / semiSampleRate
b,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,a

def butterBandStopFilter(lowcut, highcut, samplerate, order):
"生成巴特沃斯带阻滤波器"
semiSampleRate = samplerate*0.5
low = lowcut / semiSampleRate
high = highcut / semiSampleRate
b,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,a

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

控制台输出:

1
2
3
4
5
6
7
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
...

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

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

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

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

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

3. 滤波

我们应用上述滤波器对信号进行了滤波,结果如下图(局部放大):

1545886138454

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

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

相关滤波的代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
#Filter.py
iSampleRate = 2000 #采样频率,每秒2000个样本
x = np.fromfile("ecgsignal.dat",dtype=np.float32)

#进行带通滤波
b,a = butterBandPassFilter(3,70,iSampleRate,order=4)
x = signal.lfilter(b,a,x)

#进行带阻滤波
b,a = butterBandStopFilter(48,52,iSampleRate,order=2)
x = signal.lfilter(b,a,x)

... #谱分析及画图部分代码略

滤波的代码十分简单。signal.lfilter()将滤波器返回的b,a数组应用于信号x上,滤波后返回一个新数组。读者一定有点好奇,lfilter()内部到底做了些什么。对于FIR-有限脉冲响应滤波器,a组系数不适用,对于输出信号数组y的任意点y[i],其值为:
$$
y_i=b_0x_i+b_1x_{i-1}+…+b_px_{i-p}
$$
​ 这就是个简单的乘积和。对于IIR-无限脉冲响应滤波器,对于输出信号数组y的任意点y[i],其值为:
$$
y_i=\frac{b_0x_i+b_1x_{i-1}+…+b_px_{i-p}-a_1y_{i-1}-a_2y_{i-2}-…-a_qy_{i-q}}{a_0}
$$
​ 上述两个公式中,p+1为b数组的元素个数;q+1为a数组的元素个数。

巴特沃斯滤波器只是IIR滤波器的一种。 使用firwin()、remez()等函数可以设计FIR滤波器。iirdesign()函数则可以帮助设计IIR滤波器。要想彻底弄明白滤波器背后的数学原理,不是一件容易的事情。


 评论