啟發式聯合PCD快速降噪算法

何選森，徐 麗，許 瑩

(1.廣州商學院信息技術與工程學院 廣州511363；2.湖南大學信息科學與工程學院 長沙 410082)

在對信號的描述中，基于冗余字典的稀疏表示提供了一種強有力的信號模型，它是將獲得的觀測信號近似為字典中少量原子的線性組合[1]。由于噪聲無處不在，純粹的干凈信號是無法得到的，對信號降噪已經成為信號處理最重要的任務之一。稀疏表示本身具有潛在的降噪能力[2]，特別是聯合稀疏表示(joint sparse representation,JSR)[3]和同時稀疏近似(simultaneous sparse approximation,SSA)[4]在信號降噪領域都獲得了廣泛的應用。另外，收縮是建立在稀疏表示基礎上的具有吸引力的降噪技術[5]，著名的小波收縮使用正交變換實現降噪[6]。然而，新的趨勢是采用超完備變換，它對信號的分解不是唯一的。例如：幀方法、匹配追蹤法、最佳正交基法等，而基本追蹤(basis pursuit,BP)[7]法是通過優化技術將信號分解為字典元素的最佳(疊加)表示，且使表示系數的l1范數(或l0范數)最小化[8]。將BP應用于含噪聲的信號，即基本追蹤降噪(BP denoising,BPDN)[9]方法對于消除高斯白噪聲非常有效。同樣地，冗余字典的發展也催生了高效的數值技術，即迭代收縮(iterative shrinkage, IS)算法[10]。為解決BPDN問題，文獻[11]提出了4種IS算法，分別是啟發式收縮、基于迭代重加權最小二乘、階段式正交匹配追蹤、并行坐標下降(parallel coordinate descent,PCD)算法。其中PCD的搜索速度最快。因此，文獻[12]又對PCD進行了算法的收斂性分析，并引入了一種正則化函數形式，使得PCD能對坐標優化提供解析的解。

利用PCD對每個含噪聲的音頻(語音)幀進行降噪能獲得整個信號很好的降噪效果。然而，若信號很長，則分割的音頻幀數量就很大，用PCD對每幀分別進行降噪將導致繁重的計算負擔。為此，本文建立了一種以信號矩陣為操作對象的時域處理框架，并提出了基于SSA和JSR的啟發式聯合PCD算法(Joint-PCD)。在此框架中，每個音頻幀作為一個列向量生成信號矩陣，Joint-PCD算法每次是對一個信號矩陣(而不僅僅是對一個音頻幀)進行降噪，從而極大地降低了算法的運行時間開銷，提高了算法的收斂速度。

1 稀疏信號的降噪

在對信號進行表示時，最常用、最簡單的方法就是對信號進行線性變換。用s∈RN表示N維信源矢量，x∈RN表示對應的觀測矢量，T表示N×N變換矩陣，則x=Ts。由于變換后的信號數量等于信源的數量，則T稱為正交變換。信號的另一種重要表示是線性稀疏表示。假設在同一時刻只有少量的分量xi∈x是非零的，其他絕大多數分量都為零值，則稱x服從稀疏分布。若信號x被高斯噪聲v污染，通過對x的稀疏分量進行軟閾值收縮可實現x的降噪[6]。

信號x∈RN一般由幾個單位能量原子dk的線性組合近似表示，所有原子組成的集合{dk,k=1,2,···,K}構成一個字典D∈RN×K。觀測向量x可以表示為：

式中，z∈RK是稀疏的，是K個系數的列向量。當K=N時，D為正交，其表示也是唯一的；當K>N時，D為冗余，式(1)存在無限多個可能的系數集。冗余字典的優勢體現在可以從眾多的可能性中選擇最適合的一種表示，其中最吸引人的表示法是求具有最少數量的非零的系數解：

式(2)是尋找欠定線性方程組x=Dz的最稀疏解[13]。優化式(2)還可以等價地表示為[13]：

式中，ε是誤差容限；||·||p表示lp范數；l0范數是對一個向量非零元素的計數；l2范數采用歐幾里得范數。式(2)和式(3)可用正交匹配追蹤法求解[11]。式(2)還可表示為[7]：

式中，l1范數是求一個向量元素絕對值的總和。由于l1范數是l0范數的凸化[13]，因此，式(4)可以看作是式(2)的凸化。

對于含噪聲的觀測信號矢量：

式中，x∈RN是干凈的信號向量；v∈RN是能量為δ的高斯零均值白噪聲向量。對信號y降噪的優化問題表示為[8]：

優化式(6)，可由拉格朗日形式代替[8]：

式中，λ是用于調節近似誤差與系數向量稀疏性之間關系的一個參數。實際上，式(7)就是PCD算法的目標函數[11]。

一般地，從含噪聲信號矢量y中去除高斯白噪聲的過程可歸結為對下式函數的最小化[11]：

式中，第一項為對數似然，它描述了干凈信號x與含噪聲向量y之間的關系；第二項的ρ(x)表示x的先驗值，不同的ρ(x)對應于不同的f(x)展開式。將式(1)代入式(8)中，得到信號降噪的目標函數為[11]：

式中，1K表示元素為1的K×1矩陣(列向量)。

迭代收縮(IS)算法是為了解決目標函數為式(9)這樣的混合l2?lp(p≤1)優化問題而產生的新型高效數值技術。比較式(9)與式(7)可知：在PCD算法中，函數ρ(z)是l1范數。

式(7)表示的PCD算法的操作對象是向量，在音頻信號降噪過程中，PCD每次是對一個音頻幀(列向量)降噪。若音頻幀數量巨大，則PCD算法的運行時間將急劇增加。特別對于大規模數據和實時信號的降噪問題，分別用PCD對每幀降噪是不現實的，必須尋找新的方法同時(并行)對多個列向量降噪，以減少算法的計算負擔。

2 Joint-PCD啟發式算法

用冗余字典原子的線性組合近似表示一個音頻幀，稱之為簡單稀疏近似(simple sparse approximation)；對于被相同信道噪聲污染的多個幀(矢量)，用相同字典原子的不同線性組合同時近似表示，稱之為JSR。通過對多個向量同時稀疏表示的理論分析和研究[14]，從而就產生了求解SSA問題的策略[15]。對于音頻(語音)信號在加性高斯白噪聲(additive white Gaussian noise, AWGN)信道的傳輸與降噪，SSA具有重要意義。

為了對一個音頻信號的所有幀同時進行降噪，需建立新的時域框架以區別傳統的幀處理模式。這種框架的基本思想為：把信號中每個被分割的幀按順序作為一個矩陣的列向量，信號與矩陣唯一對應。由于矩陣中的各個列向量是同一個音頻信號的不同幀，因此具有相同的稀疏性?；谛碌臅r域框架，就可將簡單稀疏近似推廣為SSA，并將向量形式的PCD擴展為矩陣形式的Joint-PCD，從而利用冗余(超完備)字典原子的不同線性組合，實現對一個音頻信號所有幀同時進行降噪處理。

對于Z和D兩個變量，式(12)并不是一個凸優化問題，只有當其中一個變量給定時，它才是凸優化[3]。本文中超完備字典D是給定的，通過優化Z來構造Joint-PCD算法。式(12)還可以推廣到多個(大于2個)信號矢量的情況。

對于多個觀測信號的矩陣模型：

式中，X∈RN×P是包含了P個向量(音頻幀)x i∈RN的干凈信號矩陣；V∈RN×P是P個噪聲向量v i∈RN組成的矩陣；Y∈RN×P是含噪聲的觀測信號矩陣?；诼摵舷∈璞硎?JSR)[3]，干凈信號矩陣X可以由超完備字典D表示為X=DZ，則含噪聲的觀察信號矩陣為：

式中，||·||F是矩陣的Frobenius范數。如果系數矩陣Z中的非零行向量能正確地辯識出z i中的非零元素，就可利用線性的方法估計z i值。因此，可以把向量形式的降噪優化式(9)擴展為如下的矩陣形式：

式中，ρ(·)是任意標量函數，而ρ(Z)是衡量矩陣Z稀疏性的懲罰函數。式(16)是Joint-PCD算法的目標函數。

綜合考慮式(15)和式(16)，可以歸納出構造Joint-PCD算法的基本思想為：通過最小化Z中非零行的數量來限制參與對Y進行逼近的字典原子的總數。采用兩方面的措施來實現。首先，希望參與對Y近似的原子數量越少越好；其次，希望每個參與近似的原子要盡可能多地對矩陣Y的各列向量做出貢獻。對于系數矩陣Z，希望它的大多數行向量都是全零元素，而它的每個非零行(向量)都應包含盡可能多的非零元素?；谶@種考慮，將函數ρ(Z)設置為行?l0準范數(row-l0quasi-norm)的松弛版：式中，||Z||rx范數是把矩陣Z每一行中的最大絕對值項相加。本質上，將l∞范數應用于Z的行以促進它的非稀疏性，而將l1范數應用于所得向量z i以促進其稀疏性。若矩陣Z退化成一個列向量，則松弛范數降為l1范數，此時SSA的松弛條件與簡單稀疏近似的松弛條件相同。

把式(17)代入式(16)中，得到信號矩陣的降噪優化問題:

同樣地，從優化式(12)可知，函數ρ(Z)是矩陣Z的l1范數，本文對ρ(Z)的l1范數使用偶對稱的函數：

式中，zij是矩陣Z的第(i,j)個元素。另外，式(12)是將兩個列向量推廣到多個列向量的情形，即式(12)中的系數矩陣一般是P個列向量組成的矩陣Z∈RK×P，而不僅僅是由2個列向量組成的矩陣Z∈RK×2。

為方便起見，把||Z||rx稱為矩陣Z的rx-范數。求解式(18)的算法稱為rx-范數的Joint-PCD，求解式(12)的算法稱為l1范數的Joint-PCD。在給定系數矩陣的當前解Z k(其中k是迭代次數)條件下，定義矩陣U(D,Y,Z k)為：

式中，μ 是迭代步長；λ是收縮函數S(·)的閾值；從啟發式收縮算法[11]的推導過程可知，式(23)中的收縮函數S(·)是一個軟閾值函數[11]，其定義為：

式中，θ為軟收縮函數S(·)的閾值。

式(22)中的步長參數μ對算法的性能影響很大；對足夠小的μ >0，步長必須導致懲罰函數ρ(·)的可行下降。本文使用線搜索[10]尋找最佳μ值。線搜索作為一維優化尋求迭代步長，用于解決以下的優化問題：

式中，ρ(·)可以是rx-范數，也可以是l1-范數。

如果信號矩陣Y退化為一個列向量y，系數矩陣Z也僅僅是一個列向量z，則Joint-PCD算法將退化為PCD算法。因此，將Joint-PCD稱為一種啟發式算法。

在音頻信號處理中，經常使用兩個超完備字典[16]，即離散余弦變換(DCT)字典和Gabor字典。

一個加窗的Gabor字典為D G=[d(k,φ)](k,φ)∈Ψ，其原子以連續的參數集Ψ=[1,K]×[0,2π]為索引，其定義為：

式中，1≤k≤K，0≤t≤N?1，K是Gabor字典的大??；φ是相位；wd是加權窗口。d(k，φ)的分解可以用離散字典中的原子對表示：

以上兩式中，成對的原子是相同頻率且相位為零的余弦和正弦對。為簡化計算，詞典的列都做歸一化處理，即Gabor字典的正弦和余弦原子都使用其單位范數的版本。

輸入：觀測矩陣Y;字典D;閾值λ;函數ρ(·);迭代次數M

開始：設置k=1，迭代開始

If:k

誤差：計算E=Y-DZ k

映射：計算E T=DTE

收縮：計算E S=S(E T+Z k,λ)

線收縮：尋找μ0以獲得Zk+μ (E S-Z k)的下降

3 仿真實驗及結果分析

3.1 仿真環境和性能指標

通過仿真對矩陣形式的Joint-PCD與向量形式的PCD算法的性能進行全面比較，各算法的仿真實驗環境完全相同。

測試用的5個音頻音樂信源取自著名的音頻網站[17]。源信號類型為wav，采樣率為44 100 Hz。每個信源的樣本數量為L=216=65 536。5個音頻信源矩陣為X=[x1,x2,x3,x4,x5]T，其時域的波形如圖1所示。

圖1 5路音樂音頻源信號的波形

仿真實驗的硬件平臺：筆記本電腦，其CPU為Intel (R)Celeron (R)1007U CPU-1.50 GHz，內存為4 GB，PC機操作系統為Windows10。軟件平臺為MATLAB9(R2016a)。

降噪性能通過信噪比(signal to noise ratio,SNR)和均方根誤差(root mean square error,RMSE)來衡量：

式中，x為干凈信號；x?為降噪后信號；L為x的樣本數。

3.2 參數設置及仿真過程

對信源X施加一定強度的高斯白噪聲后形成含噪聲的觀測信號Y，噪聲強度一般采用高斯隨機變量的標準差σ來度量。顯然，信號Y中噪聲強度的大小對降噪算法的性能帶來很大影響。另外，在PCD和Joint-PCD算法中，軟收縮函數S(·)的閾值λ以及所采用的超完備字典(DCT或者Gabor)也都會對降噪性能造成一定影響。這幾個參數的最佳值分別為：噪聲強度σ=0.1；軟閾值λ=0.1；選用Gabor字典，其尺寸為N×K=64×1 024，加權窗為N＝64的對稱正弦窗。

在PCD和Join-PCD算法中，懲罰函數ρ(·)的形式、算法迭代步長μ也都需確定。不同形式的ρ(·)使用不同的矩陣(向量)范數。對于PCD，采用向量的l1范數來求解優化式(7)；對于Joint-PCD，使用矩陣的l1范數來求解優化式(12)，使用矩陣的rx-范數來求解優化式(18)。對于PCD和Joint-PCD算法的迭代步長μ (0<μ<1)，在算法的運行過程中使用線搜索來獲得最佳值，即在每次迭代中，用MATLAB函數fminbnd在線搜索μ值作為更新的步長。

音頻信源與高斯噪聲相加生成含噪聲信號，其降噪仿真過程如下。首先，每個含噪聲(時域)信號乘以長度為64的漢明窗生成語音幀，并使相鄰兩幀之間的重疊樣本數為32，即每幀的前后各有16個樣本與相鄰幀的樣本重疊，分幀使用Voicebox[18]中的MATLAB函數enframe實現。然后，用PCD和Joint-PCD算法分別對分割后的含噪信號實施降噪。執行PCD是對音頻信號的一幀進行降噪；Joint-PCD執行一次是對一個信號(而不僅僅是一幀)進行降噪。最后，重建降噪后的信號，經過降噪的音頻信號乘以逆漢明窗，在每一幀與相鄰兩幀重疊的樣本中，分別在兩邊各去除16個樣本，僅保留每一幀中間部分的樣本；將經過處理的各個音頻幀簡單串聯起來即可得到重新合成的音頻信號。

3.3 性能分析

對PCD和Joint-PCD算法的性能分析從兩個方面進行：1)算法在相同條件下所消耗的運行時間成本；2)算法降噪性能指標(SNR和RMSE值)的比較。

為了全面分析PCD和Joint-PCD算法的性能，在第一個仿真中，設置算法的最大迭代次數M=50。為便于表示，把Joint-PCD兩種范數對應的算法分別表示為Joint-PCD with rx norm和Joint-PCD withl1norm。算法Joint-PCD和PCD的運行時間指標如圖2所示。

圖2 算法的運行時間比較

從圖2可以看到，隨著迭代次數的增加，算法的運行時間幾乎都是呈線性增長，但不同算法對應的直線斜率(增長率)是不同的。對于Joint-PCD with rx norm和Joint-PCD withl1norm算法，隨著迭代次數的增加，它們的運行時間增長率幾乎是相同的；PCD算法運行時間增長率卻遠大于Joint-PCD的增長率。顯然，與PCD算法相比較，Joint-PCD算法能夠極大地降低運行時間。

為了便于比較PCD和Joint-PCD的降噪性能，在每次迭代后，分別計算每種算法對每個信號的降噪指標值。圖3給出了PCD算法的輸出SNR和RMSE波形。

圖3 PCD算法對每個信號的輸出SNR和RMSE指標

同樣地，圖4和圖5分別給出了Joint-PCD withl1norm和Joint-PCD with rx norm算法的輸出SNR和RMSE波形。

圖4 Joint-PCD with l1 norm算法對每個信號的降噪指標

比較圖3、圖4和圖5可知，對于每個信號，PCD算法的降噪性能幾乎與Joint-PCD with rx norm和Joint-PCD withl1norm算法的性能完全相同。經過25次迭代，PCD和Joint-PCD算法對每個信號的輸出SNR和RMSE值都逐漸收斂到其穩定值。這表明從降噪性能的角度來看，用Joint-PCD算法代替PCD算法是完全可行的。

圖5 Joint-PCD with rx norm算法對每個信號的降噪指標

為了更直觀地觀察PCD和Joint-PCD算法的平均降噪性能，在每次迭代后，計算并記錄了每種算法對于5個信號的平均SNR值。其結果如圖6所示。

圖6 5個信號的平均SNR指標

同樣地，PCD和Joint-PCD算法對5個信號的平均降噪性能指標RMSE的值如圖7所示。

圖7 5個信號的平均RMSE指標

從圖6和圖7可以看出，PCD和Joint-PCD算法的平均性能指標隨著迭代次數的增加具有一定的波動，特別是在算法的初始迭代階段，其波動更為明顯。這是因為每個信源的幅度值不同，而施加到各個信源上的高斯白噪聲的強度卻一樣，這就造成每個信號受到噪聲的干擾程度不同。另外，在降噪過程中，由于每個信號被噪聲污染的嚴重性不一樣，它們輸出的SNR和RMSE值(如圖3～圖5所示)就具有不同的波動范圍。當把這5個信號的降噪性能平均之后，其對應的平均SNR和RMSE值有波動就是一種很自然的現象，不可能消除這種平均性能的波動。將5個信號的降噪性能進行平均是為了研究算法穩定收斂的條件。從圖6和圖7還可看出，對于PCD和Joint-PCD with rx norm算法，平均SNR和RMSE值的波動在大約20次迭代后結束；對于Joint-PCD withl1norm算法，其波動在大約25次迭代后結束。這表明Joint-PCD和PCD算法在經過大約25次迭代后能穩定收斂。

為了更精細地分析不同算法的降噪性能，通過第二個仿真分析PCD和Joint-PCD算法的差異性。從上述的結果可知，PCD和Joint-PCD算法在經過25次迭代后能穩定收斂，因此，在這個仿真中，將算法的最大迭代次數設置為M=25。比較PCD和Joint-PCD算法達到穩定收斂所需要的運行時間，其結果如表1所示。

從表1可以看出，Joint-PCD withl1norm算法的運行時間最少，而Joint-PCD withrxnorm算法的運行時間比前者僅多了0.379 2 s，因此，兩種Joint-PCD算法的運行時間幾乎相同。與之相反，PCD算法收斂所花費的運行時間大約是Joint-PCD算法的4.6倍。這充分說明Joint-PCD算法比PCD算法具有更高的計算效率。

表1 PCD和Joint-PCD算法的運行時間s

對于PCD和Joint-PCD算法的降噪性能，采用每種算法對每個信號在整個迭代過程中輸出SNR和RMSE的平均值來度量，其結果如表2所示。表中粗體數字為相應的最佳降噪性能指標。

從表2中的數據可以看出，對于第2、第4和第5個信號，Joint-PCD withl1norm算法的平均SNR和RMSE值是最優的；對于第1和第3個信號，PCD算法的平均指標是最優的。另外，從表2還可以看出，PCD和Joint-PCD算法平均SNR值的差異僅發生在小數點后面第二位；而平均RMSE值的差異僅發生在小數點后面第四位。這表明，用Joint-PCD代替PCD，對音頻信號降噪性能幾乎不產生任何影響。

表2 每個信號在整個迭代過程的平均SNR和RMSE值

綜合表1和表2的結果可知，在降噪性能基本相同的情況下，Joint-PCD算法的計算復雜度大幅度降低。

4 結束語

本文定義了一種音頻信號的時域處理框架，即將每個被分割的音頻幀都作為一個列向量從而生成信號矩陣。在聯合稀疏表示(JSR)和同時稀疏近似(SSA)的基礎上，利用新的處理框架，提出了以信號矩陣為操作對象的Joint-PCD算法。仿真實驗的結果表明：Joint-PCD算法不僅能夠提供與向量形式的PCD算法幾乎完全相同的降噪性能，更重要的是，Joint-PCD算法提供了更高的計算效率。與傳統的PCD算法相比，Joint-PCD算法在大規模數據處理和實時信號處理應用方面具有巨大的潛在優勢。