基于全变分模型稀疏信号降噪算法的研究

王 娜,吕东澔,张 勇

(内蒙古科技大学信息工程学院,内蒙古 包头 014010)

1 引言

随着科技的进步和社会的发展,数字化的信号已经成为当前时代人们获取信息的最主要来源。由于设备内部电子元器件老化、环境中杂散的电磁波辐射干扰等因素,信号在采集、传输和处理的过程中必然会夹杂一定的噪声,噪声的出现会使原始信号和接收到的信号之间产生误差,导致在后续分析时原始信号出现失真的现象。因此想要从信号中提取有用的信息,首先需要对其进行降噪处理,信号降噪是信号处理的一个核心任务,其目的是抑制含噪信号中噪声成分并增强有用的信号成分,有效的降噪方法往往对更高层任务的成败具有重要影响。相关学者在信号去噪处理方面进行了大量的研究,例如,Gabor[1]提出了短时傅里叶变换算法,信号经过傅里叶变换得到频率信号,可在频域分析信号的频谱,但是不能提取时变、非平稳信号的局部特征;模极大值去噪[2]、小波阈值去噪[3]是基于小波思想提出的典型去噪算法,主要是通过设置小波系数来达到抑制噪声的目的,然而在小波域中执行软阈值处理时会产生伪像现象。现有的降噪算法虽具有良好的降噪性能。但是利用这些算法在对稀疏信号进行降噪时会产生一些阶梯伪像噪声,如寄生噪声尖峰和Gibbs振荡[4,5],导致原始信号中重要信号丢失。而Rudin Osher和Fatemi[6]基于变分法的思想提出基于全变分(Total Variation, TV)正则化的降噪方法,由于该方法主要是利用了信号在时域内的光滑性,在对稀疏信号进行降噪时能保留信号的边缘信息并准确还原原始信号。

TV去噪是一种基于基础信号为分段常数信号(包含信号的导数是常数这一情况)的非线性滤波方法[7],这种信号出现在地球科学、生物物理学和其它领域[8],TV去噪技术还与其它方法结合使用,以处理更一般类型的信号[9-11],但是TV去噪方法在处理稀疏信号时更为突出。TV降噪算法是通过对一个特定代价函数进行优化,目前算法研究大多是基于凸惩罚函数提出的代价函数。例如迭代软阈值算法(Iterating Soft Threshold Algorithm, ISTA)[12]用于解决1范数正则化反问题,ISTA具有线性收敛性能够实现对稀疏信号的降噪目的,但是需要经过多次迭代后才能收敛且降噪精度较差。Selesnick I W等人提出的优化最小化(Majorization Minimization, MM)[13]算法、迭代裁剪(Clipping, CLIP)[14]算法等能够有效的找到函数的最小值,达到降噪的目的,但是不能精确估计分段常数信号峰值处的幅度。这些降噪算法的代价函数都是以凸惩罚函数为前提进行的,因为凸惩罚中的1范数能够提高信号的稀疏性。

研究表明非凸罚函数可以更有效的使信号稀疏化,相关学者[15-18]利用非凸惩罚函数的特点来估计稀疏信号,但是此类算法中很少能保持代价函数的凸性。本文基于非凸惩罚项构建了算法的代建函数,它推广了标准的惩罚函数,通过选择合适的优化参数使其满足凸性条件,使去噪问题的代价函数保持凸性。新的惩罚可以更准确的估计分段常数信号中峰值处的幅值,降噪精度高于基于凸惩罚的MM优化算法以及CLIP算法。

2 全变分降噪

TV去噪算法在降噪领域中应用广泛,该方法的具体思路是把需要降噪信号的代价函数构建成一个凸函数,对这个凸函数进行求导,找到凸函数的最小值,这个最小值就是利用全变分降噪得到的去噪信号。

假设噪声数据y(n)为

y(n)=x(n)+w(n),n=0,…,N-1

(1)

其中x(n)是近似分段常数信号,w(n)是高斯白噪声。TV去噪通过解决以下优化问题来估计信号x(n)。

(2)

其中x*为降噪后信号,λ为正则化参数且λ>0用于平衡函数的权重,φn(xn)为正则项且φn(xn)>0,可以通过设计不同的正则化项设计降噪算法可以实现含噪信号的精准降噪。

N点信号x由以下矢量表示

x=[x(0),…,x(N-1)]T

(3)

矩阵D被定义为

(4)

N点信号x的一阶差分矩阵由Dx给出,其中矩阵D为(N-1)×N的一阶差分矩阵。

使用以上表示方法,TV去噪问题(2)可简洁的写成

(5)

代价函数(5)为全变分去噪模型的标准函数,本文算法是在此标准函数的基础上进行的推广设计,结合非凸惩罚项提出的全变分(Non-Convex Penalty Total Variation, NCP-TV)降噪算法。

3 本文算法

3.1 可微分凸函数

本文所提出的非凸惩罚函数是从标准的惩罚函数中减去一个可微的凸函数,其中可微的凸函数是可分离的或二元函数的和,首先定义可微的凸函数如下

令0<α<1,定义Kα:RN→R

(6)

D为式(4)中提到的一阶差分矩阵。

其中Kα为2范数,2范数可以寻求最优问题且为二范数可微易求解。通过定义可知Kα存在二阶导且二阶导数为α,又因为0<α<1,由此可知函数Kα是可微的凸函数。

命题1:.函数Kα可定义如下

K0(x)=0

(7)

(8)

证明:

若α=0:设α=0,b=0代入式(6)可得式(7)。

若0<α<1,通过对函数(6)进行最小化可得函数b的定义即

(9)

=[(Db)′Db]′+α(b-x)

=D′Db+b′D′D+α(b-x)

=2D′Db+αb-αx

(10)

令:s′(b)=0则

2D′Db+αb=αx

(2D′D+α)b=αx

(11)

b=(2D′D+αI)-1αx

3.2 非凸惩罚函数

研究证明非凸惩罚函数可以更有效的提高信号的稀疏性,所以本文改进的算法是基于非凸惩罚的全变分去噪算法,此函数是标准全变分函数的一种推广形式。

命题2:φα(x)设计为非凸函数,其定义如下

(12)

D为一阶差分矩阵。

命题3:当0<α<1,非凸惩罚函数φα(x)的非凸区间定义如下

(13)

证明:

当一个凸函数减去一个凸函数时,得到的函数很可能在某一定义域内为负值,所以定义如(13)所示的非凸区间:

由式(6)定义可知函数Kα(x)为二次函数所以

Kα≥0

(14)

又因为函数φα(x)为正则化函数所以

(15)

可得

(16)

由此可推导出式(13)。

3.3 代价函数

命题4:.由全变分降噪定义可得基于非凸惩罚的代价函数为:

(17)

对此函数进行降噪本质上是求其最优解的问题即

(18)

证明:

令n=1式(18)可写成以下形式

(19)

构成关于x的二项式函数其中c为关于b的常数项,二次函数中二次项系数大于0,函数开口向上即为凸函数由此可得

(20)

得

(21)

3.4 降噪算法

命题5:设y∈RN,λ>0,算法中产生迭代,其迭代式x(k)定义如下所示

B(k)=y+λα[x-(2D′D)+αI]-1αx

(22)

x(k+1)=tvd(zk,λ)

(23)

证明:

NCP-TV降噪算法使用前后向分裂方法(Forward Backward Splitting, FBS)[19]对代价函数进行优化,FBS方法是一种迭代优化算法,该算法最小化的函数形式如下所示

F(x)=f1(x)+f2(x)

(24)

其中f1(x)和f2(x)都是凸函数,FBS方法如下

B(k)=x(k)-∇f1(x(k))

(25)

(26)

通过x(k+1)进行迭代至收敛到F(x)的最小值。

对式(17)使用FBS算法可得

Aα(x)=Fα(x)

=f1(x)+f2(x)

(27)

其中

(28)

(29)

f1(x)的梯度为

∇f1(x)=x-y-∇Kα(x)=x-y-α(x-b)

(30)

将式(11)中b代入(27)可得

∇f1(x)=x-y-αλ[x-(2D′D+αI)-1αx]

(31)

将式(26)(28)代入上式(22)(23)可得(19)(20)形式的FBS算法。

基于全变分模型的降噪算法步骤如下:

1)输入含噪信号y,初始化x,λ和α输入参数;

2)设置k=1:Nit,其中Nit为需要迭代的次数;

3)输入中间变量B(k);

4)对x(k)进行全变分降噪;

5)k=k+1;

返回步骤2),反复迭代降噪,直到收敛到代价函数的最优值,迭代结束,最优降噪次数为Nit;

4 数据分析及实验结果

4.1 仿真信号降噪

本文实验通过对两种不同的信号进行降噪处理,第一种为仿真信号由Matlab生成的分段常数信号,在此信号的基础上人为的加入一组高斯噪声,通过使用MM、CLIP降噪算法和NCP-TV降噪算法对此含噪信号进行降噪处理,对比三种算法的降噪效果;第二种是NCP-TV算法对实际采集的肌电信号进行降噪处理,分析降噪效果并验证NCP-TV降噪算法在实际应用中的有效性。

对如图1所示的含噪分段常数信号进行降噪处理,这是带有加性高斯白噪声(σ=0.5)的信号(长度N=256)。

图1 含噪信号

如图2显示了CLIP、MM算法和NCP-TV算法的去噪结果,其中CLIP算法是基于TV降噪的一种降噪方法,由仿真结果可看出信号平滑区域降噪结果较差,MM算法是标准惩罚函数的一种全变分去噪方法,由图可以看出MM算法降噪后的信号基本可以跟随原始信号,但是峰值的降噪效果有一定偏差,从NCP-TV算法的去噪结果可以看出本文算法可以更精确的估计信号峰值及平滑区域信号的幅值,降噪精度优于上述两种算法。

图2 三种降噪算法对比

为了更加直观的分析三种降噪算法的降噪效果,计算了三种算法降噪后的误差曲线如图3所示,让噪声标准差满足0.2≤σ≤1.0,分别计算每一个σ值对应的噪声方差,方差用来度量随机变量和数学期望(即均值)之间的偏离程度。由图3可知NCP-TV算法降噪后信号的方差明显低于CLIP算法和MM算法,说明NCP-TV算法降噪后信号与原始信号之间的误差较小降噪精度较高。

图3 两种算法误差对比

分析信号降噪常用的分析方法有信噪比(SNR)、均方根误差(RMSE)、平均绝对偏差(MAE)等。

其中SNR为有用信号功率与噪声功率的比值,以dB为单位,SNR越高,代表信号中噪声值越小,故信号精度越高。其公式为

(32)

其中Ps表示有用信号功率,Pn表示噪声的有效功率。

RMSE、MAE可以反映预测模型与实际数据误差值大小,数值越低代表两者越接近精度越高其计算公式如下

(32)

(33)

其中yi表示真实值即原始信号fi表示预测值即降噪后的信号。

三种算法的分析结果如表1所示:由表可知三种降噪算法的SNR值都高于含噪信号,说明三种算法都达到了降噪的目的,其中NCP-TV算法的SNR值较其它两种算法最高说明其信号中噪声值最小降噪效果较好;由RMSE、MAE值对比可知,NCP-TV算法的RMSE和MAE值最低说明其降噪后的数据与原始信号误差值最小,降噪精度较高。

表1 两种不同算法特征指标对比

4.2 肌电信号降噪

在实际采样中,生物组织受到高压电脉冲后期的电流信号均为方波信号,求导信号值为零,通过判断零信号和求导后的零信号占总体信号一半以上,这样确定的信号是稀疏的。如图4所示为实际采集到肌电信号的部分放大数据多数为方波信号可知肌电信号也为稀疏信号。

图4 肌电稀疏信号

图5为实际采集到肌电信号的完整数据,采集过程中由于环境或设备等干扰因素肌电信号中夹杂大量干扰信号,所以在对信号进行分析之前需要对其进行降噪处理,可采用NCP-TV算法验证在实际应用中的有效性。

图5 肌电含噪信号

通过上一章节的介绍确定NCP-TV降噪算法的可行性,为了提高肌电信号的质量采用NCP-TV算法对实际采取的肌电信号进行降噪,降噪实验如图6所示。

图6 肌电信号降噪

为了更加清楚地观察降噪效果,提取信号中的前3500个数据如图7所示为原始肌电信号,通过NCP-TV算法进行降噪,降噪后的信号如图8所示,从降噪效果来看,所采用的降噪算法能够消弱噪声对肌电信号的影响。

图7 肌电含噪信号

图8 NCP-TV算法降噪

为了测试算法的适用性,对信号进行特征提取并分析,肌肉疲劳标准是随着疲劳程度的加深,即平均功率频率(MPF)值的降低,MPF反映信号频谱含量的变化,一些影响信号频谱的因素会影响该参数对疲劳的表征效果,有些疲劳过程中MPF可能下降不明显,其计算公式如下所示

(34)

式中p(f)为功率谱密度,对原始肌电信号段作频谱变换得到,f1和f2是硬件滤波的截止频率。

表2是肌电信号静止状态与疲劳状态降噪前后的MPF值,可以发现利用NCP-TV算法去噪后肌肉由静止到疲劳状态肌电信号MPF值降低的速率增加,从而说明降噪后的信号能更好的反映肌肉疲劳的标准。

表2 疲劳程度MPF值

表3为肌电信号静止状态与疲劳状态降噪前后的MAE值前后对比,随着疲劳程度的加深信号的MAE值呈增长趋势,由表可以看出信号降噪后其MAE值较降噪前增长的速率明显提升更好的反映了肌肉的疲劳标准。

表3 疲劳程度MAE值

由于肌电信号可以反映肌肉的活动特征和肌肉疲劳等信息,但因其包含多种噪声干扰无法直接进行应用,通过NCP-TV降噪算法对其进行降噪处理,降噪后信号才可以进一步进行分析处理应用于疲劳分析等领域。实验结果表明,NCP-TV算法在实际应用中具有有效性。

5 结束语

本文研究了稀疏信号的精确降噪问题,提出了NCP-TV降噪算法,该算法通过构建新的代价函数利用前向-后向方法对稀疏信号进行降噪,并将全变分模型建为最优模型,通过迭代寻优获得该模型的最优解,该算法能够有效的提高信号幅值的准确性,仿真结果表明,NCP-TV降噪算法是有效和可行的,降噪效果优于其它对比算法并且在实际应用中具有有效性。