閆曉瑾，何國榮，梅鐵民

（沈陽理工大學自動化與電氣工程學院，遼寧 沈陽 110159）

0 引言

信號解卷積是信號處理領域中一個重要的研究內容，其廣泛應用于圖像處理、語音信號處理、陣列信號處理、地震波信號處理、生物醫(yī)學、故障診斷等領域。在這些領域中，通常觀測信號是源信號與一線性系統(tǒng)沖激響應的卷積，解卷積就是從觀測信號中恢復源信號或系統(tǒng)沖激響應的過程。

已有的解卷積算法大致可分為兩類，一類是房間沖激響應未知情況下的盲解卷積算法，另一類是房間沖激響應已知的解卷積算法。盲解卷積是指在房間沖激響應未知或不需要估計的情況下，由觀測信號直接估計解卷積濾波器。Bussgang 算法[1]是自適應信號處理領域中最經(jīng)典的盲解卷積算法，該算法實現(xiàn)簡單且收斂速度快，但該算法需要源信號和解卷積噪聲的統(tǒng)計特性等先驗知識。針 對SISO（Single Input Single Output）系 統(tǒng)，文獻[2]將觀測信號轉化為復基帶信號，然后利用最小均方算法（Least Mean Square,LMS）來實現(xiàn)盲解卷積，該算法適用于高斯信號和非高斯信號，對源信號的適用面較寬，但仍然存在收斂速度和計算復雜度之間的矛盾問題。在生物醫(yī)學中，心電圖信號、脈搏信號等可以被認為是心跳信號與人體系統(tǒng)的卷積，通過對心電圖信號進行解卷積處理即可獲得病人的心跳信號，文獻[3]中利用遞歸最小二乘算法（Recursive least squares,RLS）來求解卷積濾波器，利用定點卷積核補償（Fixed-point convolution kernel compensation,FP-CKC）求出RLS算法的期望信號，利用粒子群優(yōu)化算法（Particle Swarm optimization,PSO）求出逆濾波器的初始值,該算法可用于較高的噪聲水平。在語音去混響領域中，解卷積過程是其中重要的一部分。文獻[4]根據(jù)聲源語音信號的短時傅里葉變換系數(shù)的稀疏性，將非負矩陣分解（Non-negative Matrix Factorization,NMF）應用到多通道線性預測（Multi-Channel Linear Prediction,MCLP）中來達到去混響的目的，該算法取得不錯的效果，不足之處在于明顯增加了計算量。在軸承故障診斷領域中，振動信號是由故障引起的周期性沖擊信號與機械部件的共振響應卷積的結果，傳統(tǒng)的解卷積方法有最大相關峭度解卷積算法（Maximum Correlated Kurtosis Deconvolution,MCKD），該算法的解卷積信號的包絡譜中可以呈現(xiàn)出典型的周期性脈沖的特點[5]。文獻[6]提出最小熵解卷積算法（Mimum-Entropy Deconvolution，MED），該算法可適用于非最小相位系統(tǒng)，但要求源信號具有簡單稀疏的統(tǒng)計特性。除此之外，還有基于稀疏表示的盲解卷積算法[7]、子空間分解的盲解卷積方法[8]等。在語音去混響解卷積算法中，通常先估計出房間沖激響應，然后再進行解卷積。在SISO 系統(tǒng)中，文獻[9]利用復倒譜技術估計出房間沖激響應，并在倒譜域內求其逆濾波器，然后將其逆濾波器作為改進的自然梯度算法的初始值來進行迭代解卷積。MINT（the multiple-input/output inverse-filtering theorem）定理是Miyoshi 等人提出的在房間沖激響應已知的情況下求逆濾波器的一種經(jīng)典方法[10]，文獻[11]和[12]利用觀測信號的自相關矩陣來獲得房間沖激響應的先驗知識，然后根據(jù)MINT 定理來求逆濾波器，但該算法要求源信號是平穩(wěn)的白噪聲信號。其實，利用MINT 定理求逆濾波器的方法在實際中并不可行，因為該方法對噪聲非常敏感，針對此問題，文獻[13]提出在子帶最小二乘算法中使用正則化，該算法可以降低逆濾波器對房間沖激響應的估計誤差和對觀測噪聲的敏感度。文獻[14]通過觀測信號的二階累積量估計出非最小相位系統(tǒng)的房間沖激響應，對其卷積矩陣直接求逆，該方法在高階系統(tǒng)時的計算量將非常大。文獻[15]在估計出房間沖激響應后，利用p-范數(shù)和窗函數(shù)來求逆濾波器，進而實現(xiàn)解卷積。文獻[16]給出了在估計出房間沖激響應的情況下不求逆濾波器的解卷積算法，利用卡爾曼濾波算法直接進行解卷積，該算法具有較高的噪聲穩(wěn)定性,但其缺點是計算量大。

本文在卡爾曼濾波解卷積算法的基礎上，提出了一種降低計算量的優(yōu)化算法，該算法在保證算法穩(wěn)定性的同時，又能大幅度減少計算量。該算法不僅可以用于語音去混響，還可以應用在其他領域。

1 標準的卡爾曼濾波解卷積算法

1.1 數(shù)學模型

在單輸入多輸出系統(tǒng)（Single Input Multiple Output,SIMO）中，觀測信號是聲源信號和房間沖激響應的卷積，那么，觀測信號可以表示為：

其中，s(n)表示聲源信號，n為時間，hi表示第i路房間沖激響應，xi(n)表示第i路麥克風接收的觀測信號，vi(n)表示第i路觀測噪聲，L為房間沖激響應的長度，N為麥克風的數(shù)量。

利用（1）式建立卡爾曼濾波的狀態(tài)方程和觀測方程：

狀態(tài)方程：

測量方程：

其中S(n)為n時刻的狀態(tài)矢量，由聲源信號s(n)構成狀態(tài)矢量S(n)=[s(n),s(n-1),…,s(n-L+1)]T；X(n)=[x1(n),x2(n),…,xN(n)]T為n時刻的觀測矢量；觀測矩陣由N路房間沖激響應構成，是一個N×L維矩陣（其中hi=[hi(0),hi(1),…,hi(L-1)]T）；u1(n)和u2(n)分別為均值為零的過程白噪聲和觀測白噪聲，它們的協(xié)方差矩陣分別為和（I為單位矩陣）；W(n+1,n)為L×L維狀態(tài)轉移矩陣，定義為：

并且，W(n,n+1)=WT(n+1,n)。

1.2 標準的卡爾曼濾波解卷積算法

已知房間沖激響應（h1,h2,...hN），在一個解卷積塊中的標準卡爾曼濾波解卷積算法的迭代過程列于表1 中。

表1 標準卡爾曼濾波解卷積算法

其中G(n)為卡爾曼增益，P(n)為濾波的狀態(tài)矢量估計誤差的自相關矩陣，K(n+1,n)為狀態(tài)矢量一步預測誤差的自相關矩陣，α(n)為新息?？柭鼮V波器輸出的n時刻的聲源信號的最后一個分量。

標準卡爾曼濾波算法中五個變量迭代一次的計算量統(tǒng)計于表2 中。關于算法的計算量問題，本文只統(tǒng)計乘除法的運算次數(shù)，不統(tǒng)計加減法。從表2中可以看出，標準卡爾曼濾波算法的計算量主要集中在G(n)、P(n)和K(n+1,n)上，相比之下，和α(n)的計算量可以忽略不計，所以要想減少計算量，需從G(n)、P(n)和K(n+1,n)入手。

表2 標準卡爾曼濾波算法的計算量

2 改進的卡爾曼濾波解卷積算法

仔細研究標準卡爾曼濾波解卷積算法發(fā)現(xiàn)，G(n)、P(n)和K(n+1,n)三者的迭代計算在極短的時間內即可達到收斂，并且G(n)只與H(n)有關，與觀測信號無關；而α(n)、除與H(n)有關外，還與當前時刻的觀測信號X(n)有關。對于非時變系統(tǒng)而言，G(n)一旦收斂，就不再變化，因此可以把與G(n)計算有關的項提取出來單獨計算，待其收斂后再進行α(n)和的計算，從而實現(xiàn)解卷積。因此，可以通過把標準卡爾曼濾波解卷積算法分解成兩個子循環(huán)的辦法來到達降低樣本平均計算量的目的。新算法的迭代步驟列于表3 中。

表3 改進的卡爾曼濾波解卷積算法

通過改變標準算法的迭代次數(shù)和迭代順序，使其在保證穩(wěn)定性的同時，又能明顯減少計算量。在卡爾曼增益G(n)的計算中，Bg的選取要具體問題具體分析（例如，若算法要求計算量小，則Bg盡量??；若算法要求收斂性，則Bg可以大一些）。該算法通過犧牲收斂速度來降低算法的計算量，但犧牲的收斂速度有限。

在分塊解卷積算法中，認為在一個塊內系統(tǒng)沖激響應是不變的，假設一個解卷積塊的長度為B，則從上述迭代過程中可以看出，改進算法比標準算法減少的計算量為(B-Bg)(4L3+4L2N+2LN2+N3)。

3 仿真實驗

在仿真實驗中用到的數(shù)據(jù)是實測的房間沖激響應（信道數(shù)N為4，長度為1300 點）和一段語音信號（長度為21000 點）。把已知的語音信號作為聲源信號s(n)，然后將聲源信號與實測的房間沖激響應卷積（在不同的仿真實驗中，房間沖激響應的長度會有多不同）來獲得卷積信號。由于本文只討論算法計算量的問題，所以在仿真實驗中不加入噪聲。

令解卷積信號估計誤差e(n)=se(n)-s(n)的均方值為E(n)；聲源信號的均方值為Q(n)。采用如下的滑動平均公式在線估計：

其中，滑動因子α=0.9；n是迭代次數(shù)；e(n)為算法輸出誤差；se(n)為算法的解卷積信號；s(n)為聲源信號。

定義解卷積信號的誤差-信號比作為評價算法性能的標準：

Esir曲線能夠反映出算法的收斂水平，單位是分貝（dB）。

3.1 卡爾曼增益的收斂速度

在第2 節(jié)中提到，對于非時變系統(tǒng)而言，G(n)一旦收斂，就不再變化，而卡爾曼增益G(n)的計算只與房間沖激響應H(n)有關，所以本節(jié)的仿真內容是確定卡爾曼增益的收斂速度。在本節(jié)中房間沖激響應的長度分別取L=100，150，200（從已知的房間沖激響應中任意截取）。標準卡爾曼濾波算法中對于不同長度的房間沖激響應時卡爾曼增益的收斂情況見圖1。由于卡爾曼增益在全局的收斂速度較快，所以圖1 中的卡爾曼增益的收斂情況為前400點的放大圖。從圖1 中可以看出，卡爾曼增益的收斂速度幾乎與房間沖激響應的長度L相同。

3.2 改進的卡爾曼濾波算法對解卷積結果的影響

在本節(jié)的仿真實驗中，房間沖激響應長度L=150，在上一節(jié)的仿真實驗中發(fā)現(xiàn)卡爾曼增益的收斂速度幾乎與房間沖激響應的長度L相同，所以本節(jié)中選取改進算法的卡爾曼增益的迭代次數(shù)Bg=L+10=160。圖2（a）給出了仿真卷積信號（1）聲源信號（2）改進算法的解卷積信號（3）和標準算法的解卷積信號（4）的波形比較。其中仿真卷積信號是第2 路的卷積信號。從圖2 中可以看出兩種算法都有較好的解卷積效果，改進算法只是收斂速度慢一點，一旦收斂，二者的解卷積結果一致。

圖1 不同長度的房間沖激響應對應的卡爾曼增益的收斂情況（L=100，150，200）

圖2 聲源信號、卷積信號和兩種算法的解卷積信號的比較

3.3 改進的卡爾曼濾波算法對收斂速度的影響

本節(jié)與上一節(jié)實驗相同，房間沖激響應長度L=150，B=5L=750，Bg=L+10=160。利用兩種算法的解卷積信號來繪制出Esir 曲線，如圖3 所示。

圖3 標準算法解卷積信號和改進算法解卷積信號的Esir 曲線

從圖3 可以看出，雖然改進算法的解卷積信號沒有標準算法解卷積信號那么高的信噪比，但是改進算法的解卷積信號的Esir 也能夠快速達到50dB以下。從算法的運行時間來看，在仿真實驗所用的計算機上，本節(jié)實驗的標準算法的運行時間是5.92秒，改進算法的運行時間是1.56 秒，減少了73.6%的運行時間，說明計算量減少了70%以上。

4 結語

本文提出了一種基于卡爾曼濾波解卷積的優(yōu)化算法，利用卡爾曼增益在短時間內收斂的特點，將標準卡爾曼濾波算法用兩個子循環(huán)來完成，卡爾曼增益、狀態(tài)矢量一步預測誤差的自相關矩陣和狀態(tài)估計誤差的自相關矩陣的迭代作為一個子循環(huán)，狀態(tài)矢量和新息的迭代作為一個子循環(huán)。在仿真實驗中可以看出，與標準卡爾曼濾波解卷積算法相比，本文提出的改進算法犧牲的收斂速度有限，但降低的計算量是明顯的。