基于自适应阈值小波变换的心电信号降噪方法与流程

本发明涉及心电信号降噪领域，特别涉及一种基于自适应阈值小波变换心电信号降噪方法。

背景技术：

心脏病及各种心血管疾病是世界范围内致死、致残率较高的疾病，其突发性与不可预测性使得诊治受到很大限制，随着我国人口老龄化以及生活水平的不断提高，心脏病和各种心血管疾病的发病率呈持续上升趋势，我国人口众多，医生和医疗设备相对匮乏，开展基于大数据的心电智能监测可有效降低致死致残率，其中对采集到的心电信号进行降噪预处理是关键步骤；另一方面，心电信号是一种微弱的生理信号，频率范围为0.5-150Hz，90％集中在0.5-40Hz，电压幅度范围一般在10μV-5mV，具有低频率、低幅值和人体高阻抗等特点，同一时间可受到工频干扰(50Hz/60Hz)、基线漂移(<1Hz)以及肌电干扰(5Hz-2kHz)等噪声影响，使得心电信号波形模糊不清，对随后的信号分析处理、计算机自动识别诊断造成很大影响。心电信号降噪已然成为心电监护、心脏疾病预测和诊断的核心部分。

现有技术中小波阈值去噪是一种非常有效的心电信号去噪方法，小波阈值选取算法以Donoho等人提出的通用阈值(VisuShrink)算法为主，该方法在每层上都使用同一个阈值，有“过扼杀”小波系数的倾向，信号重构后较之原始信号易产生震荡，重构平滑效果不佳。

技术实现要素：

针对上述问题，本发明目的在于提供一种基于自适应阈值小波变换的心电信号降噪方法，实现自适应各层的阈值选择方式，灵活地将噪声和信号分离，改善使用通用阈值的不足，在视觉、均方误差和信噪比三方面均有改善，具有较高的实用价值的。

为达到上述目的，本发明的提出的技术方案为：基于自适应阈值小波变换的心电信号降噪方法，其特征在于，包括如下步骤：

步骤1：采用Mallat算法，选择小波函数sym6和分解层数J，对含噪心电信号进行小波分解，得近似系数和细节系数；

步骤2：设定自适应于各层细节系数的阈值，并选择阈值函数；

步骤3：对各层的细节系数进行自适应阈值处理，去除工频干扰和肌电干扰，对近似系数处理去除基线漂移；

步骤4：对处理后的心电信号进行小波重构，得到信号的近似最优估计值。

优选的，步骤1中所述的分解层数J为根据实际需求再结合采样率和含噪信号长度确定分解最大层数为8层。

进一步的，步骤2中所述的设定自适应于各层细节系数的阈值，并通过自适应阈值计算方式计算得到第一层细节系数标准偏差、各层层细节系数标准偏差以及层参数，从而得到各层细节系数对应的阈值，具体计算如下：

阈值Th求取公式为：

其中为第一层细节系数标准偏差，层层细节系数标准偏差；

层参数β在1-3层上求取公式为：

层参数β在4-8层上求取公式为：

其中Lk表示层层细节系数长度，L表示获取的含噪信号长度，J表示小波分解最大层，j表示当前分解层；

其中|Wij|是小波变换的第一层细节系数；

|Yij|为小波变换的各层细节系数；

阈值函数采用软阈值方法，定义如下：

进一步的，步骤3所述的对近似系数处理去除基线漂移为将近似系数直接置为零。

进一步的，步骤3所述的对各层的细节系数进行自适应阈值处理，去除工频干扰和肌电干扰为将细节系数与该层下的自适应阈值使用软阈值函数进行对比，若系数小于阈值，则置为0，否则将系数的绝对值减去阈值，再将其添加正负号，最终得到原始信号系数的估计系数。

进一步的，步骤4所述的对阈值处理后的小波系数进行重构为利用Mallat算法对去噪以后的小波系数进行重构，得到原始信号的近似最优估计值。

采用上述技术方案，本发明所述的基于自适应阈值小波变换的心电信号降噪方法，充分利用小波变换多分辨率特性，采用自适应阈值选取方式，在每一级层上采用不同阈值，灵活地将噪声和信号分离，改善了使用通用阈值的不足，提高了信号特征的可分离性，在视觉、均方误差和信噪比三方面较传统方法效果更好，较好的保留了图像的细节信息，具有较高的实用价值。

附图说明

图1为本发明方法的总体流程图；

图2为MIT-BIH心律不齐数据库中118e_6号数据前1800个数据的原始信号波形图；

图3为MIT-BIH心律不齐数据库中118e_6号数据前1800个数据的原始信号频谱图；

图4为MIT-BIH心律不齐数据库中118e_6号数据前1800个数据去除基线干扰后的频谱图；

图5为MIT-BIH心律不齐数据库中118e_6号数据前1800个数据使用传统方法降噪后的波形图；

图6为MIT-BIH心律不齐数据库中118e_6号数据前1800个数据使用传统方法降噪后的频谱图；

图7为MIT-BIH心律不齐数据库中118e_6号数据前1800个数据使用本发明方法降噪后的波形图；

图8为MIT-BIH心律不齐数据库中118e_6号数据前1800个数据使用本发明方法降噪后的频谱图。

具体实施方式

下面结合附图和具体实施列，对本发明做进一步说明。

实施例1

本发明所述的基于自适应阈值小波变换的心电信号降噪方法，在Matlab中的操作包括：

步骤1：获取含噪的心电信号：选取MIT-BIH噪声数据库中的118e_6号数据的前1800个数据进行定性分析实验；在Matlab中绘制选取数据的原始波形图(图2)和频谱图(图3)，作为分别使用传统方法和本发明方法降噪后的对比参照图；

步骤2：选择小波函数φ(x)和分解层数J，对含噪心电信号进行小波分解：本发明选用与心电信号形态近似的sym6小波函数；118e_6号数据采样率360Hz，有用信号90％集中在0.5-40Hz，基线漂移干扰主要集中在1Hz以内的低频部分，所以进行8层小波分解；在Matlab中对步骤1中的数据使用sym6小波基进行8层分解，得到近似系数cA8和细节系数cD1、cD2、cD3、cD4、cD5、cD6、cD7、cD8；

步骤3：将步骤2中分解得到的代表基线漂移的cA8直接置零去除基线漂移，去除基线漂移在Matlab绘制得到去除基线干扰后的频谱图(图4)；

步骤4：设定自适应于各层细节系数的阈值，并选择阈值函数：这一步骤需要分别使用传统方法得到通用阈值，使用发明方法得到步骤2中的8层每层对应的自适应阈值；

传统方法阈值求取公式为其中n为信号样本点个数，是估计的噪声标准偏差，求取公式为其中|Yij|为小波分解的第一层细节系数；在Matlab中根据公式计算得到n＝1800，Th＝0.0302；

使用本发明方法的阈值求取方式，即本发明设定自适应于各层细节系数的阈值，阈值Th求取公式为：

其中为第一层细节系数标准偏差，各层层细节系数标准偏差；

层参数β在1-3层上求取公式为：

层参数β在4-8层上求取公式为：

其中Lk表示层层细节系数长度，L表示获取的含噪信号长度，J表示小波分解最大层，j表示当前分解层；

其中|Wij|是小波变换的第一层细节系数；

其中|Yij|为小波变换的各层层细节系数；

阈值函数采用软阈值方法，定义如下：

在Matlab中代入公式计算可得各层的阈值从cD1-cD8分别为1.5414、0.4032、0.1303、0.0818、0.0208、0.0126、0.0191、0.0174。广泛应用的阈值施加方式有硬阈值函数和软阈值函数，其中硬阈值函数在阈值点容易产生震荡，重构效果不佳，我们选择软阈值函数，

对各层的小波系数进行自适应阈值处理，去除工频干扰、基线漂移和肌电干扰：对步骤1中的数据分别使用步骤4中提到的两种阈值均结合软阈值函数进行降噪处理，得降噪后的波形图(图5)和频谱图(图6)，以对比体现出本发明的优势；

步骤5：对处理后的心电信号进行小波重构，得到信号的近似最优估计值；对步骤4中得到的两组处理后的小波系数分别进行一次小波重构操作，得到原始含噪信号的最优估计信号，得采用本发明方法降噪后的波形图(图7)和降噪后的频谱图(图8)。

由附图5-8可知，MIT-BIH噪声数据库中118e_6号数据噪声主要是基线漂移和少量肌电噪声，工频干扰不明显，改进方法处理结果(图7、图8)和传统方法处理结果(图5、图6)相比，基线漂移和肌电噪声很好的被去除，波形连续、清晰且无震荡，去除噪声效果良好。

实施例2

基于自适应阈值小波变换的心电信号降噪方法，在Matlab中的操作包括：

步骤1：获取含噪的心电信号，选取MIT-BIH心律不齐数据库中103号数据的20秒到25秒共1800个数据作为原始“纯净”信号，在其上叠加常见的三种干扰，首先叠加频率为60Hz，幅度为0.02mV的正弦信号模拟工频干扰，然后使用awgn()函数叠加指定信噪比的随机高斯白噪声模拟肌电干扰，最后叠加频率为0.3Hz,幅度为0.2mV的正弦信号模拟基线漂移，从而得到待降噪的实验数据；在叠加随机白噪声时，加入不同信噪比的白噪声，处理结果不同，本实施例以叠加噪声后信噪比为2.2819为例，分别使用传统方法和本文提到的改进方法进行阈值降噪处理，通过计算SNR(信噪比)和MSE(均方误差)来体现发明方法的改进的阈值选取方式的优势；

步骤2：选择小波函数φ(x)和分解层数J，对含噪心电信号进行小波分解；选用与心电信号形态近似的sym6小波函数；103号数据采样率360Hz，有用信号90％集中在0.5-40Hz，基线漂移干扰主要集中在1Hz以内的低频部分，所以进行8层小波分解；在Matlab中对实施例2中的数据使用sym6小波基对其进行8层分解，得到近似系数cA8和细节系数cD1、cD2、cD3、cD4、cD5、cD6、cD7、cD8；

步骤3：将步骤2中分解得到的代表基线漂移的cA8直接置零去除基线漂移；

步骤4：设定自适应于各层细节系数的阈值，并选择阈值函数，这一步骤需要分别使用传统方法得到通用阈值、使用本发明方法得到实施例2中的8层对应的自适应阈值；传统阈值求取方式与实施例1相同，在Matlab中根据公式计算得到n＝1800，Th＝1.3643；使用本发明方法的阈值求取方式，采用实施例1的计算方式，在Matlab中代入公式计算可得各层的阈值从cD1-cD8分别为1.5425、1.3790、1.0830、0.7866、0.7026、0.5217、0.3994、0.7581；阈值函数选择重构效果更佳光滑的软阈值函数，与实施例1中的相同；对各层的小波系数进行自适应阈值处理，去除工频干扰、基线漂移和肌电干扰；对步骤3的数据分别使用两种阈值均结合软阈值函数进行降噪处理，以对比体现出本发明的优势；将剩余的细节系数分别使用实施例1中计算得到通用阈值和自适应阈值进行处理；

步骤5：对处理后的心电信号进行小波重构，得到信号的近似最优估计值。对步骤4中得到的两组处理后的小波系数分别进行一次小波重构操作，得到原始含噪信号的估计信号；不同初始信噪比经过降噪处理后的信噪比对比见下表1，不同初始信噪比经过降噪处理后的均方误差对比见下表2。

表1

表2

评估降噪性能：分别计算传统方法和本发明方法阈值选取方式降噪后的SNR(信噪比)和MSE(均方误差)，SNR和MSE的计算公式如下：

其中N表示信号样本点个数为1800，s(t)为步骤1中获取的“纯净”信号，为降噪后的信号估计值。由表1可知，采用本发明方法处理后的信噪比值明显大于通用阈值方法，由表2可知，采用本发明方法降噪后的信号与“纯净”信号误差在相同的信噪比下小于采用第二列现有技术处理的结果。

尽管结合优选实施方案具体展示和介绍了本发明，但所属领域的技术人员应该明白，在不脱离所附权利要求书所限定的本发明的精神和范围内，在形式上和细节上对本发明做出各种变化，均为本发明的保护范围。