文昊翔，洪遠泉，羅 歡

（韶關學院智能工程學院，韶關，512000）

引 言

自適應算法[1?2]在系統(tǒng)辨識領域有廣泛應用。而在某些系統(tǒng)辨識應用(如回聲消除、有源噪聲控制)，因系統(tǒng)延時不確定，導致目標系統(tǒng)的沖激響應(Impulse response, IR)序列通常具有稀疏性[3]。即序列雖然包含大量系數(shù)，但其中大部分系數(shù)為零值系數(shù)或極小值以模擬延時；真正產(chǎn)生能量的大幅值系數(shù)在時域聚集，且數(shù)量只占系數(shù)總量極小一部分。

傳統(tǒng)的最小均方(Least mean square, LMS) 算法與歸一化LMS(Normalized LMS, NLMS)算法不考慮目標系統(tǒng)的稀疏性，用一個長的濾波器辨識IR 序列。過長的濾波器導致自適應系統(tǒng)收斂速度下降，計算復雜度增加[4]。

利用稀疏性可有效解決或部分緩解濾波器過長導致的各種缺陷[5]。雙濾波器結構(Dual filter struc?ture, DFS)[6?7]是其中一種有效的解決方案。它用一級濾波器W1(k)以低分辨率辨識目標系統(tǒng)以定位活躍系數(shù)位置；然后一個短的二級濾波器W2(k)在時域覆蓋所有活躍系數(shù)，并精確辨識活躍系數(shù)。W2(k)通過排除大量非活躍系數(shù)參與自適應算法以降低濾波器長度。

目前對DFS 的研究主要針對W1(k)，需解決的問題是：在不降低定位精度的條件下，如何降低W1(k)的長度與計算復雜度。傳統(tǒng)的解決辦法是對W1(k)的輸入信號進行降采樣以降低W1(k)長度[8]。但降采樣將導致信號頻譜混疊，高頻分量丟失，不利于W1(k)定位活躍系數(shù)。

傳統(tǒng)的DCT、FFT 變換等將頻域局部化而時域信息完全丟失，因此無法應用于時域定位。小波變換[9]可以同時兼顧信號的時頻特性。Haar 小波作為最基本的小波，具有支撐區(qū)間不重疊，計算復雜度低等優(yōu)良特性。通過求取Haar 小波與輸入信號的內(nèi)積，可以同時提取信號的時、頻信息。Haar 變換為延時估計、自適應辨識稀疏系統(tǒng)提供了新思路[10]。

針對回聲消除應用，Ho 等[11]較早將Haar 應用于辨識稀疏系統(tǒng)。受此啟發(fā)，Bershad 等結合DFS，提出Haar 子帶變換 (Partial Haar transform, PHT)[12]算法。PHT 先在某個尺度對W1(k)的輸入信號進行PHT 變換以壓縮信號長度，然后再進行自適應濾波，并且Bershad 等從理論上證明W1(k)將最終收斂于IR 序列在該尺度下的Haar 變換系數(shù)。

以PHT 為基礎，Kechichian 等[13]引入模糊邏輯判定活躍系數(shù)位置；Ribas 等[14]則嘗試進一步降低算法的計算復雜度。在眾多學者的努力下，PHT 已經(jīng)成為一類自適應辨識稀疏系統(tǒng)的有效算法。

本文提出改進PHT (Improved?PHT, I?PHT)算法。新算法以PHT 為基礎，通過為W1(k)引入時變收斂步長并優(yōu)化該步長使后驗誤差最小化以加快W1(k)收斂速度；然后通過分時保存、計算W1(k) 的歸一化因子以降低算法的計算復雜度；最后以回聲消除為應用背景對新算法進行實驗仿真以驗證新算法的有效性。

1 問題描述與相關工作

1.1 自適應系統(tǒng)辨識

目標系統(tǒng)IR 序列記為H=[h1,h2,…,hN]T，其中N為序列長度。記k時刻輸入信號向量為X(k)=[x(k-1),x(k-2),…,x(k-N)]T，向量X(k)包含最近N個時刻的輸入信號。X(k)與H卷積，卷積結果加上背景噪聲v(k)形成系統(tǒng)輸出d(k)。

系統(tǒng)辨識需解決的問題為：如何根據(jù)X(k)與d(k)估計H。根據(jù)自適應理論，可以用一個FIR 濾波器W(k)=[w1(k),w2(k),…,wN(k)]以自適應算法更新其系數(shù)以逼近、辨識目標系統(tǒng)的IR 序列。其中LMS 算法系數(shù)更新方程為

式中：μ為收斂因子，用于平衡算法的收斂速度與辨識精度。濾波器W(k)按自適應算法更新其系數(shù)，經(jīng)過若干次迭代，W(k)將收斂于目標系統(tǒng)的IR 序列。若記E[]為取數(shù)學期望，則有

式(4)為基本的自適應理論，與之相關的詳細內(nèi)容可參考文獻[1]。

圖1 典型稀疏系統(tǒng)的IR 序列Fig.1 Impulse response of typical sparse sys?tem

1.2 稀疏性與DFS

為保證充分建模，W(k)的長度必須大于等于H的長度。但當目標系統(tǒng)具備稀疏性時，若能獲得系統(tǒng)的延時估計，只需精確辨識活躍系數(shù)即可完成對整個目標系統(tǒng)的自適應辨識。一個典型稀疏系統(tǒng)的IR 序列如圖1 所示[6]。它是電子回聲路徑的IR 序列，采樣頻率8 kHz。

DFS 可利用稀疏性提高算法效率，其原理如圖2[6]所示。W1(k)以低分辨率辨識目標系統(tǒng)并進行延時估計；然后一個短的W2(k)在時域覆蓋所有活躍系數(shù)，并精確辨識活躍系數(shù)。

1.3 Haar 子帶變換

離散Haar 小波母函數(shù)定義為

離散Haar 小波定義為

圖2 雙濾波器結構Fig.2 Dual filter structure

式中：t僅取整數(shù)值；m為尺度因子，取值范圍1≤m≤log2N，N為輸入信號長度，m越大，對應的Haar 時間窗越寬，以利于提取信號低頻分量，m越小，則時間窗越窄以提取高頻分量；n為時移因子，通過調(diào)整n平移小波母函數(shù)以覆蓋不同的時域范圍。

為便于闡述，以輸入信號長度N=8 為例說明。Haar 變換矩陣定義如下

式中：第1 列對應m=3，與輸入信號求內(nèi)積后提取信號低頻分量；第2、3 列對應m=2，分別提取信號時域前半、后半部分頻率信息；第4～7 列對應m=1，分別提取信號對應4 部分的高頻分量。

由式(7)可得Haar 子帶變換矩陣定義如下

利用Haar 子帶變換，Bershad 等[12]提出PHT 算法，W1(k)在Haar 域?qū)π盘栠M行自適應濾波，以降低濾波器長度。系數(shù)更新過程為

同樣以N=8 為例，Haar 尺度因子M定義為M=log2N=3。當取m=3 時，按照式(11) PHT 變換的定義，式(8)的變換矩陣將與輸入信號向量X(k)進行矩陣相乘，W1(k)長度N1壓縮為N1=1；若m=2，則N1=2；若m=1 時，則N1=4?？傻?，式(11)的 PHT 變換可大幅壓縮W1(k)的長度。

經(jīng) Bershad 等證明[12]，W1(k)按式(11—13)進行自適應算法，收斂后有

因此PHT 算法既可壓縮W1(k)的長度，又可以獲得與目標系統(tǒng)H相關的信息。以W1(k)的信息為基礎，再利用峰值系數(shù)位置[13]、模糊邏輯[14]或移動窗積分法[15]即可進行延時估計以定位活躍系數(shù)。

式(5—14)為PHT 算法的基本原理，與之相關的詳細理論推導可參考文獻[12]。

2 改進PHT 算法

以PHT(式(11—13))為基礎，提出I?PHT 算法提高W1(k)收斂速度。新算法先為PHT 引入時變收斂步長β(k)，然后以后驗誤差最小化為目的優(yōu)化β(k)。

2.1 推導過程

以式(13)為基礎引入時變收斂步長β(k)可得

自適應算法的先驗誤差e1(k)定義如式(12)。后驗誤差ε1(k)定義為

新算法目的為求取β(k)最優(yōu)值，使式（17）極大化。

其物理意義為在先驗誤差e1(k)確定的條件下，使后驗誤差ε1(k)能量最小化。以下為化簡f(β(k))過程。由式(12)可得

同理，由式(16)可得

將式(15)代入式(19)，化簡、整理后得

式(18)與式(20)兩式相減可得

由式(21)可得，f(β(k))是關于變量β(k)的二次函數(shù)。對函數(shù)f(β(k))關于β(k)求導，當導數(shù)為零時，f(β(k))取得極大值，此時β(k)即為最優(yōu)解。即

解得

將式(23)代入式(15)可得

再次引入收斂步長μi平衡算法收斂速度與穩(wěn)態(tài)失調(diào)。最后，新算法的系數(shù)更新方程為

2.2 計算復雜度優(yōu)化

I?PHT(式(24))與 PHT(式(13))相比，主要是式(24)的分母 [ZTm(k)Zm(k)](即歸一化因子)，引入額外的計算復雜度。

在NLMS 算法中，因X(k)與X(k-1)是簡單的線性移位關系，可以利用式(25)快速計算歸一化因子。

但在I?PHT 中，因Zm(k)與Zm(k-1)不再滿足線性移位關系，因此[(k)Zm(k)]不能直接套用式(25)快速計算獲得，需另外尋找新的高效計算方法。

為方便討論，記Zm(k)=[z1(k),z2(k),…,zN1(k)]。定義壓縮比為r=N/N1，式中N為X(k)長度，N1為Zm(k)的長度。通過PHT 變換(式(11))，X(k)中相鄰的r個信號被壓縮為Zm(k)中的1 個信號。因此雖Zm(k)與Zm(k-1)不再保持線性移位關系，但Zm(k)與Zm(k-r)仍保持線性移位關系。針對該關系，建立向量P= [p1(k),p2(k),…,pr(k)]。記i=1+(kmodr)，式中mod 為求余運算。將(k)Zm(k)值保存在pi(k)中，有pi(k)=(k)Zm(k)。即把算法歸一化因子按不同時刻分別保存在向量P中，并以r為循環(huán)。此時有

至此，新算法的實現(xiàn)可歸納為按步驟執(zhí)行式(11,12)與式(27—29)。

與 PHT 算法相比，I?PHT 每次迭代僅額外引入 3 次加/減法、2 次乘法。

3 實驗仿真

以回聲消除為應用背景進行實驗仿真以檢驗新算法的性能。目標系統(tǒng)H采用圖1 所示的實際回聲路徑，長度為N=512?；钴S系數(shù)定位的判斷邏輯采用移動窗積分法[15]。

壓縮比r決定W1(k)的長度N1。若r取值過小，則不能有效縮短W1(k)的長度；若r取值過大，則分辨率過低，無法有效定位活躍系數(shù)。本實驗按文獻[13]取r=4，即N1=128。

實驗仿真分兩部分進行，實驗1 檢驗新算法對W1(k)性能的影響，實驗2 檢驗新算法對整個DFS 的影響。

3.1 實驗1

以方差為1 的高斯白噪聲作為輸入信號X(k)。X(k)與H卷積，再以另一高斯白噪聲作為背景噪聲v(k)，XT(k)H與v(k)二者信噪比為 20 dB。

實驗 1 主要比較 PHT 與 I?PHT 兩種算法。對 I?PHT，μi=0.5；對 PHT，μp=μi/E[(k)Zm(k)]。其中E[(k)Zm(k)]根據(jù)輸入信號的先驗知識統(tǒng)計獲得。對W1(k)的性能評價采用兩個指標：失調(diào)與有效定位活躍系數(shù)所需迭代次數(shù)。

(1)失調(diào)，單位為dB，計算公式為

式(30)中H1的計算公式為H1=ψM,m H。

兩算法的W1(k)失調(diào)曲線如圖3 所示。由圖3 可得，I?PHT 與PHT 相比，穩(wěn)態(tài)失調(diào)相似，但收斂速度有明顯提高。

（2)定位活躍系數(shù)所需迭代次數(shù)。當W2(k)開始覆蓋目標系統(tǒng)的峰值系數(shù)時，可認為W1(k)已準確定位活躍系數(shù)。由圖4 所展示的實驗得，為準確定位活躍系數(shù)，PHT 約需408 次迭代，而 I?PHT 僅需約 101 次迭代。

圖3 一級濾波器失調(diào)曲線對比Fig.3 Misalignment curve comparison of the first-order filter

3.2 實驗2

因W1(k)僅用于低精度辨識目標系統(tǒng)以定位活躍系數(shù)，精確辨識目標系統(tǒng)還需W2(k)完成，所以W2(k)的性能改善才是整個自適應系統(tǒng)的最終目標，對W2(k)性能的評估與W1(k)相比，更顯重要。

在計算W2(k)算法失調(diào)時，DFS 需先將短的W2(k)在首尾補零擴展至長度為N=512，記擴展后的序列為W(k)，再按式（31）計算失調(diào)。

重復實驗 1 的參數(shù)設置，PHT 與 I?PHT 的W2(k)均采用LMS 算法更新系數(shù)，W2(k)的長度N2=64。μ=0.077。兩算法的W2(k)失調(diào)曲線如圖4 所示。

結合圖4 與實驗 1 可得，因 I?PHT 可使W1(k)獲得更快的收斂速度，在更短的時間內(nèi)準確定位活躍系數(shù)，所以最終導致W2(k)的收斂速度亦有大幅提高。因PHT 的W1(k)約需400 次迭代才可獲得穩(wěn)定的延時估計，因此其W2(k)只有在400 次迭代之后才進入快速收斂階段。而I?PHT 的W2(k)僅需約100 次迭代即進入快速收斂階段，獲得比PHT 更快的整體收斂速度。

為更好地模擬回聲消除系統(tǒng)工作環(huán)境，以實際語音信號激勵目標系統(tǒng)，并再次進行實驗仿真。系統(tǒng)信噪比為20 dB。為檢驗算法跟蹤能力，在20 s 時，目標系統(tǒng)IR 序列整體左移 200 個單位。針對W1(k)的收斂步長，I?PHT 設為μi=0.1；PHT 設為μp=μi/E[(k)Zm(k)]。W2(k)的收斂步長，I?PHT 與 PHT 均統(tǒng)一設為μ=0.03。單濾波器 LMS算法作為經(jīng)典的自適應算法亦納入對比，濾波器長度設為N=512，步長設置為μ=0.03/8。其他設置與實驗1 一致。3 種算法的W2(k)失調(diào)如圖5 所示。

圖4 二級濾波器失調(diào)曲線對比Fig.4 Misalignment curve comparison of two-stage filter

圖5 3 種自適應算法失調(diào)曲線對比Fig.5 Misalignment curve comparison of three adaptive algorithms

由于DFS 能有效降低濾波器長度，其W2(k)長度僅為LMS 算法W(k)長度的1/8。隨著濾波器長度的下降，算法的收斂速度迅速提高。由圖5 得，I?PHT、PHT 的收斂速度顯著高于LMS 算法。

而同為 DFS，針對W1(k)而言，I?PHT 優(yōu)于 PHT 的地方主要有以下 3 點：

(1) 因I?PHT 的W1(k)能更快地定位活躍系數(shù)，因此其整體收斂速度明顯優(yōu)于PHT。圖3—5 可證明該結論。

(2) 實驗發(fā)現(xiàn)當輸入信號是真實語音信號時，I?PHT 的W1(k)收斂步長μi取值范圍大于PHT 的μp。取μi<0.8 仍能保證W1(k)收斂，但 I?PHT 必須滿足μp<0.1/E[ZTm(k)Zm(k)]才能保證W1(k)收斂。因此 I?PHT 的穩(wěn)定性高于 PHT 算法。

(3) I?PHT 的步長μi取值范圍有較清晰的指導值，其理論取值范圍為0<μi<1。當輸入信號為高斯噪聲時，該理論范圍可保證W1(k)收斂；當輸入為強相關、非平穩(wěn)的語音信號時，各種自適應算法的穩(wěn)定性均會減弱，此時I?PHT 的μi實驗取值范圍約為0<μi<0.8。但PHT 算法需先獲得E[ZTm(k)Zm(k)]的先驗知識，才能確定μp的取值范圍。因此，與PHT 算法相比，I?PHT 更易于實現(xiàn)。

4 結束語

針對稀疏系統(tǒng)辨識應用，DFS 是一類有效的解決方案。它用W1(k)定位活躍系數(shù)位置，用一個短的W2(k)精確辨識活躍系數(shù)。它通過降低濾波器有效長度，以達到提高收斂速度，降低計算復雜度的目的。

本文主要針對W1(k)進行討論，提出I?PHT 算法。新算法先對W1(k) 的輸入信號進行PHT 變換以壓縮信號長度；然后引入時變步長，并以后驗誤差最小化為目標函數(shù)優(yōu)化時變步長；最后將新算法的歸一化因子分時、循環(huán)保存到一個向量中，并通過自回歸方式維護歸一化因子以降低計算復雜度。

以回聲消除為應用背景對I?PHT 進行實驗仿真，仿真結果表明，DFS 能有效降低濾波器長度以提高收斂速度。而與PHT 相比，I?PHT 算法收斂速度與穩(wěn)定性均有明顯提高。