CN105046654B - 一种基于粒子群优化的心电信号自适应非局部均值降噪方法 - Google Patents

一种基于粒子群优化的心电信号自适应非局部均值降噪方法 Download PDF

Info

Publication number
CN105046654B
CN105046654B CN201510348834.6A CN201510348834A CN105046654B CN 105046654 B CN105046654 B CN 105046654B CN 201510348834 A CN201510348834 A CN 201510348834A CN 105046654 B CN105046654 B CN 105046654B
Authority
CN
China
Prior art keywords
signal
particle
search window
formula
fit
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201510348834.6A
Other languages
English (en)
Other versions
CN105046654A (zh
Inventor
张旭明
张靖阳
李永帅
任金霞
李想
王顺顺
周欢
丁明跃
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Huazhong University of Science and Technology
Original Assignee
Huazhong University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Huazhong University of Science and Technology filed Critical Huazhong University of Science and Technology
Priority to CN201510348834.6A priority Critical patent/CN105046654B/zh
Publication of CN105046654A publication Critical patent/CN105046654A/zh
Application granted granted Critical
Publication of CN105046654B publication Critical patent/CN105046654B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
  • Soundproofing, Sound Blocking, And Sound Damping (AREA)

Abstract

本发明公开了基于粒子群优化的心电信号自适应非局部均值降噪方法，包括以下步骤：(1)对原始心电信号进行非局部均值滤波，得到滤波结果；(2)采用粒子群算法对每一信号点i自适应确定非局部均值滤波方法中的两个关键参数，即衰减参数和搜索窗口；(3)利用自适应确定的上述两个参数，对原始心电信号进行非局部均值滤波，获得最终的降噪信号。

Description

一种基于粒子群优化的心电信号自适应非局部均值降噪方法
技术领域
本发明属于医学信号降噪领域，更具体地，涉及一种医学超声图像降噪方法。
背景技术
心电信号采集自病人体表，是一种无创性的检测手段。因此，心电信号采集过程中受到基线漂移、工频干扰、肌电信号等多种噪声的干扰，影响后续临床诊断的准确性。因此这些噪声能够被适当的去除，对心电信号诊断具有决定性的意义。
现有的心电信号滤波算法主要包括有限冲击响应(Finite Impulse Response,FIR)滤波、无限冲击响应(Infinite Impulse Response,IIR)滤波、中值滤波、小波滤波和非局部均值降噪等。其中，非局部均值方法因充分利用了信号的冗余性，较其它方法能更好地抑制噪声的不利影响。但传统非局部均值降噪方法，对不同的心电信号波段，都采用相同的衰减系数，不能较好地适应心电信号具有明显周期性和区域性的特点，导致难以兼顾均匀区域的平滑和细节信息的保护。
发明内容
针对现有技术中存在的问题，本申请提供的是一种基于粒子群优化的心电信号自适应非局部均值降噪方法，简称ANLM。其中通过对衰减参数和搜索窗口进行研究和设计，与现有降噪方法相比，ANLM具有更好的恢复性能。
为实现上述目的，按照本发明的一个方面，提供了一种基于粒子群优化算法的心电信号自适应非局部均值降噪方法,该方法包括以下步骤：
步骤1对原始心电信号进行非局部均值滤波，得到预滤波结果。
对原始信号O进行非局部均值滤波，得到预滤波信号M，该信号中每个信号点Mi按式(1)计算
其中，Oj代表原始信号O中位于点j处信号的幅值，Ω表示搜索窗，ω(i,j)是衡量原始信号O中分别以i和j为中心的两信号块O(Ni)和O(Nj)之间的相似度，其计算方法如(2)式所示：
其中，O(Ni)和O(Nj)分别代表信号O中分别以i和j为中心的两信号块，Ci为归一化参数，表示L2范数，G表示高斯内核，*表示卷积操作符号，h代表固定的衰减参数。
步骤2利用粒子群算法来对每一信号点i优化衰减参数和搜索窗口,将信号的偏差的平方和其方差的和作为适应度函数，使其自适应于不同的信号点，进而得到每一个信号点对应的最优化的衰减参数和搜索窗口。
(2-1)初始化粒子群。
以搜索窗口Ωi和衰减参数hi为优化对象，在二维目标搜索空间中初始化N个粒子组成一个种群，其中第k个粒子表示为一个二维的位置向量xk＝(Sk,hk),(k＝1,2,…,N)，Sk为搜索窗Ωi的半径。而第k个粒子的移动速度对应的二维的向量为vk＝(vk1,vk2),(k＝1,2,…,N)。
(2-2)对于信号O，用逐信号点的均方误差的期望E[MSE(Mi)]作为优化衰减参数hi和搜索窗Ωi的适应度函数去计算每个粒子的适应度Fit(k)。考虑到均方误差的期望E[MSE(Mi)]越小，即表示对应粒子的搜索窗大小Ωi和衰减参数hi降噪效果越好。故适应度Fit(k)采用下式(3)计算：
Fit(k)＝-E[MSE(Mi)] (3)
其中均方误差的期望E[MSE(Mi)]如下式(4)计算：
E[MSE(Mi)]＝bias2(Mi)+var(Mi) (4)
其中偏差平方项bias2(Mi)和方差项var(Mi)计算公式分别如式(5)和式(6)所示：
其中Mi是预滤波信号M中信号点i对应的幅度值,σ为信号O噪声标准差。
(2-3)用下式表示每个粒子迄今为止搜索到的局部最优位置pbest，并利用(3)式计算局部最优位置的适应度即个体最优fpbest(k)：
pbest＝(Sk,hk),(k＝1,2,…,N) (7)
用对应每个粒子的适应度Fit(k)和个体极值fpbest(k)比较，如果Fit(k)大于fpbest(k)，则用Fit(k)替换掉fpbest(k)；
(2-4)用下式表示整个粒子群迄今为止搜索到的全局最优位置gbest，并利用(3)式计算其适应度即全局极值fgbest：
gbest＝(Sg,hg) (8)
对每个粒子，用它的适应度值Fit(k)和全局极值fgbest比较，如果Fit(k)大于fgbest，则用Fit(k)替换掉fgbest；
(2-5)根据下式(9)和(10)更新粒子的速度和位置；
vk＝T*vk+c1r1(pbest-xk)+c2r2(gbest-xk) (9)
xk＝xk+vk (10)
其中c1和c2是学习因子，r1和r2是介于[0,1]之间的随机数，T是惯性权重。
如果满足结束条件(误差足够好或到达最大循环次数)退出，否则返回(2-2)。
(2-6)粒子群优化完之后，会得到群体最优的位置gbest＝(Sg,hg)，于是就得到了逐信号点最优化的搜索窗S＝Sg，最优化的衰减参数h＝hg。
步骤3用步骤2粒子群算法所得逐信号最优的搜索窗半径Sg和衰减参数hg对原始信号O进行逐信号点的非局部均值降噪，得到降噪之后的信号M′。该降噪信号中每一信号点i对应的幅度Mi′按式(11)计算：
其中Ωi′表示信号中以信号点i为中心的搜索窗，其窗口大小为步骤3得到的对应信号点i最优的Sg。ω′(i,j)是衡量预处理信号M中分别以i和j为中心的两信号块O′(Ni)和O′(Nj)之间的相似度，其计算公式如下(12)：
O'(Ni)和O'(Nj)分别代表预处理信号M中以信号点i和j为中心的信号块，hg是对应信号点i最优的衰减参数。
优选地，所述心电信号降噪方法，步骤1非局部均值预滤波搜索窗Ω窗口大小在800至1200之间，衰减参数h在15至20之间。为了实现较好的噪声抑制效果，得到噪声少的预滤波信号，搜索窗Ω窗口大小优选1000，衰减参数h优选为16。
优选地，所述心电信号降噪方法，步骤(2-1)所述的N取值在15至30之间，为了兼顾粒子群算法优化效果和时间效率,优选N为20。
优选地，所述心电信号降噪方法，其特征在于，考虑到心电信号的细节尖峰区和信号平坦区的特征，步骤(2-1)所述搜索窗大小Sk的搜索范围为500-2000，衰减参数hk的搜索范围为5-30。
优选地，所述心电信号降噪方法，步骤(2-5)考虑到粒子群算法优化迭代的效果和效率，所述惯性权重优选T＝0.7298，学习因子c1和c2优选为c1＝c2＝1.4962，最大迭代次数优选为50。
优选地，为了提高步骤3自适应逐信号非局部均值降噪效果，得到噪声较少的最终降噪效果，所述信号块O′(Ni)和O′(Nj)大小优选10。
总体而言，按照本发明的上述技术构思与现有技术相比，主要具备以下的技术优点：
1、在细节保留和噪声滤除的平衡性方面具有先天优势，同时考虑到心电信号具有明显周期性和区域性的特点，对于心电信号的不同波段区域，利用粒子群算法自适应的调整非局部均值降噪的核心参数，即衰减系数和搜索窗大小；
2、能在有效抑制心电信号中噪声的同时更好地保护信号细节信息。
附图说明
图1是本发明基于粒子群优化的心电信号自适应非局部均值降噪方法的流程图；
图2是本发明粒子群算法步骤2的流程图；
图3(a)是模拟心电信号；
图3(b)是噪声方差为60时的原始信号；
图3(c)是对噪声方差为60的原始信号实施例1ANLM方法的降噪结果；
图3(d)是对噪声方差为60的原始信号实施NLM方法的降噪结果；
图3(e)是对噪声方差为60的原始信号实施IIR滤波方法的降噪结果；
图3(f)是对噪声方差为60的原始信号实施FIR滤波方法的降噪结果；
图3(g)是对噪声方差为60的原始信号实施中值滤波方法的降噪结果；
图3(h)是对噪声方差为60的原始信号实施小波变换滤波方法的降噪结果；
图4(a)是真实心电信号；
图4(b)是实施例2 ANLM方法的降噪结果；
图4(c)是实施NLM方法的降噪结果；
图4(d)是实施IIR滤波方法的降噪结果；
图4(e)是实施FIR滤波方法的降噪结果；
图4(f)是实施中值方法的降噪结果；
图4(g)是实施小波变换滤波方法的降噪结果。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
本发明提供的心电信号降噪方法，简称ANLM，如图1所示，包括以下步骤：
步骤1对原始心电信号进行非局部均值滤波，得到滤波结果。
对原始信号O进行非局部均值滤波，得到预滤波信号M，该信号中每个信号点Mi按式(1)计算
其中，Oj代表原始信号O中位于点j处信号的幅值，Ω表示搜索窗，其大小选择为1000，ωi,j是衡量信号O中分别以i和j为中心的两信号块O(Ni)和O(Nj)之间的相似度，其计算方法如(2)式所示：
其中，Ci为归一化参数，表示L2范数，G表示高斯内核，*表示卷积操作符号，h代表固定的衰减参数，其优选为16。
步骤2利用粒子群算法来对每一信号点i优化衰减参数和搜索窗口,将信号的偏差的平方和其方差的和作为适应度函数，使其自适应于不同的信号点，进而得到每一个信号点对应的最优化的衰减参数和搜索窗口。
(2-1)初始化粒子群。
以搜索窗口Ωi和衰减参数hi为优化对象，在二维目标搜索空间中初始化N个粒子组成一个种群，其中第k个粒子表示为一个二维的位置向量xk＝(Sk,hk),(k＝1,2,…,20)，Sk为搜索窗Ωi的半径。而第k个粒子的移动速度对应的二维的向量为vk＝(vk1,vk2),(k＝1,2,…,20)。
优选的，Sk的搜索范围为[500,2000]，hk的搜索范围为[5,30]。
(2-2)对于信号O，用逐信号点的均方误差的期望E[MSE(Mi)]作为优化衰减参数hi和搜索窗Ωi的适应度函数去计算每个粒子的适应度Fit(k)。考虑到均方误差的期望E[MSE(Mi)]越小，即表示对应粒子的搜索窗大小Ωi和衰减参数hi降噪效果越好。故适应度Fit(k)采用下式(3)计算：
Fit(k)＝-E[MSE(Mi)] (3)
其中均方误差的期望E[MSE(Mi)]如下式(4)计算：
E[MSE(Mi)]＝bias2(Mi)+var(Mi) (4)
其中偏差平方项bias2(Mi)和方差项var(Mi)计算公式分别如式(5)和式(6)所示：
其中Mi是预滤波信号M中信号点i对应的幅度值,σ为信号O噪声标准差。
(2-3)用下式表示每个粒子迄今为止搜索到的局部最优位置pbest，并利用(3)式计算局部最优位置的适应度即个体最优fpbest(k)：
pbest＝(Sk,hk),(k＝1,2,…,20) (7)
用对应每个粒子的适应度Fit(k)和个体极值fpbest(k)比较，如果Fit(k)大于fpbest(k)，则用Fit(k)替换掉fpbest(k)；
(2-4)用下式表示整个粒子群迄今为止搜索到的全局最优位置gbest，并利用(3)式计算其适应度即全局极值fgbest：
gbest＝(Sg,hg) (8)
对每个粒子，用它的适应度值Fit(k)和全局极值fgbest比较，如果Fit(k)大于fgbest，则用Fit(k)替换掉fgbest；
(2-5)根据下式(9)和(10)更新粒子的速度和位置；
vk＝T*vk+c1r1(pbest-xk)+c2r2(gbest-xk) (9)
xk＝xk+vk (10)
其中c1和c2是学习因子，优选为1.4962，r1和r2是介于[0,1]之间的随机数，T是惯性权重，优选为0.7298。
如果满足最大循环次数退出，否则返回(2-2)。其最大循环次数设置为50。
(2-6)粒子群优化完之后，会得到群体最优的位置gbest＝(Sg,hg)，于是就得到了逐信号点最优化的搜索窗S＝Sg，最优化的衰减参数h＝hg。
步骤3用步骤2粒子群算法所得逐信号最优的搜索窗半径Sg和衰减参数hg对原始信号O进行逐信号点的非局部均值降噪，得到降噪之后的信号M′。该降噪信号中每一信号点i对应的幅度Mi′按式(11)计算：
其中Ωi′表示信号中以信号点i为中心的搜索窗，其窗口大小为步骤3得到的对应信号点i最优的Sg。ω′i,j是衡量预处理信号M中分别以i和j为中心的两信号块O′(Ni)和O′(Nj)之间的相似度，其计算公式如下(12)：
O′(Ni)和O′(Nj)分别代表预处理信号M中以信号点i和j为中心的信号块，其信号块大小为10，hg是对应信号点i最优的衰减参数。
以下为实施例：
实施例1
对图3(a)所示的模拟心电信号加入方差为20、30、40、50、60的高斯噪声，得到图3(b)原始信号O,接着采用本专利提出的心电信号降噪方法进行降噪处理：
步骤1对原始心电信号图3(b)进行非局部均值滤波，得到滤波结果。
对原始信号O进行非局部均值滤波，得到预滤波信号M，该信号中每个信号点Mi按式(1)计算
其中，Oj代表原始信号O中位于点j处信号的幅值，Ω表示搜索窗，其大小选择为1000，ωi,j是衡量信号O中分别以i和j为中心的两信号块O(Ni)和O(Nj)之间的相似度，其计算方法如(2)式所示：
其中，Ci为归一化参数，表示L2范数，G表示高斯内核，*表示卷积操作符号，h代表固定的衰减参数，对于模拟信号优选为16。Ci为归一化参数，以确保成立，采用下式计算
步骤2利用粒子群算法来对每一信号点i优化衰减参数和搜索窗口,将信号的偏差的平方和其方差的和作为适应度函数，使其自适应于不同的信号点，进而得到每一个信号点对应的最优化的衰减参数和搜索窗口。
(2-1)初始化粒子群。
以搜索窗口Ωi和衰减参数hi为优化对象，在二维目标搜索空间中初始化N个粒子组成一个种群，其中第k个粒子表示为一个二维的位置向量xk＝(Sk,hk),(k＝1,2,…,20)，Sk为搜索窗Ωi的半径。而第k个粒子的移动速度对应的二维的向量为vk＝(vk1,vk2),(k＝1,2,…,20)。
优选的，Sk的搜索范围为[500,2000]，hk的搜索范围为[5,30]。
(2-2)对于信号O，用逐信号点的均方误差的期望E[MSE(Mi)]作为优化衰减参数hi和搜索窗Ωi的适应度函数去计算每个粒子的适应度Fit(k)。考虑到均方误差的期望E[MSE(Mi)]越小，即表示对应粒子的搜索窗大小Ωi和衰减参数hi降噪效果越好。故适应度Fit(k)采用下式(3)计算：
Fit(k)＝-E[MSE(Mi)] (3)
其中均方误差的期望E[MSE(Mi)]如下式(4)计算：
E[MSE(Mi)]＝bias2(Mi)+var(Mi) (4)
其中偏差平方项bias2(Mi)和方差项var(Mi)计算公式分别如式(5)和式(6)所示：
其中Mi是预滤波信号M中信号点i对应的幅度值,σ为信号O噪声标准差。
(2-3)用下式表示每个粒子迄今为止搜索到的局部最优位置pbest，并利用(3)式计算局部最优位置的适应度即个体最优fpbest(k)：
pbest＝(Sk,hk),(k＝1,2,…,20) (7)
用对应每个粒子的适应度Fit(k)和个体极值fpbest(k)比较，如果Fit(k)大于fpbest(k)，则用Fit(k)替换掉fpbest(k)；
(2-4)用下式表示整个粒子群迄今为止搜索到的全局最优位置gbest，并利用(3)式计算其适应度即全局极值fgbest：
gbest＝(Sg,hg) (8)
对每个粒子，用它的适应度值Fit(k)和全局极值fgbest比较，如果Fit(k)大于fgbest，则用Fit(k)替换掉fgbest；
(2-5)根据下式(9)和(10)更新粒子的速度和位置；
vk＝T*vk+c1r1(pbest-xk)+c2r2(gbest-xk) (9)
xk＝xk+vk (10)
其中c1和c2是学习因子，优选为1.4962，r1和r2是介于[0,1]之间的随机数，T是惯性权重，优选为0.7298。
如果满足最大循环次数退出，否则返回(2-2)。其最大循环次数设置为50。
(2-6)粒子群优化完之后，会得到群体最优的位置gbest＝(Sg,hg)，于是就得到了逐信号点最优化的搜索窗S＝Sg，最优化的衰减参数h＝hg。
步骤3用步骤2粒子群算法所得逐信号最优的搜索窗半径Sg和衰减参数hg对原始信号O进行逐信号点的非局部均值降噪，得到降噪之后的信号M′。该降噪信号中每一信号点i对应的幅度Mi′按式(11)计算：
其中Ωi′表示信号中以信号点i为中心的搜索窗，其窗口大小为步骤3得到的对应信号点i最优的Sg。ω′i,j是衡量预处理信号M中分别以i和j为中心的两信号块O′(Ni)和O′(Nj)之间的相似度，其计算公式如下(12)：
O′(Ni)和O′(Nj)分别代表预处理信号M中以信号点i和j为中心的信号块，其信号块大小为10，hg是对应信号点i最优的衰减参数。
为了说明本专利提出的ANLM方法的优越性，我们将传统NLM方法、IIR滤波、FIR滤波、中值滤波、小波变换滤波这五种算法作为对比算法。
为客观评价各种方法的滤波性能，采用信噪比(SNR)和均方误差(RMSE)作为评价指标，这两个指标其定义如下：
其中I为未被噪声污染的ECG信号，In为加高斯噪声的ECG信号
其中n为未加噪声ECG信号的长度。
表1给出了各种滤波算法对应仿真心电信号降噪结果的均方误差和信噪比。从表1可看出，本专利提出的ANLM方法较其它比较方法具有更低的RMSE和更高的SNR，这表明ANLM方法在所有的比较方法中具有最优的恢复性能。
表1不同噪声水平下各种滤波算法对应仿真心电信号降噪结果的信噪比和均方误差
从图3可以看出：IIR、FIR和中值滤波对应的心电降噪信号中残留了较多的噪声，如图3(e)、图3(f)、图3(g)中包含平滑区的长方形窗口所示；NLM和小波算法对心电图中的细节造成一定的破坏，如图3(d)、图3(h)中包含尖峰的正方形窗口所示。相比之下，我们提出的ANLM算法既能很好地保护边缘等信号细节，又能有效地滤除平滑区的噪声，具有较其它比较算法更优的恢复性能。
实施例2
为了说明本算法的实用性，我们将其用于对儿科头部受伤病人在500Hz采样率下采集到的实际心电图4(a)的降噪处理，采用本专利提出的ANLM方法进行降噪处理的步骤如下：
步骤1对图4(a)原始心电信号进行噪声评估，得到噪声方差。
对原始信号O进行噪声评估，该信号的噪声方差σ′按下式计算
σ′＝1.4826*median(|R-median(R)|)
其中R＝{R1,…,Rx}是通过手动选择的噪声信号O中同质区域的局部残差集，并且
步骤2对实际心电信号图4(a)进行非局部均值滤波，得到滤波结果。
对图4(a)进行非局部均值滤波，得到预滤波信号M，该信号中每个信号点Mi按式(1)计算
其中，Oj代表原始信号O中位于点j处信号的幅值，Ω表示搜索窗，其大小选择为1000，ωi,j是衡量信号O中分别以i和j为中心的两信号块O(Ni)和O(Nj)之间的相似度，其计算方法如(2)式所示：
其中，Ci为归一化参数，表示L2范数，G表示高斯内核，*表示卷积操作符号，h代表固定的衰减参数，对于实际信号其优选为16。Ci为归一化参数，以确保成立，采用下式计算
步骤3利用粒子群算法来对每一信号点i优化衰减参数和搜索窗口,将信号的偏差的平方和其方差的和作为适应度函数，使其自适应于不同的信号点，进而得到每一个信号点对应的最优化的衰减参数和搜索窗口。
(3-1)初始化粒子群。
以搜索窗口Ωi和衰减参数hi为优化对象，在二维目标搜索空间中初始化N个粒子组成一个种群，其中第k个粒子表示为一个二维的位置向量xk＝(Sk,hk),(k＝1,2,…,20)，Sk为搜索窗Ωi的半径。而第k个粒子的移动速度对应的二维的向量为vk＝(vk1,vk2),(k＝1,2,…,20)。
优选的，Sk的搜索范围为[500,2000]，hk的搜索范围为[5,30]。
(3-2)对于信号O，用逐信号点的均方误差的期望E[MSE(Mi)]作为优化衰减参数hi和搜索Ωi的适应度函数去计算每个粒子的适应度Fit(k)。考虑到均方误差的期望E[MSE(Mi)]越小，即表示对应粒子的搜索窗大小Ωi和衰减参数hi降噪效果越好。故适应度Fit(k)采用下式(3)计算：
Fit(k)＝-E[MSE(Mi)] (3)
其中均方误差的期望E[MSE(Mi)]如下式(4)计算：
E[MSE(Mi)]＝bias2(Mi)+var(Mi) (4)
其中偏差平方项bias2(Mi)和方差项var(Mi)计算公式分别如式(5)和式(6)所示：
其中Mi是预滤波信号M中信号点i对应的幅度值,σ′是步骤1所得预估的信号噪声标准差。
(3-3)用下式表示每个粒子迄今为止搜索到的局部最优位置pbest，并利用(3)式计算局部最优位置的适应度即个体最优fpbest(k)：
pbest＝(Sk,hk),(k＝1,2,…,20) (7)
用对应每个粒子的适应度Fit(k)和个体极值fpbest(k)比较，如果Fit(k)大于fpbest(k)，则用Fit(k)替换掉fpbest(k)；
(3-4)用下式表示整个粒子群迄今为止搜索到的全局最优位置gbest，并利用(3)式计算其适应度即全局极值fgbest：
gbest＝(Sg,hg) (8)
对每个粒子，用它的适应度值Fit(k)和全局极值fgbest比较，如果Fit(k)大于fgbest，则用Fit(k)替换掉fgbest；
(3-5)根据下式(9)和(10)更新粒子的速度和位置；
vk＝T*vk+c1r1(pbest-xk)+c2r2(gbest-xk) (9)
xk＝xk+vk (10)
其中c1和c2是学习因子，优选为1.4962，r1和r2是介于[0,1]之间的随机数，T是惯性权重，优选为0.7298。
如果满足最大循环次数退出，否则返回(3-2)。其最大循环次数设置为50。
(3-6)粒子群优化完之后，会得到群体最优的位置gbest＝(Sg,hg)，于是就得到了逐信号点最优化的搜索窗S＝Sg，最优化的衰减参数h＝hg。
步骤4用步骤3粒子群算法所得逐信号最优的搜索窗半径Sg和衰减参数hg对原始信号O进行逐信号点的非局部均值降噪，得到降噪之后的信号M′。该降噪信号中每一信号点i对应的幅度Mi′按式(11)计算：
其中Ωi′表示信号中以信号点i为中心的搜索窗，其窗口大小为步骤3得到的对应信号点i最优的Sg。ω′(i,j)是衡量预处理信号M中分别以i和j为中心的两信号块O′(Ni)和O′(Nj)之间的相似度，其计算公式如下(12)：
O′(Ni)和O′(Nj)分别代表预处理信号M中以信号点i和j为中心的信号块，其信号块大小为10，hg是对应信号点i最优的衰减参数。
为了说明本专利提出的ANLM方法的优越性，我们将传统NLM方法、FIR滤波、IIR滤波、中值滤波、小波变换滤波这五种算法作为对比算法。从图4对比可以明显看出：ANLM滤波后的实际心电图更平滑，在边缘区保留了最多的细节信息。上述对比证实了本专利提出的算法在实际心电图降噪中的有效性和优越性。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种基于粒子群优化的心电信号自适应非局部均值降噪方法，其特征在于，包括以下步骤：
步骤1：对原始心电信号进行非局部均值滤波，得到滤波结果，即预滤波信号；
步骤2：利用粒子群算法来对所述预滤波信号的每一信号点的衰减参数和搜索窗口进行优化,将所述预滤波信号的偏差的平方、方差的和来构建适应度函数，针对所述预滤波信号中的每个信号点，利用该适应度函数对每个点进行处理，进而得到每一个信号点对应的最优化的衰减参数和搜索窗口；
具体包括如下子步骤：
(2-1)初始化粒子群
以搜索窗口Ωi和衰减参数hi为优化对象，在二维目标搜索空间中初始化N个粒子组成一个种群，其中第k个粒子表示为一个二维的位置向量xk＝(Sk,hk),其中k＝1,2,L,N，Sk为搜索窗Ωi的半径， 而第k个粒子的移动速度对应的二维的向量为vk＝(vk1,vk2)，其中k＝1,2,L,N；
(2-2)对于所述原始心电信号O，用逐信号点的均方误差的期望E[MSE(Mi)]作为优化衰减参数hi和搜索窗Ωi的适应度函数去计算每个粒子的适应度Fit(k)，考虑到均方误差的期望E[MSE(Mi)]越小，即表示对应粒子的搜索窗大小Ωi和衰减参数hi降噪效果越好， 故适应度Fit(k)采用下式(3)计算：
Fit(k)＝-E[MSE(Mi)] (3)
其中均方误差的期望E[MSE(Mi)]如下式(4)计算：
E[MSE(Mi)]＝bias2(Mi)+var(Mi) (4)
其中偏差平方项bias2(Mi)和方差项var(Mi)计算公式分别如式(5)和式(6)所示：
其中Mi是预滤波信号M中信号点i对应的幅度值,σ为信号O噪声标准差；
(2-3)用下式表示每个粒子迄今为止搜索到的局部最优位置pbest，并利用(3)式计算局部最优位置的适应度即个体最优fpbest(k)：
pbest＝(Sk,hk),k＝1,2,L,N(7)
用对应每个粒子的适应度Fit(k)和个体极值fpbest(k)比较，如果Fit(k)大于fpbest(k)，则用Fit(k)替换掉fpbest(k)；
(2-4)用下式表示整个粒子群迄今为止搜索到的全局最优位置gbest，并利用(3)式计算其适应度即全局极值fgbest：
gbest＝(Sg,hg)(8)
对每个粒子，用它的适应度值Fit(k)和全局极值fgbest比较，如果Fit(k)大于fgbest，则用Fit(k)替换掉fgbest；
(2-5)根据下式(9)和(10)更新粒子的速度和位置；
vk＝T*vk+c1r1(pbest-xk)+c2r2(gbest-xk) (9)
xk＝xk+vk (10)
其中c1和c2是学习因子，r1和r2是介于[0,1]之间的随机数，T是惯性权重；
如果满足结束条件，即误差足够好或到达最大循环次数，退出，否则返回(2-2)；
(2-6)粒子群优化完之后，会得到群体最优的位置gbest＝(Sg,hg)，于是就得到了逐信号点最优化的搜索窗S＝Sg，最优化的衰减参数h＝hg；
步骤3：用步骤2得到的对于每个信号点最优的搜索窗半径和衰减参数对所述原始心电信号的所有信号点进行逐点的非局部均值降噪，得到降噪之后的信号，由此完成心电信号自适应非局部均值降噪处理。
2.如权利要求1所述的降噪方法，其特征在于：在所述步骤1中，对所述原始心电信号O进行非局部均值滤波，得到预滤波信号M，该信号中每个信号点Mi按式(1)计算
其中，Oj代表所述原始心电信号O中位于点j处信号的幅值，Ω表示搜索窗，ω(i,j)是衡量信号O中分别以i和j为中心的两信号块O(Ni)和O(Nj)之间的相似度，其计算方法如(2)式所示：
其中，O(Ni)和O(Nj)分别代表信号O中分别以i和j为中心的两信号块，Ci为归一化参数，表示L2范数，G表示高斯内核，*表示卷积操作符号，h代表固定的衰减参数。
3.如权利要求1所述的降噪方法，其特征在于：在步骤3中，所述降噪之后的信号中每一信号点i对应的幅度Mi′按式(11)计算：
其中Ωi′表示信号中以信号点i为中心的搜索窗，其窗口大小为步骤3得到的对应信号点i最优的Sg， ω′(i,j)是衡量预处理信号M中分别以i和j为中心的两信号块O′(Ni)和O′(Nj)之间的相似度，其计算公式如下(12)：
O'(Ni)和O'(Nj)分别代表预处理信号M中以信号点i和j为中心的信号块，hg是对应信号点i最优的衰减参数。
CN201510348834.6A 2015-06-23 2015-06-23 一种基于粒子群优化的心电信号自适应非局部均值降噪方法 Active CN105046654B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510348834.6A CN105046654B (zh) 2015-06-23 2015-06-23 一种基于粒子群优化的心电信号自适应非局部均值降噪方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510348834.6A CN105046654B (zh) 2015-06-23 2015-06-23 一种基于粒子群优化的心电信号自适应非局部均值降噪方法

Publications (2)

Publication Number Publication Date
CN105046654A CN105046654A (zh) 2015-11-11
CN105046654B true CN105046654B (zh) 2017-11-10

Family

ID=54453174

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510348834.6A Active CN105046654B (zh) 2015-06-23 2015-06-23 一种基于粒子群优化的心电信号自适应非局部均值降噪方法

Country Status (1)

Country Link
CN (1) CN105046654B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106067783B (zh) * 2016-06-13 2018-11-09 电子科技大学 基于粒子群算法的fir滤波器设计方法
CN109785246B (zh) * 2018-12-11 2021-04-30 奥比中光科技集团股份有限公司 一种非局部均值滤波的降噪方法、装置及设备
CN110507320B (zh) * 2019-09-04 2022-09-23 杭州回车电子科技有限公司 一种脑电信号滤波方法及设备

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101882304A (zh) * 2010-06-24 2010-11-10 西北工业大学 一种sar图像自适应去噪和特征增强方法
CN102609904A (zh) * 2012-01-11 2012-07-25 云南电力试验研究院（集团）有限公司电力研究院 双变量非局部平均滤波x射线图像消噪方法
CN102903084A (zh) * 2012-09-25 2013-01-30 哈尔滨工程大学 一种α稳定模型下的小波域图像噪声方差估计方法
CN102915735A (zh) * 2012-09-21 2013-02-06 南京邮电大学 一种基于压缩感知的含噪语音信号重构方法及装置
CN103576060A (zh) * 2013-10-11 2014-02-12 华南理工大学 基于小波自适应阈值的局部放电信号去噪方法
CN103945091A (zh) * 2014-04-22 2014-07-23 苏州大学 一种基于fpga进化学习的数字图像滤波电路设计方法
CN104706350A (zh) * 2014-11-13 2015-06-17 华中科技大学 一种心电信号降噪方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7856150B2 (en) * 2007-04-10 2010-12-21 Arcsoft, Inc. Denoise method on image pyramid

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101882304A (zh) * 2010-06-24 2010-11-10 西北工业大学 一种sar图像自适应去噪和特征增强方法
CN102609904A (zh) * 2012-01-11 2012-07-25 云南电力试验研究院（集团）有限公司电力研究院 双变量非局部平均滤波x射线图像消噪方法
CN102915735A (zh) * 2012-09-21 2013-02-06 南京邮电大学 一种基于压缩感知的含噪语音信号重构方法及装置
CN102903084A (zh) * 2012-09-25 2013-01-30 哈尔滨工程大学 一种α稳定模型下的小波域图像噪声方差估计方法
CN103576060A (zh) * 2013-10-11 2014-02-12 华南理工大学 基于小波自适应阈值的局部放电信号去噪方法
CN103945091A (zh) * 2014-04-22 2014-07-23 苏州大学 一种基于fpga进化学习的数字图像滤波电路设计方法
CN104706350A (zh) * 2014-11-13 2015-06-17 华中科技大学 一种心电信号降噪方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Region-based non-local means algorithm for noise removal;Zeng W L等;《Electronics Letters》;20111231;第47卷(第20期);第1125-1127页 *
采用结构自适应块匹配的非局部均值去噪算法;钟莹等;《电子与信息学报》;20131231;第35卷(第12期);第2908-2915页 *

Also Published As

Publication number Publication date
CN105046654A (zh) 2015-11-11

Similar Documents

Publication Publication Date Title
Guo et al. A review of wavelet analysis and its applications: Challenges and opportunities
CN105125204B (zh) 一种基于esmd方法的心电信号降噪方法
CN107495959A (zh) 一种基于一维卷积神经网络的心电信号分类方法
CN110680308B (zh) 基于改进emd与阈值法融合的心电信号去噪方法
CN105046654B (zh) 一种基于粒子群优化的心电信号自适应非局部均值降噪方法
CN104706350B (zh) 一种心电信号降噪方法
CN110163825B (zh) 一种人类胚胎心脏超声图像去噪和增强方法
CN109816599A (zh) 一种基于小波分解卷积神经网络的图像条带噪声抑制方法
CN114693561A (zh) 一种基于卷积神经网络的核磁共振图像处理方法与系统
WO2024125321A1 (zh) 一种基于分割熵的油液磨粒特征信号提取方法
Malghan et al. A review on ECG filtering techniques for rhythm analysis
CN107341769A (zh) 一种心电信号除噪方法及系统
Qureshi et al. Analysis of ECG signal processing and filtering algorithms
Yadav et al. Denoising and SNR improvement of ECG signals using wavelet based techniques
Hong et al. FFA-DMRI: A network based on feature fusion and attention mechanism for brain MRI denoising
Rakshit et al. Hybrid approach for ECG signal enhancement using dictionary learning‐based sparse representation
CN117173038A (zh) 一种核磁共振图像增强方法与系统
Wang et al. Research on denoising algorithm for ECG signals
Talbi A new ECG denoising technique based on LWT and TVM
CN109658357A (zh) 一种面向遥感卫星图像的去噪方法
KR101048763B1 (ko) 신호 검출 장치 및 방법
Radhika et al. An adaptive optimum weighted mean filter and bilateral filter for noise removal in cardiac MRI images
CN109242797A (zh) 基于均质和异质区域融合的图像去噪方法、系统及介质
Barmase et al. Wavelet transform-based analysis of QRS complex in ECG signals
CN110263760B (zh) 一种基于病理自适应的心音全变分滤波方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant