基于联合稀疏表示和同时稀疏近似的并行坐标下降去噪算法

何选森 徐 丽

1(广州商学院信息技术与工程学院 广东 广州 511363)2(湖南大学信息科学与工程学院 湖南 长沙 410082)

0 引 言

信号的稀疏表示(Sparse representation)[1]是将信号分解为超完备字典(Over-complete dictionary)[2]少量原子的线性组合。由于观测信号存在噪声,信号处理的任务之一是去噪(Denoising)。而稀疏表示具有潜在的降噪能力[3],联合稀疏表示(Joint Sparse Representation,JSR)[4]和同时稀疏近似(Simultaneous Sparse Approximation,SSA)[5]已被广泛地应用于音频和图像去噪。收缩(Shrinkage)则是基于稀疏表示的著名降(去)噪技术[6],如经典的小波收缩是利用正交变换实现降噪[7],而新的趋势是采用超完备变换替代正交变换。另外,基本追踪(Basis Pursuit,BP)[8]是通过优化将信号分解为字典原子的最佳叠加,从而使表示系数的l1范数(或l0范数)最小化[9],而基本追踪去噪(BP Denoising,BPDN)[10]对于消除高斯白噪声是最佳的。随着超完备字典的发展,迭代收缩(Iterative Shrinkage,IS)算法[11]已成为一种高效的数值技术。Elad[12]提出了四种IS算法以解决BPDN问题,其中的并行坐标下降(Parallel Coordinate Descent,PCD)算法[12]具有最快的搜索速度。文献[13]对PCD进行了收敛性分析,并引入了一种正则化函数形式,使得PCD能对坐标优化提供解析的解。

对于音频信号,用PCD分别处理每一帧能够获得很好的去噪性能。然而,当信号很长时,分割的音频帧数量会很大,用PCD对每帧进行降噪使得计算量急剧增加,造成极大的运行时间成本。为此,本文通过构建以信号(矩阵)为操作对象的时域处理框架,提出一种启发式矩阵型PCD算法(Joint-PCD)。在新的时域框架中,把信号中的每个语音帧看作一个列向量生成信号矩阵,Joint-PCD算法是对一个信号(矩阵)实施降噪,而不仅仅只对一帧处理,因而可以降低算法的运行时间成本。

1 稀疏信号的降噪

若在同一时刻,信号向量x∈RN中只有很少分量是非零的,其他绝大多数分量都为零值,则称x服从稀疏分布,即x可由几个单位能量的原子dk的线性组合来近似[14]。由所有原子构成的集合{dk,k=1,2，…,K}称为字典D∈RN×K。设向量z∈RK是稀疏的,则信号x可以表示为：

x=Dz

(1)

式中:z是由x的K个重构系数组成的列向量。当字典D是正交的,K=N,则表示是唯一的。当字典D是冗余(超完备)的,K>N,则重构每个信号,就有无限多个可能的系数集。正是这种非唯一性提供了表示方法的可选择性。一种更好的表示法是选择具有最少的非零系数的解[15]。

s.t.x=Dz

式(2)是寻找欠定(K>N)线性方程组x=Dz的最稀疏解[15]。类似地,式(2)也可以表示为[15]：

s.t.x=Dz

式中:l1范数是求向量各元素绝对值之和。由于l1范数是l0范数的凸化[15],因此,式(4)也可看作是式(2)的凸化。

对于含噪声的观测信号向量:

y=x+v

(5)

式中:x∈RN是不含噪声的干净信号；v∈RN是有限能量(表示为δ)的高斯零均值白噪声向量。对含噪声信号y进行去噪处理的优化问题可表示为[16]:

式(6)还可以表示为拉格朗日乘数形式[9]:

式中:λ是正的参数,它用于平衡近似误差与系数向量稀疏性二者之间的关系。本质上,式(7)就是基于含噪声信号向量y和超完备字典D的PCD算法的目标函数[12]。

从y中去除高斯白噪声是对以下函数最小化[12]:

式中:等号右边的第一项为干净信号x与观测信号y之间的对数似然；第二项的ρ(x)表示在x上的先验值。f(x)的展开式与ρ(x)有关。将式(1)代入式(8)中得到[12]:

式中:1K表示元素为1的K×1列向量。

IS就是解决式(9)所对应的混合l2-lp(p≤1)优化问题的算法。比较式(9)和式(7)可知,在PCD算法中,先验函数ρ(z)是l1范数,而且PCD的操作对象是向量。

2 新的时域处理框架

用超完备字典原子的线性组合来表示信号向量,就称为简单稀疏近似(Simple sparse approximation),在PCD算法中使用的就是简单稀疏近似。对于很长的音频信号,其分割的帧数就很多,利用PCD分别处理每一帧将造成算法运行时间成本的剧增。为此,需要将以向量为操作对象的模式拓展成以信号(矩阵)为操作对象的模式。

对多个向量同时稀疏近似理论的分析和研究[17],产生出求解SSA问题的策略[18]。在通信系统中,对于被同一信道噪声污染的音频信号,其分割的(帧)向量也都具有相同的分布特性。对不同的帧向量,可以使用同一个字典的不同原子的线性组合来近似计算,即以联合稀疏的模式促进它们非零分量的耦合。因此,为了同时(或并行)地对一个信号的所有帧(向量)进行降噪操作,需要将传统的帧处理模式转变成如下新的时域处理框架。

根据分割规则,将音频信号分割成一系列的帧,每个帧作为一个列向量生成矩阵,使得信号与矩阵唯一的对应。由于一个矩阵中所有的列向量对应于同一个音频信号的不同帧,因此这些列向量具有相同的稀疏性。另外,由于每个列向量都受到同一个信道噪声的污染,因而不同噪声分量之间也不需要统计独立。基于这样的信号矩阵,可以同时处理多个列向量以节省降噪的时间开销。

3 Joint-PCD启发式算法

在新的时域处理框架之下,将简单稀疏近似的概念推广到同时稀疏近似(SSA),并以此为基础,将PCD扩展为矩阵型的Joint-PCD算法,其基本原理如下。

令y1=x1+v1、y2=x2+v2分别为干净信号x1,x2∈RN的含噪声观测量,其中v1,v2∈RN为加性噪声。假设信号x1、x2可以由超完备字典D∈RN×K近似表示为x1=Dz1、x2=Dz2,其中z1,z2∈RK为x1、x2的重建系数。由式(7)可得:

联合考虑以上这两个优化问题,则有：

(11)

再令X=[x1,x2]、Y=[y1,y2]、Z=[z1,z2],则将式(11)简化为如下标准的稀疏编码问题:

对Z和D二者来说,式(12)不是一个凸优化,但当固定其中一个时,式(12)就成为凸优化问题[4]。本文中,字典D是给定的,而通过优化系数矩阵Z来构造Joint-PCD算法。

式(12)是对两个信号y1、y2同时去噪的优化问题。显然,优化式(12)可以推广到多个信号向量的情况。考虑矩阵形式的观测信号模型为[19]：

Y=X+V

(13)

式中:X∈RN×P是包含P个向量xi∈RN的干净信号矩阵；V∈RN×P是P个向量vi∈RN组成的噪声矩阵；Y∈RN×P是含噪声观测信号矩阵。基于联合稀疏表示(JSR)[4],信号X可以由字典D表示为X=DZ,于是观察信号矩阵Y为[19]：

Y=DZ+V

(14)

式中:ρ(·)是任意标量函数,而ρ(Z)是矩阵Z稀疏性的惩罚函数。式(16)就是Joint-PCD算法的目标函数。

由式(15)和式(16)可归纳出构造Joint-PCD算法的基本思想。为了限制对观测信号Y进行逼近的字典原子的总数,需要使Z中非零行向量的数量最小化。为此,从两方面来考虑。首先,参与对Y近似的字典原子总数越少越好;其次,每个参与近似的原子要对尽可能多的Y的列向量做出贡献。也就是说,系数矩阵Z是稀疏的,即它的大多数行向量都是零向量,而Z的每个非零行向量都应该有尽可能多的非零元素。基于此,将稀疏性惩罚函数ρ(Z)设置为行l0准范数(Row-l0quasi-norm)的松弛版:

把式(17)代入到式(16)中,则观测信号矩阵Y的降噪任务变成为混合范数的优化问题:

同样地,由式(12)可知,惩罚函数ρ(Z)是矩阵Z的l1范数。在本文中,对ρ(Z)的l1范数使用偶对称函数:

式中:zij是矩阵Z的第(i,j)个元素。

这里需要提示,在式(12)中Z和Y都是由两个向量构成的矩阵,然而,在下文中,将两个向量情况推广到多个向量的情况,即式(12)中的Z和Y是由P(P≤2)个列向量组成的矩阵Z∈RK×P和Y∈RN×P,而不仅仅是由2个列向量组成的矩阵Z∈RK×2和Y∈RN×2。

U(D,Y,Zk)=diag-1{DTD}DT(Y-DZk)+Zk

(20)

一般地,超完备字典D的列是归一化的,即diag-1(

DTD)=I,其中I是单位矩阵,则可将式(20)简化为:

U(D,Y,Zk)=DT(Y-DZk)+Zk

(21)

矩阵Z的更新迭代值Zk+1为:

Zk+1=Zk+μQk

(22)

式中：μ是迭代步长参数。

Qk=S[U(D,Y,Zk),λ]-Zk

(23)

式中：λ是收缩函数S(·)的阈值。从四种IS算法的推导过程[12]可知,式(23)中的收缩函数S(·)是一个软阈值函数,其定义为：

式中:θ为软收缩函数S(·)的阈值。

另外,Joint-PCD算法的性能会受到式(22)中步长参数μ影响。对足够小的μ>0,步长要保证惩罚函数ρ(·)的可行下降。本文使用线搜索[11]来寻找最佳的步长μ值。线搜索作为一维优化,它用于解决以下的问题:

式中:ρ(·)可以是矩阵的rx范数,也可以是矩阵的l1范数。

若信号矩阵Y退化成一个列向量y,系数矩阵Z也退化成一个列向量z,则Joint-PCD算法退化为PCD算法。因此,将Joint-PCD称为启发式算法。

在音频信号处理中,经常使用离散余弦变换(DCT)字典和Gabor字典。一个窗口式DCT字典为DDCT=[d1,d2,…,dK],它的第k个原子定义为[20]：

一个窗口式Gabor字典为DG=[d(k,φ)](k,φ)∈Ψ,其原子以连续的参数集Ψ=[1,K]×[0,2π]为索引,其定义为[20]：

式中:1≤k≤K,0≤t≤N-1,N是音频帧的长度,K是字典的尺寸；wd是加权窗口；φ是相位。Gabor字典连续索引原子d(k,φ)的分解可以用离散字典中的原子对来表示[20]。

式中：原子对是相同频率且相位为零的余弦和正弦对。为简化计算,词典的列都做归一化处理,如Gabor字典的正弦和余弦原子都使用对应的单位范数版本。

综上所述,对启发式Joint-PCD算法描述见算法1。

输入:观测矩阵Y，字典D，阈值λ，函数ρ(·)，迭代次数M。

Begin:设置k=1,迭代开始。

ifk

误差:计算E=Y-DZk;

映射:计算ET=DTE;

收缩:计算ES=S(ET+Zk,λ);

线收缩:寻找μ0以获得Zk+μ(ES-Zk)的下降;

放松:更新迭代Zk+1=Zk+μ0(ES-Zk);

Next:k←k+1;

End

4 仿真实验及结果分析

4.1 仿真环境和性能指标

为了验证Joint-PCD算法降噪的有效性,对Joint-PCD和PCD在音频信号降噪性能上进行对比分析。不同算法的仿真环境是完全一样的,以便得到公平的结论。

仿真平台为笔记本电脑,其CPU为Intel(R) Celeron(R) 1007U CPU-1.50 GHz,内存4 GB,操作系统为Windows 10。仿真软件平台为MATLAB 9(R2016a)。

降噪的性能用信噪比(signal to noise ratio,SNR)和均方根误差(root mean square error,RMSE)来衡量:

测试用的五个音频音乐信源[21]的类型为.wav,信号采样率为44.1 MHz。每个信源的样本数为L=216=65 536。五个信源均服从超高斯分布,它们之间的区别体现在具有不同的峭(峰)度(kurtosis)值。信源X=[x1,x2,x3,x4,x5]T的归一化峰度值分别为2.400 8、4.776 6、1.052 8、4.690 2、3.133 7,其对应的时域波形如图1所示。

图1 五路音乐音频源信号的波形

降噪的仿真过程如下。音频信源与高斯噪声相加生成含噪声信号。首先,每个含噪声时域信号乘以长度为64的汉明窗生成语音帧,连续两帧之间的重叠样本数为32,即每帧的前后各有16个样本与相邻帧的样本重叠。分帧过程是通过使用Voicebox[22]中的MATLAB函数enframe来实现。然后,将PCD和Joint-PCD算法分别应用于分割后的信号以进行降噪。每次使用PCD是对音频信号的一帧进行处理;而使用Joint-PCD算法每次是对一个信号(而不仅是一帧)进行处理。最后,重建降噪后的信号,即每个降噪的帧乘以逆汉明窗,并且在相邻两帧重叠的样本中,分别在两边去除50%的样本,仅保留每一帧的中间部分的样本,将经过处理的各个片段简单串联起来即可得到重新合成的音频信号。

4.2 算法参数设置

在PCD算法中,需要给定软收缩函数S(·)的阈值λ。常用的阈值规则为全局阈值(universal threshold)[7]和最小最大阈值(minimax threshold)[7]。本文采用全局阈值的设置,满足以下条件:

式中:L是信号的样本数量；σ是高斯噪声的标准差。显然,要确定阈值λ,则必须先确定噪声的强度σ。

为了确定噪声强度对降噪性能的影响,采用经典的小波降噪方法进行仿真。对5个信源X加噪声形成观测信号Y,当Y的信噪比分别为10、15、20时,利用小波重复进行100次降噪实验,并计算在每次实验中5个信号的平均SNR值和RMSE值,其结果如图2所示。

图2 5个信号的平均SNR和RMSE值

从图2可知,随着输入SNR的降低(噪声强度的增大),小波的降噪性能变差。当输入SNR为10时,经过降噪后5个信号的平均SNR小于9.8,即降噪后的信噪比小于输入的信噪比。对音频信号来说,这是一个非常恶劣的环境。

为了便于分析,对不同的输入SNR,计算出每个信源被施加高斯噪声的强度σ。其结果如表1所示。

表1 输入SNR与噪声强度σ的对应关系

从表1中的数据可知,由于每个信源的幅度(功率)不同,对同一个输入SNR值,其对应于不同信源的噪声强度也是不同的。对于输入SNR为10的情况,信号x2、x3、x4、x5上施加的噪声强度接近0.1,但信号x1上施加的噪声强度为0.2。为获得一个统一的噪声强度,在下面的仿真中,5个信源都被施加σ=0.1的高斯噪声形成观测信号。

给定噪声强度后,再来确定软阈值λ。由于信号采样点数为L,则(2log2L)1/2=4.709 6,由文献[7]可知,在给定σ=1的情况下,其参考的阈值为λ*=3.31。将真实的σ=0.1代入可得λ≤σ(3.31=0.331,这也只是λ的取值范围,还需要通过仿真来寻找PCD算法最佳的阈值λ。仿真条件为：每个含噪信号都被分割成长度为64的帧,PCD的迭代次数设置为M=5,阈值λ从0.050到0.331变化,更新步长为0.01,每个λ值做一次实验,超完备字典是DCT和Gabor。各个信号的性能指标SNR和RMSE如图3所示。

(a) DCT

(b) Gabor图3 每个信号对不同λ的SNR和RMSE值

从图3可知,无论是DCT还是Gabor词典,当λ=0.1时,每个信号的SNR达到最大值,而RMSE也达到对应的最小值。因此在PCD算法中,取λ=0.1作为最佳的阈值。另外,从图3还可看出,采用不同的超完备字典,PCD算法的降噪性能是有差异的。为了比较DCT和Gabor的性能,将迭代次数设置为M=15,并分别计算PCD对每个信号在每次迭代中的SNR和RMSE值。其结果如图4所示。

图4 每个信号对不同字典的SNR和RMSE值

为便于分析,仅观察信号1(signal1)的情况。从图4可知,使用Gabor字典的SNR和RMSE值都优于DCT字典对应的指标。其他信号的SNR和RMSE值也有相似的关系。为此,在PCD和Joint-PCD算法中,使用Gabor字典。

另外,在算法的参数中,还需要确定惩罚函数ρ(·)及迭代步长μ的值。在下面的仿真中,对于PCD算法,采用向量的l1范数求解式(7);对于Joint-PCD算法,采用矩阵的l1范数求解式(12),采用矩阵的rx范数求解式(18)。为了提高算法的性能,各算法的步长μ(0<μ<1)采用MATLAB函数fminbnd通过在线搜索最优的μ值作为其迭代步长。这样就完成了PCD和Joint-PCD算法中所有参数的设置。

4.3 算法性能比较

在对算法去噪性能的比较过程中,为方便起见,把rx范数和l1范数的Joint-PCD分别表示为Joint-PCD withrxnorm和Joint-PCD withl1norm。

关于算法的性能,从两个方面对PCD和Joint-PCD进行比较。(1) 在相同条件下算法所消耗的运行时间成本;(2) 去噪的性能指标值SNR和RMSE。在这个仿真实验中,将算法的最大迭代次数M设置为50以进行全面的考察。其运行时间如图5所示。

图5 算法的运行时间比较

可以看出,算法的运行时间随迭代次数的增加基本上呈线性增长。显然,rx范数和l1范数的Joint-PCD算法的运行时间增长率几乎是相同的;而PCD算法的运行时间增长率却远远大于Joint-PCD的增长率。

图6分别给出了PCD算法、Joint-PCD withl1norm和Joint-PCD withrxnorm算法在该仿真中去噪性能指标值SNR和RMSE的波形图。

(a) PCD

(b) Joint-PCD with l1 norm

(c) Joint-PCD with rx norm图6 PCD和Joint-PCD算法对每个信号的SNR和RMSE值

可以看出,PCD与Joint-PCD算法的降噪性能几乎是完全相同的,而且经过大约25次迭代,三种算法的降噪指标SNR和RMSE值都逐渐收敛到其稳定值。

为了比较算法的平均降噪性能,将每种算法对5个信号的平均SNR和平均RMSE值绘制在图7中。

(a)

(b)图7 5个信号的平均SNR和RMSE指标

可以看出,PCD和Joint-PCD算法的平均SNR和RMSE值在算法的初始阶段随着迭代次数的增加具有一定的波动,这是因为每个信号的SNR和RMSE值具有不同的波动范围。对于PCD和Joint-PCD withrxnorm算法,平均降噪指标大约在经过20次迭代后结束波动,趋于平稳;而Joint-PCD withl1norm算法大约在经25次迭代后其降噪指标趋于平稳。也就是说,PCD和Joint-PCD的平均降噪性都能在25次迭代后收敛到稳定的指标值。

为了从数值上直观地比较PCD和Joint-PCD在降噪性能上的差异,在以下仿真中,将算法的最大迭代次数M设置为25。表2给出PCD和Joint-PCD的运行时间。

表2 PCD和Joint-PCD算法的运行时间 单位:s

可以看出,l1范数的Joint-PCD运行时间是最少的,而rx范数的Joint-PCD运行时间也仅仅是多了0.379 2 s,因此,两种Joint-PCD算法的运行时间几乎是相同的。相反地,PCD算法的运行时间大约是Joint-PCD的4.6倍,这说明在音频信号的降噪过程中,Joint-PCD算法降低了计算负担。显然,与PCD算法相比较,Joint-PCD算法的收敛速度得到了有效的提高。

对于PCD和Joint-PCD算法的降噪性能,采用每个信号在整个迭代过程中SNR和RMSE的平均值来度量。表3给出了每个信号的平均降噪性能指标。

表3 每个信号在整个迭代过程的平均SNR和RMSE值

可以看出,Joint-PCD withl1norm算法对于第2、第4和第5个信号的平均SNR和RMSE值是最优的;而PCD算法仅对于第1和第3个信号的平均指标值是最优的。另外,PCD算法和Joint-PCD算法的平均SNR值的差别仅体现在小数点后面第二位;而算法的平均RMSE值的差别仅发生在小数点后面第四位上。因此,从去噪性能指标的差异性来看,利用Joint-PCD算法代替PCD算法对音频信号进行降噪处理,其降噪的效果几乎完全相同。

5 结 语

为了解决PCD算法在实际应用中的计算成本问题,本文定义一种新的音频信号时域处理框架,将每个被分割的音频帧作为一个列向量生成信号矩阵。在联合稀疏表示(JSR)和同时稀疏近似(SSA)的基础上,利用新的处理框架,提出以矩阵为操作对象的Joint-PCD算法。仿真的结果表明,Joint-PCD算法不仅能够提供与PCD算法几乎完全相同的降噪性能,而且Joint-PCD算法提供了更高的计算效率,节约了算法的运行时间成本，从而加快了算法的收敛速度。更为重要的是,与PCD算法相比较,Joint-PCD算法在大规模数据处理和实时信号处理应用方面具有巨大的潜在优势。