小波变换在心电图ECG中的应用

1.算法介绍

首先来看看什么是小波变换？

要回答这个问题，首先得从傅里叶变换FFT开始说起，了解过FFT的人都知道，对一段信号进行傅里叶变换，可以知道这段信号有哪些频率成分存在，比如：

fs = 250; N = 10000;   %采样频率和数据点数
n = 0:N-1;t = n/fs;   %时间序列
x = sin(2*pi*50*t)+sin(2*pi*60*t); 
y = fft(x,N);    %对信号进行快速Fourier变换
mag = abs(y);    %求得Fourier变换后的振幅
n1 = 0:N-1;
f = n*fs/N;    %频率序列
plot(f,mag);    %绘出随频率变化的振幅

运行结果：

图1 FFT变换

模拟一个50Hz和60Hz叠加的信号源，对信号源做FFT变换，可以看到在50Hz和60Hz都分别出现了高峰，代表该信号主要的频率分量为50Hz和60Hz。这也是著名的频域分析。

这里就引出一个问题，对信号做FFT可以知道整个信号包含哪些频率，和频谱能量分布情况，对于平稳信号来说是没问题，自然界大部分都是非平稳信号，这个时候我想知道某个频率出现的时间点，FFT就办不到了，FFT只是频域分析，没有时间的概念，当然有人就会说，对某段信号再进行FFT不就好了，也就是短时傅里叶变换，当然这种也是可以的，但是窗口设计大小很难去设定，有兴趣的可以去了解一下，这里不展开讲。

另外一种更好的方法就是本章的重点，也就是小波变换，先来看小波变换的公式：

小波变换有两个变量，一个是伸缩 a、一个是平移τ，伸缩反映的是频率情况，平移反映的是时间，假定用小波不同尺度（伸缩）下都与源信号进行运算，那最终能得到的是任意位置（时刻）的频率情况，即时频图，这样我们就可以知道每个位置（时刻）的频率成分了，这就是非常著名的小波变换了，小波变换就像信号处理里面的一个显微镜，能很好的分析出信号在时间、频率上的信息。

2.算法实现

实现小波变换，主要分为两块，一是小波基的选取，二是小波变换算法。

常见的小波基，如Haar、Daubechies、Coiflets、Symlets、Meyer等等，小波基的选取主要考虑支撑长度、对称性、正则性等，这里不展开讲，需要根据实际情况进行选取，有兴趣的朋友可以了解下。

小波变换算法，目前最经典的莫过于Mallat快速算法，Mallat快速算法主要分为两大块，一是对信号进行分解，二是对分解后的信号进行重构，这里借助一个例程分别讲分解和重构的过程。

假设信号为{{420.197,423.526,423.520,423.351,424.520,428.001,430.793,428.925}}，利用db2小波基对信号进行分解

要得到db2小波基，用Matlab获得：

>> [Lo_D,Hi_D,Lo_R,Hi_R] = wfilters('db2')

Lo_D =
   -0.1294    0.2241    0.8365    0.4830

Hi_D =
   -0.4830    0.8365   -0.2241   -0.1294

Lo_R =
    0.4830    0.8365    0.2241   -0.1294

Hi_R =
   -0.1294   -0.2241    0.8365   -0.4830

信号分解的流程如下：

1.对信号进行延拓，这里采用对称延拓方式延拓源信号，得到信号：

这时信号长度变为 n +(filterlen-1)*2 = 8 + (4 - 1)*2 = 14；

2. 得到近似分量

2.1将第1步得到的结果与低通滤波器卷积，得到结果y1：

这时信号长度变为 n+(filterlen-1)*2 +(filterlen-1) = 14 + （4-1） = 17；

2.2  去掉延拓部分，即前后数据各去掉 (filterlen-1)个数据，这时数据长度变为 17 - (filterlen-1)*2 = 17 - （4 - 1）*2 = 11；

2.3 对数据进行抽二采样，最终得到信号的近似分量app：

3 得到细节分量

3.1将第1步得到的结果与高通滤波器卷积，得到结果y2：

这时信号长度变为 n+(filterlen-1)*2 +(filterlen-1) = 14 + （4-1） = 17；

3.2  去掉延拓部分，即前后数据各去掉 (filterlen-1)个数据，这时数据长度变为 17 - (filterlen-1)*2 = 17 - （4 - 1）*2 = 11；

3.3 对数据进行抽二采样，最终得到信号的细节分量coef：

为了验证准确性，我们直接采用Matlab验证：

信号的重构流程如下：这里我们直接使用分解流程得到的近似分量app和细节分量coef作为输入

1.选择近似分量app,隔点插零，这时数据长度变为 2n + 1 = 11;

2.信号延拓，数据长度变为：   2n + 1 + （filterlen -1）*2 = 17；

3.与重构的低通滤波器卷积，得到信号y1:

y1 长度变为 2n + 1 + （filterlen -1）*2 + （filterlen -1） = 20；

4.选择细节分量coef,隔点插零，这时数据长度变为 2n + 1 = 11;

5.信号延拓，数据长度变为：   2n + 1 + （filterlen -1）*2 = 17；

6.与重构的高通滤波器卷积，得到信号y2:

y2 长度变为 2n + 1 + （filterlen -1）*2 + （filterlen -1） = 20；

7.将y1与y2相加，得到信号长度依然是20，去除前面及后面（filterlen -1）*2个无关数据，得到最终数据，数据长度为 2n + 1 - （filterlen -1） = 8

可以看到，通过db2小波基对原始信号进行分解之后，再重构回来，得到的结果和源数据一样（重构后的数据长度可能会不一致，这个不能用来判定算法是否有效，我这里是用Matlab验证过）。

3.算法应用

前面我们已经对小波算法有大概的了解，这里已医疗领域，处理心电图数据的基线漂移为例来展开。

由于各种原因，测人体心电图时，原始数据总会带来一些基线漂移，频率一般在1Hz以内，为了去除这种基线漂移，我们常常会使用到小波变换，对信号进行分析。

假如心电图ECG波的采样率为250Hz,使用小波变化对信号进行分解时，不同的分解层，含有不同频率的信号分量，如下：

图2 小波分解

试验表明，基线漂移集中在CA7上（这个可能需要根据自己的系统实测得到），这时我们可以通过将CA7重构，得到信号L，信号L即为心电图的基线漂移部分，这时我们只需要将原始信号减去这部分重构得到的信号，即可得到不漂移的心电图信号，

这里我直接用真实的基线偏移的心电图数据，看看滤波效果，如下：

图3 提取基线漂移

这个用Matlab工具，直接对一段严重基线漂移的心电信号分解，在第7次分解的近似分量已经很接近了基线漂移了，把第7次分解的近似分量进行重构，即可提取出干扰信号。

另外附上在实际产品中的滤波效果，人体实测效果：

图4 人体实测小波转换滤波效果