说明书一种心电信号降噪方法
技术领域
本发明属于医学信号降噪领域，更具体地，涉及一种心电信号降噪方法。
背景技术
心电图为心血管疾病的诊断提供了一种无创、安全、高效的手段。但是,心电信号极其微弱，在采集过程中容易受到基线漂移、工频信号、肌电信号和运动伪影等多种噪声的干扰，这些噪声会降低心电信号的信噪比，不利于疾病的诊断。心电信号降噪对提升心电信号质量，从而提高心血管疾病诊断的正确率具有理论价值和实际意义。
现有的心电信号滤波算法主要包括有限冲击响应(Finite Impulse Response,FIR)滤波、无限冲击响应(Infinite Impulse Response,IIR)滤波、中值滤波、小波滤波和非局部均值降噪等。其中非局部均值降噪在细节保留和背景噪声滤除取得了较好的平衡性，降噪效果较好。
然而现有的心电信号非局部均值降噪方法，不能较好的适应心电信号具有明显周期性和区域性的特点，对不同的心电信号波段，都采用相同的衰减系数，导致难以兼顾均匀区域的平滑和细节信息的保护，滤波效果扔有待提高。
发明内容
针对现有技术的以上缺陷或改进需求，本发明提供了一种心电信号降噪方法，其目的在于在现有非局部均值降噪方法的基础上，对不同的心电信号波段，自适应调整衰减系数，由此解决现有心电信号滤波算法滤波效果不佳的技术问题。
为实现上述目的，按照本发明的一个方面，提供了一种心电信号降噪方法，其特征在于，包括以下步骤：
(1)对原始心电信号O进行均值滤波，计算每一信号点的梯度值gi及衰减系数hi，并判断该信号点为边缘信号点或非边缘信号点；
(1-1)对原始心电信号O进行均值滤波，得到均值滤波心电信号I；
(1-2)对于均值滤波心电信号I中的每一点i，计算其梯度值gi并判断该点是否为边缘区域信号点，并根据判断结果计算衰减系数hi，其中：
梯度值gi为其左邻域信号点与该信号点差值的平均值；
判断该点是否为边缘区域信号点，方法如下：根据预先设定的梯度阈值T，判断该信号点为边缘区信号点或非边缘区信号点，当gi大于T，判断该信号点为边缘区信号点，否则判断该信号点位非边缘区信号点；
根据判断结果计算衰减系数hi，按照式(I)计算：
hi=h01+(gigmax)2gi>Th0gi≤T---(I)]]>
其中gmax是预先设定的梯度值的最大值，h0为固定衰减系数；
(2)对原始心电信号O的进行非局部均值降噪，得到非局部均值预滤波结果M，其中每一点Mi按照式(II)计算，
Mi=Σj∈Ωω(i,j)·Oj---(II)]]>
其中，Ω为搜索窗，ω(i,j)为信号块与相似窗权重，采用式(III)计算：
ω(i,j)=1Ciexp(-G*||O(Ni)-O(Nj)||22hi2)---(III)]]>
其中，G是高斯核函数，Ci为归一化参数，O(Ni)为信号块，O(Nj)为相似 窗，hi为步骤(1)获得的衰减系数；
(3)对步骤(2)中获得的非局部均值预滤波结果M，根据步骤(1)获得的梯度值gi结果修正，得到降噪后的心电信号；
(3-1)对于非局部均值预滤波结果M每一点Mi，计算其梯度值g′i，梯度值g′i为其左邻域信号点与该信号点差值的平均值，根据预先设定的梯度阈值T，判断该信号点为边缘区信号点或非边缘区信号点；如果判断结果与根据步骤(1-2)判断的结果相同，则该点不需要修正，将Mi作为降噪后的心电信号，否则执行步骤(3-2)；
(3-2)对非局部均值预滤波结果中需要修正的点Mi进行修正：将式(I)中的gi替换为g′i，得到该点的修正后的衰减参数h′i，将式(III)中的hi替换为h′i，得到该点修正后的信号块与相似窗权重ω(i,j)'，将式(II)中的ω(i,j)替换为ω(i,j)'，得到该点修正后的非局部均值预滤波结果Mi'，将Mi'作为降噪后的心电信号。
优选地，所述心电信号降噪方法，其步骤(1-1)所述的均值滤波窗口大小在3×3至6×6之间，优选5×5。
优选地，所述心电信号降噪方法，步骤(1-2)所述的gmax在15至20之间，优选16。
优选地，所述心电信号降噪方法，步骤(1-2)所述的T在9至14之间，优选10。
优选地，所述心电信号降噪方法，其步骤(1-2)所述左领域为邻域信号点i-3,i-2和i-1。
优选地，所述心电信号降噪方法，其步骤(2)所述搜索窗的大小为7×7。
优选地，所述心电信号降噪方法，其步骤(2)所述信号块与相似窗大小相同，为3×3。
总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：
本发明提供的心电信号降噪方法，是基于非局部均值的降噪方法，在细节保留和噪声滤除的平衡性方面具有先天优势。同时考虑到心电信号具有明显周期性和区域性的特点，对于心电信号的不同波段区域，自适应的调整非局部均值降噪的核心参数，即衰减系数。相对于现有技术，本发明提供的心电信号的降噪方法，能在有效抑制心电信号中噪声的同时更好地保护信号细节信息。
附图说明
图1是本发明自适应非局部均值心电信号滤波的流程图；
图2是不同噪声水平下各种滤波算法对应仿真心电信号降噪结果的信噪比；
图3是实施例1噪声方差为20时各种滤波算法对应的仿真心电图滤波结果，其中图3(a)是模拟心电信号，图3(b)是原始信号，图3(c)实施例1提供的降噪结果，图3(d)是NLM方法降噪结果，图3(e)是IIR滤波方法降噪结果，图3(f)是FIR滤波方法降噪结果，图3(g)是中值滤波方法降噪结果，图3(h)是小波变换滤波方法降噪结果；
图4是实施例1噪声方差为60时各种滤波算法对应的仿真心电图滤波结果，其中图4(a)是模拟心电信号，图4(b)是原始信号，图4(c)实施例1提供的降噪结果，图4(d)是NLM方法降噪结果，图4(e)是IIR滤波方法降噪结果，图4(f)是FIR滤波方法降噪结果，图4(g)是中值滤波方法降噪结果，图4(h)是小波变换滤波方法降噪结果；
图5是实施例2各种滤波算法对实际心电信号进行降噪的结果，其中图5(a)是真实心电信号，图5(b)实施例2提供的降噪结果，图5(c)是NLM方法降噪结果，图5(d)是IIR滤波方法降噪结果，图5(e)是FIR滤波方法降噪结果，图5(f)是中值滤波方法降噪结果，图5(g)是 小波变换滤波方法降噪结果。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
本发明提供的心电信号降噪方法，简称ANLM，如图1所示，包括以下步骤：
(1)对原始心电信号O进行均值滤波，计算每一信号点的梯度值gi及衰减系数hi，并判断该信号点为边缘信号点或非边缘信号点；
(1-1)对原始心电信号O进行均值滤波，得到均值滤波心电信号I；考虑到降噪速度和降噪效果，所述的均值滤波窗口大小在3×3至6×6之间，优选5×5。
(1-2)对于均值滤波心电信号I中的每一点i，计算其梯度值gi并判断该点是否为边缘区域信号点，并根据判断结果计算衰减系数hi，其中：
梯度值gi为其左邻域信号点与该信号点差值的平均值；所述左领域为邻域信号点i-3,i-2和i-1。
判断该点是否为边缘区域信号点，方法如下：根据预先设定的梯度阈值T，判断该信号点为边缘区信号点或非边缘区信号点，当gi大于T，判断该信号点为边缘区信号点，否则判断该信号点位非边缘区信号点；
根据判断结果计算衰减系数hi，按照式(I)计算：
hi=h01+(gigmax)2gi>Th0gi≤T---(I)]]>
其中gmax是预先设定的梯度值的最大值，h0为固定衰减系数。
gmax在15至20之间，优选16；T在9至14之间，优选10；h0取值为0.1至25之间。
h0按照如下方法确定：
h0由经验公式h0＝Cσ2得到，C是常数，取值范围是0.3到1.5，σ2是噪声的方差，取值优选在1至60之间，根据信号质量越差噪声方差越大的原则来选取。
(2)对原始心电信号O的进行非局部均值降噪，得到非局部均值预滤波结果M，其中每一点Mi按照式(II)计算，
Mi=Σj∈Ωω(i,j)·Oj---(II)]]>
其中，Ω为搜索窗，ω(i,j)为信号块与相似窗权重，采用式(III)计算：
ω(i,j)=1Ciexp(-G*||O(Ni)-O(Nj)||22hi2)---(III)]]>
其中，G是高斯核函数，Ci为归一化参数，O(Ni)为信号块，O(Nj)为相似窗，hi为步骤(1)获得的衰减系数，*代表卷积运算，为L2范数，Ci为归一化参数以确保成立，并采用下式计算
Ci=Σj∈Ωexp(-G*||O(Ni)-O(Nj)||22hi2)]]>
为兼顾降噪时间和降噪效果，优选地，所述搜索窗的大小为7×7，所述 信号块与相似窗大小相同，为3×3。
(3)对步骤(2)中获得的非局部均值预滤波结果M，根据步骤(1)获得的梯度值gi结果修正，得到降噪后的心电信号；
(3-1)对于非局部均值预滤波结果M每一点Mi，计算其梯度值g′i，梯度值g′i为其左邻域信号点与该信号点差值的平均值，根据预先设定的梯度阈值T，判断该信号点为边缘区信号点或非边缘区信号点；如果判断结果与根据步骤(1-2)判断的结果相同，则该点不需要修正，将Mi作为降噪后的心电信号，否则执行步骤(3-2)；
(3-2)对非局部均值预滤波结果中需要修正的点Mi进行修正：将式(I)中的gi替换为g′i，得到该点的修正后的衰减参数h′i，将式(III)中的hi替换为h′i，得到该点修正后的信号块与相似窗权重ω(i,j)'，将式(II)中的ω(i,j)替换为ω(i,j)'，得到该点修正后的非局部均值预滤波结果Mi'，将Mi'作为降噪后的心电信号。
以下为实施例：
实施例1
对图3(a)、图4(a)所示的模拟心电信号，我们分别加入方差为20、30、40、50、60的高斯噪声，得到原始信号O，如图3(b)、图4(b)所示，采用本专利提出的心电降噪方法进行降噪处理的步骤如下：
(1)对原始心电信号O进行均值滤波，计算每一信号点的梯度值gi及衰减系数hi，并判断该信号点为边缘信号点或非边缘信号点；
(1-1)对原始心电信号O进行均值滤波，得到均值滤波心电信号I；考虑到降噪速度和降噪效果，所述的均值滤波窗口大小为5×5。
(1-2)对于均值滤波心电信号I中的每一点i，计算其梯度值gi并判断该点是否为边缘区域信号点，并根据判断结果计算衰减系数hi，其中：
梯度值gi为其左邻域信号点与该信号点差值的平均值；所述左领域为邻域信号点i-3,i-2和i-1。
判断该点是否为边缘区域信号点，方法如下：根据预先设定的梯度阈值T，判断该信号点为边缘区信号点或非边缘区信号点，当gi大于T，判断该信号点为边缘区信号点，否则判断该信号点位非边缘区信号点；
根据判断结果计算衰减系数hi，按照式(I)计算：
hi=h01+(gigmax)2gi>Th0gi≤T---(I)]]>
其中gmax是预先设定的梯度值的最大值，h0为固定衰减系数。
gmax选取为16；T选取为10；h0对于仿真心电信号为25。
(2)对原始心电信号O的进行非局部均值降噪，得到非局部均值预滤波结果M，其中每一点Mi按照式(II)计算，
Mi=Σj∈Ωω(i,j)·Oj---(II)]]>
其中，Ω为搜索窗，ω(i,j)为信号块与相似窗权重，采用式(III)计算：
ω(i,j)=1Ciexp(-G*||O(Ni)-O(Nj)||22hi2)---(III)]]>
其中，G是高斯核函数，Ci为归一化参数，O(Ni)为信号块，O(Nj)为相似窗，hi为步骤(1)获得的衰减系数，*代表卷积运算，为L2范数，Ci为归一化参数以确保成立，并采用下式计算
Ci=Σj∈Ωexp(-G*||O(Ni)-O(Nj)||22hi2)]]>
为兼顾降噪时间和降噪效果，优选地，所述搜索窗的大小为7×7，所述信号块与相似窗大小相同，为3×3。
(3)对步骤(2)中获得的非局部均值预滤波结果M，根据步骤(1)获得的梯度值gi结果修正，得到降噪后的心电信号；
(3-1)对于非局部均值预滤波结果M每一点Mi，计算其梯度值g′i，梯度值g′i为其左邻域信号点与该信号点差值的平均值，根据预先设定的梯度阈值T，判断该信号点为边缘区信号点或非边缘区信号点；如果判断结果与根据步骤(1-2)判断的结果相同，则该点不需要修正，将Mi作为降噪后的心电信号，否则执行步骤(3-2)；
(3-2)对非局部均值预滤波结果中需要修正的点Mi进行修正：将式(I)中的gi替换为g′i，得到该点的修正后的衰减参数h′i，将式(III)中的hi替换为h′i，得到该点修正后的信号块与相似窗权重ω(i,j)'，将式(II)中的ω(i,j)替换为ω(i,j)'，得到该点修正后的非局部均值预滤波结果Mi'，将Mi'作为降噪后的心电信号。
为了说明本专利提出的ANLM方法的优越性，我们将传统NLM方法、IIR滤波、FIR滤波、中值滤波、小波变换滤波这五种算法作为对比算法。这些方法采用的参数如下：传统NLM方法中，搜索窗和相似窗的大小分别为7×7和3×3，方法中衰减系数h0为25；FIR和IIR滤波，我们对其进行了迭代处理，共运行算法两次；中值滤波窗口选为5×5；小波变换层数为5，选择的基波为sym8，其阈值采用Minimax阈值获得。
为客观评价各种方法的滤波性能，采用信噪比(SNR)和均方误差(RMSE)作为评价指标，这两个指标其定义如下：
SNR=10log(I2(I-In)2)]]>
其中I为未被噪声污染的ECG信号，In为加高斯噪声的ECG信号。
RMSE=Σ(I-In)2n]]>
其中n为未加噪声ECG信号的长度。
表1和图2分别给出了各种滤波算法对应仿真心电信号降噪结果的均方误差和信噪比。从表1和图2可看出，本专利提出的ANLM方法较其它比较方法具有更低的RMSE和更高的SNR，这表明ANLM方法在所有的比较方法中具有最优的恢复性能。
表1不同噪声水平下各种滤波算法对应仿真心电信号降噪结果的均方误差
从图3、图4可以看出：IIR、FIR和中值滤波对应的心电降噪信号中残留了较多的噪声，如图3(e)、图3(f)、图3(h)、图4(e)、图4(f)、图4(h)长方形窗口所示；NLM和小波算法对心电图中的细节造成一定的破坏，如图3(d)、图3(g)、图4(d)、图4(g)正方形窗口所示。相比之下，我们提出的ANLM算法既能很好地保护边缘等信号细节，又能有效地滤除平滑区的噪声，具有较其它比较算法更优的恢复性能。
实施例2
为了说明本算法的实用性，我们将其用于对儿科头部受伤病人在500Hz采样率下采集到的实际心电图的降噪处理，采用本专利提出的ANLM方法进行降噪处理的步骤如下：
(1)对原始心电信号O进行均值滤波，计算每一信号点的梯度值gi及衰减系数hi，并判断该信号点为边缘信号点或非边缘信号点；
(1-1)对原始心电信号O进行均值滤波，得到均值滤波心电信号I；考虑到降噪速度和降噪效果，所述的均值滤波窗口大小选5×5。
(1-2)对于均值滤波心电信号I中的每一点i，计算其梯度值gi并判断该点是否为边缘区域信号点，并根据判断结果计算衰减系数hi，其中：
梯度值gi为其左邻域信号点与该信号点差值的平均值；所述左领域为邻域信号点i-3,i-2和i-1。
判断该点是否为边缘区域信号点，方法如下：根据预先设定的梯度阈值T，判断该信号点为边缘区信号点或非边缘区信号点，当gi大于T，判断该信号点为边缘区信号点，否则判断该信号点位非边缘区信号点；
根据判断结果计算衰减系数hi，按照式(I)计算：
hi=h01+(gigmax)2gi>Th0gi≤T---(I)]]>
其中gmax是预先设定的梯度值的最大值，h0为固定衰减系数。
gmax选取为16；T选取为10；h0对于临床心电信号为0.1。
(2)对原始心电信号O的进行非局部均值降噪，得到非局部均值预滤波结果M，其中每一点Mi按照式(II)计算，
Mi=Σj∈Ωω(i,j)·Oj---(II)]]>
其中，Ω为搜索窗，ω(i,j)为信号块与相似窗权重，采用式(III)计算：
ω(i,j)=1Ciexp(-G*||O(Ni)-O(Nj)||22hi2)---(III)]]>
其中，G是高斯核函数，Ci为归一化参数，O(Ni)为信号块，O(Nj)为相似窗，hi为步骤(1)获得的衰减系数，*代表卷积运算，为L2范数，Ci为归一化参数以确保成立，并采用下式计算
Ci=Σj∈Ωexp(-G*||O(Ni)-O(Nj)||22hi2)]]>
为兼顾降噪时间和降噪效果，优选地，所述搜索窗的大小为7×7，所述信号块与相似窗大小相同，为3×3。
(3)对步骤(2)中获得的非局部均值预滤波结果M，根据步骤(1)获得的梯度值gi结果修正，得到降噪后的心电信号；
(3-1)对于非局部均值预滤波结果M每一点Mi，计算其梯度值g′i，梯度值g′i为其左邻域信号点与该信号点差值的平均值，根据预先设定的梯度阈值T，判断该信号点为边缘区信号点或非边缘区信号点；如果判断结果与根据步骤(1-2)判断的结果相同，则该点不需要修正，将Mi作为降噪后的心电信号，否则执行步骤(3-2)；
(3-2)对非局部均值预滤波结果中需要修正的点Mi进行修正：将式(I)中的gi替换为g′i，得到该点的修正后的衰减参数h′i，将式(III)中的hi替换为h′i，得到该点修正后的信号块与相似窗权重ω(i,j)'，将式(II)中的ω(i,j)替换为ω(i,j)'，得到该点修正后的非局部均值预滤波结果Mi'，将Mi'作为降噪后的心电信号。
为了说明本专利提出的ANLM方法的优越性，我们将传统NLM方法、FIR滤波、IIR滤波、中值滤波、小波变换滤波这五种算法作为对比算法，这些方法采用的参数和实施例1中相同。从图5对比可以明显看出：ANLM滤波后的实际心电图更平滑，在边缘区保留了最多的细节信息。上述对比 证实了本专利提出的算法在实际心电图降噪中的有效性和优越性。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。