SciPy实践 - ECG信号的谱分析、数字滤波及相位校正

  阅读本部分内容需要傅里叶变换以及数字滤波器方面的相关知识。

  本书随书代码CH10目录下有一个ecgsignal.dat文件，其中存储了一段人体心电信号（electrocardiogram，简称ECG）。该文件以4字节浮点数格式存储样本，单位为μV，采样总数等于文件尺寸/ 4，采样频率2000样本/秒。需要说明的是，这个心电信号不是标准的医用心电信号，作者在一台其他用途的医用电生理设备上，用左手拿着正电极，右手拿着负电极，简单记录了上述信号。作者故意没有涂抹降低皮肤电极阻抗的导电膏，以便引入“工频干扰”，即由交流市电网络辐射进入信号采集系统的频率为50Hz的周期性干扰波。

版权声明

本文节选自《Python编程基础及应用》第二版， 高等教育出版社，作者陈波，刘慧君。

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

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

1. 信号的频谱分析

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

  下述程序从ecgsignal.dat文件读出了原始信号，然后交给scipy.signal子模块进行快速傅里叶变换，得到信号的频谱并绘图。

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, matplotlib.pyplot as plt

sampleRate = 2000					#采样频率,每秒2000个样本
x = np.fromfile("ecgsignal.dat",dtype=np.float32) 
x -= np.mean(x)                     #去直流分量
sampleCount = x.shape[0]			#采样数
t = np.linspace(0,sampleCount/sampleRate,sampleCount)  
xFFT = np.abs(np.fft.rfft(x)/sampleCount)  #快速傅里叶变换
xFreqs = np.linspace(0, sampleRate/2, sampleCount//2+1)

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

ax1 = plt.subplot(2,1,2)             #下方子图画频谱
ax1.set_xlabel("Freq(Hz)")
ax1.set_ylabel("Power")
ax1.plot(xFreqs, xFFT)
plt.subplots_adjust(0.1,0.1,0.95,0.95,hspace=0.2)
plt.show()

  上述程序的绘图结果见图10-E1。

  用鼠标单击下方工具栏中的放大镜，然后按住左键框选放大，可以帮助我们帮察信号及频谱的细节，请见图10-E2。

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

🚩第7行：将x数组中的每个元素都减去全部元素的均值np.mean(x)，从而去除信号中的直流分量。

🚩第8行：np.linspace(0,sampleCount/sampleRate,sampleCount)则生成了与信号x相同维度的一维数组t，其中的数据对应每个样本的采样时间。其中，sampleCount为采样总数，sampleRate为采样频率。本例中，采样总数为10942，采样频率2000，信号的总时间长度为10942/2000等于5.471秒。代码第16行ax0.plot(t,x)以时间t为横轴坐标，信号电压x为纵轴坐标，在ax0子图上画出时域信号。

🚩第9行：np.fft模块支持快速傅里叶变换，其中的rfft()函数对实数信号x进行FFT运算。根据计算公式，还需要将返回结果除以iSampleCount以便正确显示波形能量。rfft()函数的返回iSampleCount//2+1个复数，分别表示从0(Hz)到iSampleRate/2(Hz)的频率能量。np.abs()对rfft()返回数组中的复数求模(abs)。受限于香农采样定理，采样频率为2000Hz的信号，有效的信号频率最高为1000Hz。

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

🚩第13、18行：plt.subplot(2,1,1)将当前图按2行1列划分，在编号1的单元格（上方）上创建子图ax0；同理，plt.subplot(2,1,2)则在2号单元格（下方）位置创建子图ax1。

  在下方频谱图中，我们看到在50Hz处存在一个巨大的峰，这些能量对应着那些叠加在原始信号上的有规律的周期小波，这就是所谓的工频干扰。读者可以数一数上方时域信号中从3.0至4.0秒的周期小波的周期（峰到峰）数，应该是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的成分。

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

相关滤波器的设计代码如下：

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
38

#filterResponse.py
import numpy as np,scipy.signal as signal
import matplotlib.pyplot as plt

def butterFilter(lowCut, highCut, sampleRate, order, filterType):
    semiSampleRate = sampleRate*0.5
    low,high = lowCut/semiSampleRate, highCut/semiSampleRate
    b,a = signal.butter(order,[low,high],btype=filterType)
    print(f"{filterType} filter, order = {order}")
    print("b.shape:",b.shape,"a.shape:",a.shape)
    print("b=",b, "\na=",a)
    return b,a

sampleRate = 2000
plt.figure(figsize=(12,5))
ax0 = plt.subplot(121)
for k in [2, 3, 4]:
    b, a = butterFilter(3,70,sampleRate, k,"bandpass")
    w, h = signal.freqz(b, a, worN=2000)
    ax0.plot(w*sampleRate*0.5/np.pi,np.abs(h),label=f"order {k}")

ax1 = plt.subplot(122)
for k in [2, 3, 4]:
    b, a = butterFilter(48, 52,sampleRate,k,"bandstop")
    w, h = signal.freqz(b, a, worN=2000)
    ax1.plot(w*sampleRate*0.5/np.pi,np.abs(h),label=f"order {k}")
ax1.set_xlim(35,65)

for a in [ax0,ax1]:
    a.plot([0, 0.5 * sampleRate], [np.sqrt(0.5), np.sqrt(0.5)],
             '--', label='sqrt(0.5)')
    a.set_xlabel('Frequency (Hz)')
    a.set_ylabel('Gain')
    a.grid(True)
    a.legend(loc='best')

plt.subplots_adjust(0.05,0.05,0.95,0.95,wspace=0.2)
plt.show()

上述程序的输出结果为：

1
2
3
4
5
6
7
8
9

bandpass filter, order =  2
b.shape: (5,) a.shape: (5,)
b= [ 0.00961483  0.         -0.01922967  0.          0.00961483] 
a= [ 1.         -3.70024346  5.14358108 -3.18588885  0.74255496]
bandpass filter, order =  3
......   #此处有多行节略
bandpass filter, order =  4
b.shape: (9,) a.shape: (9,)
......   #此处有多行节略

  下述讨论中，f0等于采样频率的一半，即代码中的sampleRate*0.5，本例中，f0=2000 *0.5=1000。

🚩第5 ~ 12行：生成巴特沃斯滤波器。从控制台输出可以看到，signal.butter()函数生成滤波器的结果为两个一维数组b和a。数组里存储一些实数，滤波器的阶（order）越大，b、a包含的实数个数越多。b、a是IIR（无限冲击响应）滤波器的系数。btype/filterType参数指明了滤波器的类型，bandpass意为带通，bandstop意为带阻。[low,high]则指明了滤波器的截止频率。上述代码说明，这里的low和high应等于截止频率/f0。本例中，3Hz的截止频率对应low = 3/1000 = 0.003，70Hz的截止频率对应high = 70 / 1000 = 0.07。

🚩第17 ~ 20行：生成2、3、4阶的巴特沃斯带通滤波器，计算滤波器的频率响应，并在左子图ax0中绘制频响曲线。freqz()函数用于计算滤波器的频率响应，返回w、h两个数组。其中，w是圆频率数组，通过w*f0/π可以计算出与其对应的频率，h是w对应频率点的响应，它是一组复数，其幅值表示滤波器的增益特性，相角表示滤波器的相位特性。参数worN表明了要计算的频率的项数，该值越大，计算越精细。

🚩第23 ~ 26行：生成2、3、4阶的巴特沃斯带阻滤波器，计算滤波器的频率响应，并在右子图ax1中绘制频响曲线。

🚩第27行：设置右子图ax1的x轴坐标显示范围，以便更好地观察不同阶带阻滤波器的频响差异。

🚩第29 ~ 31行：在左右子图中画sqrt(0.5)虚线，其对应所谓的-3dB增益。在该增益处，对应频率的信号的功率缩减为最高功率的一半。

3. 滤波

  应用上述带通和带阻滤波器对上述ECG信号进行滤波后，结果如图10-E4所示（局部放大）。

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

  在48 ~ 52Hz带阻滤波器的作用下，50Hz附近工频干扰几乎完全消失。上图中，我们看到了基线稳定不飘移，50Hz工频干扰波完全去除后的“干净”的ECG信号。这个信号来自于心脏，常在医用心电监护仪上看到。每个“尖波”对应着一次心跳，用60秒除以两个尖波之间的时间间隔，即得“病人”的心率。

  与滤波相关的代码片段如下：

1
2
3
4
5
6
7
8
9
10
11
12

#iirfilter.py
...    #butterFilter()函数定义略，同filterResponse.py

sampleRate = 2000                     #采样频率,每秒2000个样本
x = np.fromfile("ecgsignal.dat",dtype=np.float32)
x -= np.mean(x)                       #去直流分量
b,a = butterFilter(3,70,sampleRate,4,"bandpass")   
x = signal.lfilter(b,a,x)           #进行带通滤波
b,a = butterFilter(48,52,sampleRate,2,"bandstop")
x = signal.lfilter(b,a,x)           #进行带阻滤波
	
...   #谱分析及画图部分代码略

🚩第4 ~ 6行：读入原始信号样本序列，并去除其中的直流分量。

🚩第7 ~ 8行：首先调用butterFilter()函数生成4阶3 ~ 70Hz带通滤波器的b、a参数数组，然后应用signal.lfilter()函数进行滤波。

🚩第9 ~ 10行：对前述带通滤波后的信号进行2阶48 ~ 52Hz的带阻滤波。

  滤波是由signal.lfilter()函数完成的，该函数将b、a数组应用于信号数组x上，计算并返回一个新数组。巴特沃斯滤波器属于IIR（无限冲击响应）滤波器，对于输出信号数组y的任意样本yi，lfilter()函数按下述公式计算其值：

  其中N为b、a参数数组的元素个数。

  回顾filterReponse.py程序的输出，可以发现a0恒为1.0。

  使用signal.firwin()函数可以生成FIR（有限冲击响应）滤波器。下述代码应用一个3 ~ 30Hz的FIR带通滤波器对上述ECG信号进行了数字滤波，绘图结果如图10-E5所示（局部放大）。

1
2
3
4
5
6
7
8
9
10
11

#firfilter.py
import numpy as np, scipy.signal as signal 
from matplotlib import pyplot as plt

sampleRate = 2000     #采样频率,每秒2000个样本
x = np.fromfile("ecgsignal.dat",dtype=np.float32)
x -= np.mean(x)                       #去直流分量
b = signal.firwin(1001,[3.0,30.0],fs=sampleRate,pass_zero=False)  
x = signal.lfilter(b,1,x)           #进行带通滤波

...   #频谱分析及绘图代码略

🚩第8行：使用signal.firwin()生成3 ~ 30Hz的FIR带通滤波器，滤波器参数数组长度为1001，关键字参数fs代表原始信号的采样频率。与IIR滤波器不同，firwin()所生成的FIR滤波器仅返回了b参数数组，本例中，其元素个数为1001。

  对于FIR滤波器，参数数组长度 = 滤波器阶数 + 1。本例中，长度为1001的参数数组对应1000阶的滤波器。

🚩第9行：使用FIR带通滤波器的b参数数组对ECG信号进行数字滤波。

  当使用FIR滤波器进行数字滤波时，对于输出信号数组y的任意样本yi，signal.lfilter()函数按下述公式计算其值：

  其中N为b参数数组的元素个数。

  上述公式中，a参数数组被忽略，因此，本行代码中使用整数1代替a数组。

  如图10-E5所示，在应用3 ~ 30Hz的FIR带通滤波器后，程序得到了稳定不漂移，50Hz工频干扰已被移除的“漂亮”的ECG信号（上）。下方频谱图证实，原始信号中30Hz以上的高频成分，几乎被完全移除。

4. 滤波中的相位变化

  使用signal.lfilter()函数进行数字滤波会产生相位移动/时延。对于FIR滤波器，其时延为定值，可根据滤波器阶数计算。对于IIR滤波器，时延与原始信号的频率有关，即原始信号中不同的频率分量，其时延亦不同，这会使得滤波后的信号发生“扭曲”，形状发生变化。

  如图10-E6展示了应用前述3 ~ 30Hz带通FIR滤波器所带来的信号时延。图中最上方为原始信号，中间为FIR滤波器的输出信号。容易看出，相对于原始信号，FIR滤波器的输出信号在去除50Hz工频干扰的同时，还使得ECG各个峰发生了时延，每个尖峰均向后推迟了约250ms。这个250ms对应FIR滤波器阶数1000的一半，即500样本，由于采样频率是2000Hz，500个样本的时间长度即为250ms。

  无论是FIR还是IIR滤波器，均可改用signal.filtfilt()函数以消除滤波过程中的相移/时延。下述程序分别使用signal.lfilter()和signal.filtfilt()函数对信号进行了滤波处理，并绘图对时延差异进行比较，绘图结果即图10-E6。

1
2
3
4
5
6
7
8
9
10
11

#firphase.py
...    #代码有节略
sampleRate = 2000      #采样频率,每秒2000个样本
x = np.fromfile("ecgsignal.dat",dtype=np.float32)
x -= np.mean(x)                      #去直流分量

b = signal.firwin(1001,[3.0,30.0],fs=sampleRate,pass_zero=False) 
y1 = signal.lfilter(b,1,x)           #滤波，有相位延迟
y2 = signal.filtfilt(b,1,x)          #滤波，无相位延迟   

...    #绘制原始信号x、有时延滤波信号y1、无时延滤波信号y2，代码略

🚩第7行：生成3 ~ 30Hz的FIR带通滤波器，得滤波器系数数组b。

🚩第8行：使用signal.lfilter()进行滤波，如前所述，滤波结果信号y1有时延。

🚩第9行：使用signal.filtfilt()使用同一滤波器进行滤波，如前所述，滤波结果信号y2无时延。

  图10-E6中最下方的曲线即为signal.filtfilt()函数所对应的滤波结果信号，其中所有的信号峰均与最上方的原始信号完全对齐，滤波的时延被完美消除。

  signal.filtfilt()函数是通过一次顺向滤波加一次反向滤波来消除相移/时延的。简单地说，signal.filtfilt()函数首先按照与signal.lfilter()相同的计算公式，从前往后进行一轮滤波，产生了正向的时延；然后，再用相同的计算公式从后往前进行一轮滤波，产生了反向的时延；两个时延相互抵消，从而得到零相移/零时延的滤波结果。

  图10-E7则展示了应用前述3 ~ 70Hz带通IIR滤波器和48 ~ 52Hz带阻IIR滤波器进行有相移及无相移滤波后的结果差异。图上方为原始信号，中间为signal.lfilter()的有相移滤波信号，下方则为signal.filtfilt()的无相移滤波信号。看起来，上中下三个信号各个ECG峰的时间差异不大，这是因为IIR滤波器的阶相对较小。但如果仔细观察，在应用相同的IIR带通和带阻滤波器后，中间、下方两个结果信号的波形形状差异很大，这是因为IIR滤波器对于不同的原始信号分量有不同的相移/时延，导致结果信号“扭曲”。显然，signal.filtfilt()函数输出的下方无相移滤波信号更忠实于原始信号。

  图10-E7由程序iirphase.py绘制，请参见随书源代码CH10子目录。

本例数据及程序包下载链接：

http://codelearn.club/download/ecgfilter.zip