一种基于小波能量的心电信号去噪算法的制作方法

本发明涉及信号噪声处理技术领域，涉及一种基于小波能量的心电信号去噪算法。

背景技术：

心电图(ecg)信号记录了心脏的电活动，可以反映心脏状态。因此，心电图信号在医学上被广泛运用于心脏病诊断。然而，心电信号在采集过程中容易受到各种噪声的影响而产生退化，降低心电图的应用价值。导致心电信号退化的原因有很多，如来自生物体内部的肌电干扰、呼吸干扰；来自生物体外部的工频干扰、不良接地带来的其他干扰等。所以在对心电图进行分析诊断之前，对心电信号进行预处理去噪是个至关重要的过程。

由于心电信号噪声源的种类不同，对其进行去噪处理成为一个难题，在目前提出的众多方法中，小波阈值去噪法是一种行之有效的方法。比如，donoho和johnstone提出了一种适用于非平稳信号的基于离散小波变换的小波阈值降噪方法(d.l.donohoandi.m.johnstone,“idealspatialadaptationviawaveletshrinkage”,biometrika,1994,vol.81,pp.425-455.anddonoho,“de-noisingbysoftthresholding”,ieeetrans.trans.inform.theory,vol.41,pp.613-627,1995.)。harishchandrat.patil提出了一种基于小波分解的心电信号去噪阈值估计新方法，其中阈值计算采用每级的最大和最小小波系数(harishchandrat.patil,r.s.holambe,newapproachofthresholdestimationfordenoisingecgsignalusingwavelettransform,indiaconference(indicon),2013annualieeeconference,mumbai,india,13-15dec.2013.)。awal提出了一种自适应相关小波阈值技术(adwt)(awalma,mostafass,ahmadm,etal.anadaptiveleveldependentwaveletthresholdingforecgdenoising[j].biocyberneticsandbiomedicalengineering,2014,34(4):238-249.)。sayadi提出了一种多自适应仿生小波变换方法(mabwt)(sayadio,shamsollahimb.multiadaptivebionicwavelettransform:applicationtoecgdenoisingandbaselinewanderingreduction[j].eurasipjournalonadvancesinsignalprocessing,2007,2007(1):041274.)。

但是，上述方法需要对每一层的小波系数进行阈值处理，计算量大，处理效率低。

技术实现要素：

针对现有技术的上述不足，本发明实施例提供一种基于小波能量的心电信号去噪算法，根据小波能量确定需要去噪的小波系数，而无需对所有小波系数进行处理，有效缩减了计算量。

为实现上述目的，本发明实施例的技术方案为：

一种基于小波能量的心电信号去噪算法，包括：

s1、信号均值化处理及公差值的确定；

s2、求出最大分解层数，对均值化处理后的信号进行小波分解，实际的分解层数不小于该最大分解层数；

s3、根据小波分解后的各层细节系数能量，计算小波能量曲线除边界点以外的第一个极大值和第一个极小值的所在层数，根据该所在层数及其与所述公差值的大小关系，判断需要去噪的细节阶数变化点，若小波能量曲线的第一个极大值和第一个极小值均不存在，则预设所述细节阶数变化点；

s4、将细节阶数变化点之前的各层细节系数进行阈值去噪，根据细节阶数变化点所在层及之后层的细节系数、最后一层的近似系数以及经过阈值去噪后的各层细节系数对信号进行小波重构。

优选的，小波分解中使用的小波基是db5。

优选的，设待求的所述细节阶数变化点所在层数为chg_locs，小波能量曲线第一个极大值所在层数为max_locs(1)，第一个极小值所在层数为min_locs(1)，步骤s3为：

s3a、当所述第一个极大值和第一个极小值同时存在时：

当max_locs(1)＜min_locs(1)：若max_locs(1)大于所述公差值，计算公式为：否则，计算公式为：chg_locs＝max_locs(1)；或者，

当max_locs(1)＞min_locs(1)时：若min_locs(1)大于所述公差值，计算公式为：否则，计算公式为：chg_locs＝min_locs(1)；或者，

s3b、当所述第一个极小值存在，但所述第一个极大值不存在时：若min_locs(1)大于所述公差值，计算公式为：否则，计算公式为：chg_locs＝min_locs(1)；或者，

s3c、当所述第一个极大值存在，但所述第一个极小值不存在时：若max_locs(1)大于所述公差值，计算公式为：否则，计算公式为：chg_locs＝max_locs(1)；或者，

s3d、当所述第一个极大值和第一个极小值均不存在时，定义为：chg_locs＝3。

更优选的，所述公差值不小于3，可以为3～5。

优选的，还包括对重构后的信号进行平滑滤波处理的步骤。

进一步，设重构信号极小值中的最大值为lmax，重构信号极大值中的最小值为upmin，所述平滑滤波处理的步骤包括：

依据lmax和upmin分割子带；

对重构后的信号从第二个点开始依次进行逐点平滑处理：设待处理点为sig(i)，其前一个点为sig(i-1)，其后一个点为sig(i+1)，若sig(i)位于lmax和upmin之间且sig(i-1)和sig(i+1)中的至少一个位于lmax和upmin之间，对sig(i)进行平滑滤波处理；否则，跳过sig(i)，取下一个点按相同方式处理，直至遍历完信号所有的点。

进一步，平滑滤波处理方法为：若sig(i)、sig(i-1)和sig(i+1)均位于lmax和upmin之间，平滑滤波公式为:sig(i)＝(sig(i-1)+sig(i)+sig(i+1))/3；若sig(i)和sig(i-1)均位于lmax和upmin之间且sig(i+1)不在lmax和upmin之间，平滑滤波公式为:sig(i)＝(sig(i-1)+sig(i))/2；若sig(i)和sig(i+1)均位于lmax和upmin之间且sig(i-1)不在lmax和upmin之间，平滑滤波公式为:sig(i)＝(sig(i)+sig(i+1))/2。

本发明一种实施例的有益效果：

该方法根据小波能量选择出需要进行阈值去噪的小波系数，其他小波系数保持不变，与以往采用自适应阈值技术的方法相比，计算量有了显著的减少，提高了去噪效率。进一步对信号进行平滑滤波处理，能改善心电信号的质量。测试结果表明，该方法去噪效果优异，具有较高信噪比。

附图说明

图1为本发明具体实施例提供的心电信号去噪算法总流程图；

图2为本发明具体实施例提供的子带平滑滤波处理方法流程图；

图3为本发明具体实施例对mit-bih数据库中114号记录的去噪效果图；

图4对本发明具体实施例对mit-bih数据库中113号记录加入50hz工频干扰后的信号的去噪效果图。

具体实施方式

下面通过附图及具体实施方式对本发明进行详细的说明。

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。

图1为本发明具体实施例的一种基于小波能量的心电信号去噪算法的流程图，包括：

步骤一：信号均值化处理及公差值的确定。

在对信号进行小波分解之前，首先对信号进行零均值化处理，目的是去掉信号中的直流成分，使得零频信号幅值很小，在后续的功率谱分析中可以减少直流能量对低中频信号谱分析的影响。零均值化处理方法为：对输入信号x，进行x＝x-mean(x),mean(x)代表x的均值。

所述公差值用于在后续的处理中判断去噪的层数。

步骤二：根据matlab中的wmaxlev函数来求出最大分解层数，对均值化处理后的信号进行小波分解，实际的分解层数宜不小于该最大分解层数。为了获得更多的小波系数，可以适当增加分解层数，但不宜过高，否则易造成假频。

本实施例中，小波分解中使用的小波基是db5，与其他小波基相比，db5与心电信号的形态相似，计算结果最佳。

步骤三：根据小波分解后的各层细节系数能量，计算出小波能量曲线的极大和极小点位置，以此来判断出需要去噪的细节阶数变化点。变化点之前的层次认为需要进行去噪，变化点及之后所在的层次，认为无需去噪。

具体步骤如下：

(1)根据matlab中自带的detcoef函数，计算出各层细节系数，接着求出各层细节系数包含的能量大小，计算公式为：其中ej表示第j层的细节系数能量，n表示该层细节系数个数，cofi表示该层的第i个细节系数。

(2)根据计算出来的各层细节系数能量，可以得到小波能量曲线，接着计算出该能量曲线除边界点以外的第一个极大值和第一个极小值的所在层数，以此判断需要去噪的细节阶数变化点。

设待求的所述细节阶数变化点所在层数为chg_locs，小波能量曲线第一个极大值所在层数为max_locs(1)，第一个极小值所在层数为min_locs(1)，所述细节阶数变化点的判断分为以下几种情形：

a.当所述第一个极大值和第一个极小值同时存在时：

a1、当max_locs(1)＜min_locs(1)时：若max_locs(1)大于所述公差值，计算公式为：否则，计算公式为：chg_locs＝max_locs(1)；或者，

a2、当max_locs(1)＞min_locs(1)时：若min_locs(1)大于所述公差值，计算公式为：否则，计算公式为：chg_locs＝min_locs(1)；或者，

b.当所述第一个极小值存在，但所述第一个极大值不存在时：若min_locs(1)大于所述公差值，计算公式为：否则，计算公式为：chg_locs＝min_locs(1)；或者，

c.当所述第一个极大值存在，但所述第一个极小值不存在时：若max_locs(1)大于所述公差值，计算公式为：否则，计算公式为：chg_locs＝max_locs(1)；或者，

d.当所述第一个极大值和第一个极小值均不存在时，定义为：chg_locs＝3。

其中，所述公差值在步骤一中设定，发明人发现，对于该方法，公差值不宜小于3，因为当公差值过小时，得到的去噪层数过少，去噪效果不明显。对于一般的信号，随着公差值增加到一定值，去噪层数不再随着公差值的增大而增加，公差值取3～5即可。当信号的噪声水平非常高时，可以适当增加公差值。

步骤四：将细节阶数变化点之前的各层细节系数进行阈值去噪，根据细节阶数变化点所在层及之后层的细节系数、最后一层的近似系数以及经过阈值去噪后的各层细节系数对信号进行小波重构。

去噪方法可以采用本领域公知的软阈值或硬阈值法，其中，阈值计算公式为：其中t为阈值，n为样本数，σ为噪声标准差。

步骤五：为了进一步减少噪声的影响，还可以对重构后的信号进行平滑滤波处理，处理方法的一种实例如下：

设重构信号极小值中的最大值为lmax，重构信号极大值中的最小值为upmin，依据lmax和upmin分割子带；

对重构后的信号从第二个点开始依次进行逐点平滑处理：设待处理点为sig(i)，其前一个点为sig(i-1)，其后一个点为sig(i+1)，若sig(i)位于lmax和upmin之间且sig(i-1)和sig(i+1)中的至少一个位于lmax和upmin之间，对sig(i)进行平滑滤波处理；否则，跳过sig(i)，取下一个点按相同方式处理，直至遍历完信号所有的点。具体分为以下几种情况：

当sig(i)、sig(i-1)和sig(i+1)均位于lmax和upmin之间时，平滑滤波公式为:sig(i)＝(sig(i-1)+sig(i)+sig(i+1))/3；

当sig(i)和sig(i-1)均位于lmax和upmin之间且sig(i+1)不在lmax和upmin之间时，平滑滤波公式为:sig(i)＝(sig(i-1)+sig(i))/2；

当sig(i)和sig(i+1)均位于lmax和upmin之间且sig(i-1)不在lmax和upmin之间时，平滑滤波公式为:sig(i)＝(sig(i)+sig(i+1))/2；

当sig(i-1)和sig(i+1)均不位于lmax和upmin之间时，平滑滤波公式为:sig(x)＝sig(x)；

当信号sig(i)不位于lmax和upmin之间时，平滑滤波公式为:sig(x)＝sig(x)。

在实际应用中，可采用如图2所示的方法进行平滑滤波处理。

为证明本实施例的去噪效果，选取mit-bih数据库中114号记录进行去噪，去噪过程中，设定公差值设为5，计算的最大分解层数为8，实际分解层数为11，采用软阈值小波去噪及平滑滤波处理，去噪效果如图3，可以看出，小波阈值去噪后的信号与原信号相比，明显减少了大量噪声，而进一步平滑滤波后的信号去除了更多的噪声。

对mit-bih数据库的113号记录加入50hz工频干扰，然后对干扰后的信号进行去噪，设定公差值为5，计算的最大分解层数为8，实际分解层数为11，采用软阈值小波去噪及平滑滤波处理，去噪效果如图4所示。可以看出，本方法去噪效果明显且原始信号的失真较小。

另外对添加了5db白高斯噪声的mit-bih数据库的部分信号进行去噪，利用改善信噪比(snrimp)分析去噪信号的性能，研究重构信号的质量。改善信噪比的公式定义如下:

其中x(i)为原始信号，y(i)为添加噪声后的信号，s(i)为最终去噪并平滑后的信号。

将本实施例的snrimp与自适应相关小波阈值技术(adwt)和多自适应仿生小波变换(mabwt)的snrimp进行比较，结果如表1所示。可以看出，本文提出的方法改善信噪比较高。同时由于本方案小波阈值去噪只使用了部分细节系数，与以往的小波阈值去噪方法相比，计算量有了显著的减少，提高了去噪效率。

表1

以上实施例是对本发明的解释，但是，本发明并不局限于上述实施方式中的具体细节，本领域的技术人员在本发明的技术构思范围内进行的多种等同替代或简单变型方式，均应属于本发明的保护范围。