摘 要： 麥克風(fēng)陣列信號(hào)處理已經(jīng)廣泛引用于語(yǔ)音處理、通信等領(lǐng)域。時(shí)延的精確估計(jì)決定了麥克風(fēng)陣列聲源定位系統(tǒng)的性能。將重要抽樣的方法應(yīng)用到麥克風(fēng)陣列的時(shí)延估計(jì)中，可以保證時(shí)延估計(jì)精確性的前提之下，有效地回避了網(wǎng)格搜索和迭代計(jì)算，使結(jié)果不依賴于初始值，并且降低了計(jì)算的復(fù)雜度。仿真實(shí)驗(yàn)結(jié)果表明，重要抽樣時(shí)延估計(jì)方法在計(jì)算量降低一個(gè)數(shù)量級(jí)時(shí)，仍保持著與ML網(wǎng)格搜索法相近的性能，更適合麥克風(fēng)陣列時(shí)延的實(shí)時(shí)估計(jì)。

關(guān)鍵詞： 麥克風(fēng)陣列； 聲源定位； 時(shí)延估計(jì)； 重要抽樣

中圖分類號(hào)： TN911?34 文獻(xiàn)標(biāo)識(shí)碼： A 文章編號(hào)： 1004?373X（2015）24?0036?04

Time delay estimation of microphone array based on important sampling

YANG Bo1， LIU Zhicheng1， LI Demei2

（1. Engineering University of CAPF， Xi’an 710086， China； 2. Inner Mongolia PAP Corps of CAPF， Hohhot 010000， China）

Abstract： Microphone array signal processing has been widely used in the fields of speech processing and communication. The performance of microphone array sound source localization system is determined by the accurate estimation of time delay. If the important sampling method is applied to the time delay estimation of microphone array， the grid search and iterative calculation can be effectively avoided in the precondition of guaranteeing the accuracy of time delay estimation， so that the results do not depend on the initial value， and the computational complexity is reduced. The simulation results show that the time delay estimation method based on the important sampling keeps the performance similar to ML grid search method while the calculated amount is reduced by an order of magnitude， and is more suitable for the real?time estimation for the time delay of the microphone array.

Keywords： microphone arrays； acoustic source localization； time delay estimation； important sample

麥克風(fēng)陣列的聲源位置估計(jì)的問(wèn)題。實(shí)質(zhì)上就是對(duì)陣元之間接收信號(hào)時(shí)延的估計(jì)；因此，能否準(zhǔn)確和快速的估計(jì)出時(shí)延，決定了麥克風(fēng)陣列是否能夠?qū)崟r(shí)的、準(zhǔn)確的估計(jì)出聲源的位置。時(shí)延的估計(jì)在很多領(lǐng)域當(dāng)中都是一個(gè)重要問(wèn)題，例如雷達(dá)、聲納和無(wú)線通信網(wǎng)絡(luò)[1]。傳統(tǒng)方法采用匹配濾波和互相關(guān)的技術(shù)來(lái)進(jìn)行時(shí)延估計(jì)。眾所周知的最大似然時(shí)延估計(jì)ML方法，在時(shí)延估計(jì)方面展示出了良好的性能，尤其在高信噪比和信號(hào)特性緩慢變化時(shí)可以方差逼近著名的克拉美羅界（CRLB）下限。盡管如此，ML方法在計(jì)算的復(fù)雜度上仍存在一定缺陷，基于網(wǎng)格搜索的ML估計(jì)計(jì)算復(fù)雜度會(huì)隨著參數(shù)數(shù)量增加而急劇上升。當(dāng)然，采用迭代的方法成功的解決了ML的問(wèn)題，但是估計(jì)結(jié)果的準(zhǔn)確度過(guò)分依賴于初始值，如果沒有一個(gè)好的初始值，迭代的結(jié)果就會(huì)出現(xiàn)較大偏差[2]。采用ML方法對(duì)多個(gè)未知參數(shù)進(jìn)行估計(jì)時(shí)，引入重要抽樣[3]的概念可以有效地避開網(wǎng)格搜索和迭代所產(chǎn)生的計(jì)算復(fù)雜度問(wèn)題。

1 系統(tǒng)模型

假設(shè)麥克風(fēng)陣列總共含有P個(gè)陣元，每個(gè)陣元都將接收信號(hào)傳輸?shù)揭粋€(gè)中心節(jié)點(diǎn)進(jìn)行處理，麥克風(fēng)陣列的接收信號(hào)可以表示為[4]：

[yt=i=1Pαixt-τi+w（t）] （1）

式中：[xt]是聲源信號(hào)；[αi]是衰減因子，即第i個(gè)麥克風(fēng)的接收信號(hào)相對(duì)聲源信號(hào)的衰減系數(shù)，在[0，1]區(qū)間內(nèi)取值（近場(chǎng)模型中，[αi]各不相同；遠(yuǎn)場(chǎng)模型中，[αi]均相同，可以歸一化為1）；[τi]為第i個(gè)陣元接收信號(hào)同聲源信號(hào)之間的時(shí)延；[w（t）]是加性噪聲，為了簡(jiǎn)化模型，假設(shè)不同陣元間噪聲相互獨(dú)立與信號(hào)不相關(guān)。

在接收端，以[Fs]的采樣頻率對(duì)信號(hào)進(jìn)行采樣，采樣率[Fs=1Ts]，采樣信號(hào)可以表示為：

[ynTs=i=1PαixnTs-τi+wnTs， n=0，1，2，…，N-1] （2）

式中N為樣本總數(shù)。

將接收端樣本轉(zhuǎn)換到頻域，模型就可以用矩陣的形式表達(dá)。對(duì)式（2）進(jìn)行離散傅里葉變換，可以得到：

[Yk=i=1PαiXkexp-j2πkτiN+Wk] （3）

其中：[{Y（k）}N-1k=0]，[{X（k）}N-1k=0]和[{W（k）}N-1k=0]分別是[{y（nTs）}N-1k=0]，[{x（nTs）}N-1k=0]和[{w（nTs）}N-1k=0]的離散傅里葉變換。

將第1個(gè)陣元當(dāng)作參考陣元，其接收信號(hào)采樣的頻域表示為：

[Y1k=α1Xkexp-j2πkτiN+W1k] （4）

那么陣列接收信號(hào)與參考陣元接收信號(hào)之間的時(shí)延[Δτi=τi-τ1]，i=1，2，…，P，顯然[Δτ1=0]。為了突顯出待估計(jì)的參數(shù)，對(duì)式（3）進(jìn)行處理得到：

[Yk=i=1PβiY1kexp-j2πkΔτiN+Wpk] （5）

式中：

[βi=αiα1， i=1，2，…，P] （6）

[Wpk=Wk-βiW1（k）exp-j2πkΔτiN] （7）

然而式（5）滿足一般線性模型（GLM），可以將其轉(zhuǎn)化為如下形式：

[Y=ΦpΔτβ+Wp] （8）

式中：矩陣[Y]是麥克風(fēng)陣列所接收到的信號(hào)；矩陣[ΦpΔτ]僅取決于未知的時(shí)延[Δτ=[Δτ1，Δτ2，…，][Δτi…，Δτp]]。

[ΦpΔτ=ΦpΔτ1，ΦpΔτ2，ΦpΔτi，…，ΦpΔτp] （9）

向量[{ΦpΔτi}Pi=1]可以定義成：

[ΦpΔτi=Y10，Y11e-j2πΔτiN，Y12e-j2π2?τiN，…， Y1N-1e-j2π（N-1）?τiNT] （10）

其中：[Y1=[Y10，Y11，Y12，…，Y1N-1]T]是與參考麥克風(fēng)接收脈沖相一致包含傅里葉系數(shù)的樣本，[Wp=[Wp0，Wp1，Wp2，…，WpN-1]T]是加性噪聲成分。在這里并沒有必要保證噪聲嚴(yán)格服從高斯分布，噪聲樣本從連續(xù)時(shí)間過(guò)程得到之后，根據(jù)中心極限定理，每一個(gè)傅里葉系數(shù)[W（k）]都可以看作對(duì)這些樣本的一個(gè)簡(jiǎn)單加權(quán)，如果接收樣本數(shù)量足夠多，噪聲就可以被模型化為一個(gè)高斯隨機(jī)變量。那么，似然函數(shù)[5]可以表示如下：[ΛΔτ，β∝pY；Δτ，β= 1πMσ2MWexp -1σ2W×Y-Φp（Δτ）βH（Y-Φp（Δτ）β）] （11）

式中：[pY；Δτ，β]是[Δτ]和[β]參數(shù)化的概率密度函數(shù)；[σW]是[Wk]的方差，即噪聲功率。

傳統(tǒng)的ML估計(jì)方法是通過(guò)對(duì)式（11）搜索似然函數(shù)最大時(shí)所對(duì)應(yīng)的[Δτ]來(lái)得到[ΔτML]。然而，[ΛΔτ，β]由[Δτ]和[β]共同決定，進(jìn)行它們的聯(lián)合估計(jì)計(jì)算復(fù)雜度較大。因此，[有必要得出一個(gè)僅由未知量Δτ]決定的似然函數(shù)。為此可以考慮，將[β]表示為與[Δτ]相關(guān)的形式[β（Δτ）]。[β（Δτ）]通過(guò)由一個(gè)給定的[Δτ]來(lái)對(duì)對(duì)數(shù)似然函數(shù)[LΔτ，β=log {ΛΔτ，β}]搜索最大值所對(duì)應(yīng)的[β]得到。

將式（11）中[β]用[β（Δτ）]替換，然后將與[Δτ]無(wú)關(guān)的部分略去，就得到了如下的系統(tǒng)壓縮似然函數(shù)：[LcΔτ=1σ2WYHΦpΔτ（ΦpΔτHΦpΔτ）-1ΦpΔτHY] （12）

2 時(shí)延的估計(jì)

理論上通過(guò)找出[LcΔτ]最大值所對(duì)應(yīng)的[Δτ]即可估計(jì)出來(lái)未知時(shí)延。若直接用多維網(wǎng)格搜索的方法來(lái)求解這個(gè)問(wèn)題，計(jì)算復(fù)雜度也會(huì)隨著未知參數(shù)數(shù)量的增加呈指數(shù)方式增加。一些用來(lái)解決多維最優(yōu)化問(wèn)題的可選途徑展示出了相對(duì)好的結(jié)果，但是大部分途徑涉及到了迭代運(yùn)算。因此需要有一個(gè)良好的初始值，否則就無(wú)法保證收斂于全局最優(yōu)解。在這個(gè)情況之下，Pincus提出了一個(gè)可以確保得到多維函數(shù)全局最優(yōu)解的方法。根據(jù)[6]中提到的定理，代價(jià)函數(shù)[LcΔτ]的多維全局最優(yōu)解[Δτ=[Δτ1，Δτ2，…，][Δτi…，Δτp]T]，可以表示如下：

[Δτi=limρ→∞∫J…∫JΔτiexp {ρLc（Δτ）}dΔτ∫J…∫Jexp{ρLc（Δu）}dΔu] （13）

式中J是待估計(jì)時(shí)延[Δτ]被限定的取值區(qū)間。

定義一個(gè)偽概率密度函數(shù)：

[L‘c，ρ0Δτ=exp {ρ0LcΔτ}∫J…∫Jexp{ρ0Lc（Δu）}dΔu] （14）

那么，[Δτi]的最優(yōu)值為：

[Δτi=∫J…∫JΔτiL‘c，ρ0ΔτdΔτ] （15）

當(dāng)[ρ0]取值趨近無(wú)窮時(shí)，[L‘c，ρ0Δτ]就變?yōu)榱艘粋€(gè)P維狄拉克（Dirac）函數(shù)。它的極大值所在位置會(huì)趨近[LcΔτ]的極大值所在位置，即使[Δτ]并不是真正的隨機(jī)變量，[L‘c，ρ0Δτ]包含了概率密度函數(shù)的所有性質(zhì)。 然而式（15）需要計(jì)算多維積分，直接進(jìn)行求解存在很大困難，不難想到[Δτi]是[Δτi]的估計(jì)值，是依據(jù)[L‘c，ρ0Δτ分布的Δτ]當(dāng)中的第i個(gè)元素。那么，就可以用蒙特卡洛（Monte?Carlo）方法有效估計(jì)出[Δτ]的期望值[Δτ]，可以表示如下：

[Δτ=1Rk=1RΔτk] （16）

[{Δτk}Rk=1]是根據(jù)[L‘c，ρ0Δτ]生成向量[Δτ]的R次實(shí)現(xiàn)。這樣就可以把式（15）中的多重積分用一個(gè)簡(jiǎn)單的樣本平均替換。此時(shí)新的問(wèn)題又出現(xiàn)了，[L‘c，ρ0Δτ]相對(duì)于待估計(jì)參數(shù)[Δτ]是非線性的，并不適用于簡(jiǎn)單的生成[Δτ]。這時(shí)需要構(gòu)造與[L‘c，ρ0Δτ]相似并且簡(jiǎn)單易于生成實(shí)現(xiàn)[Δτ]的函數(shù)。重要抽樣方法能夠有力地解決多重積分的問(wèn)題，觀察如下的關(guān)系式：

[J…JfΔτL‘c，ρ0ΔτdΔτ=J…JfΔτL‘c，ρ0Δτg′（Δτ）g′（Δτ）dΔτ ] （17）

式中：[g′（Δτ）]是另一個(gè)偽概率密度函數(shù)，稱為歸一化重要性函數(shù)（IF）；[f·]是一個(gè)給定的函數(shù)。

如此可得式（17）左邊就是[fΔτL‘c，ρ0Δτg'（Δτ）]的均值，其中[Δτ]根據(jù)[g′（Δτ）]生成。因此，在構(gòu)造[g′（Δτ）]的時(shí)候，既要保持與[L‘c，ρ0Δτ]的相似性，又要保證易于生成[Δτ]。根據(jù)蒙特卡洛（Monte?Carlo）方法，可以將式（17）寫成如下形式：

[∫J…∫JfΔτL‘c，ρ0Δτg′Δτg′ΔτdΔτ=1Rk=1RfΔτkL‘c，ρ0Δτkg′（Δτk）] （18）

其中[Δτk]是根據(jù)[g′Δτ]對(duì)[Δτ]的第k次實(shí)現(xiàn)。

由式（12）可以看出，[LcΔτ]主要由[（ΦpΔτHΦpΔτ）-1]構(gòu)成，為了避免矩陣的轉(zhuǎn)置，可以用對(duì)角陣[k=1NY1k2-1Ip]來(lái)近似表示[（ΦpΔτHΦpΔτ）-1]。[Ip]是[p×p]確定矩陣。這樣IF函數(shù)[gρ1（Δτ）]可以定義為：

[gρ1Δτ=expρ1σ2Wk=1NY1k2YHΦpΔτΦpΔτHY]

[=i=1Pexp ρ1σ2Wk=1NY1k2I（Δτi）] （19）

式中：[ρ1]是常數(shù)；[I（Δτi）]是在頻域求[Δτi]值的數(shù)據(jù)周期圖，如下所示：

[IΔτi=k=1NY1k*Y（k）exp j2π（k-1）ΔτiN2] （20）

將[gρ1Δτ]歸一化處理得到歸一化重要性函數(shù)（IF）：

[gρ1′′Δτ=i=1Pexp {ρ1′I（Δτi）}∫J…∫Ji=1Pexp {ρ1′I（Δui）}dΔui] （21）

[ρ1′=ρ1σ2Wk=1NY1k2] （22）

理論上歸一化的IF函數(shù)[gρ1′′Δτ]的圖像含有P個(gè)峰，每個(gè)峰出現(xiàn)在不同時(shí)延所對(duì)應(yīng)的位置。因此每實(shí)現(xiàn)一次[gρ1′′Δτ]，就會(huì)生成一組[Δτk]，[Δτk（i）]分別在每一個(gè)峰所對(duì)應(yīng)的位置附近生成。然而由于噪聲影響，會(huì)出現(xiàn)一些次級(jí)波峰對(duì)生成時(shí)延造成影響。增大參數(shù)[ρ1′]可以減少次級(jí)波峰，但是同時(shí)也可能毀壞有用的波峰。最優(yōu)的[ρ1′]是當(dāng)IF函數(shù)剛好出現(xiàn)P個(gè)波峰的時(shí)候，所對(duì)應(yīng)的最大的[ρ1′]值。

式（21）呈現(xiàn)出很多優(yōu)勢(shì)，比如可以將每個(gè)時(shí)延在[gρ1Δτ]中的聯(lián)合貢獻(xiàn)分離成了每個(gè)時(shí)延的個(gè)體貢獻(xiàn)的乘積形式[7]。因此，通過(guò)[gρ1Δτ]對(duì)[Δτ]進(jìn)行實(shí)現(xiàn)也就可以轉(zhuǎn)換成對(duì)[Δτ]中的每一個(gè)[Δτi]的實(shí)現(xiàn)，同時(shí)也就將求[Δτ]的問(wèn)題轉(zhuǎn)化成求[Δτi]的問(wèn)題。

根據(jù)（18）可以得到：

[Δτi=1Rk=1RΔτk（i）FΔτk] （23）

[FΔτk=L‘c，ρ0Δτkgρ1′′Δτk] （24）

根據(jù)式（23）可以估計(jì)出[Δτ=[Δτ1，Δτ2，…，Δτi…，Δτp]T]。

[FΔτk]作為一個(gè)加權(quán)系數(shù)，盡可能減小偏離[Δτi]的實(shí)現(xiàn)對(duì)估計(jì)[Δτi]時(shí)產(chǎn)生的影響。為了避免指數(shù)運(yùn)算可能造成的計(jì)算數(shù)值溢出，可以將[FΔτk]用[F′Δτk]進(jìn)行替換[8]：

[F′Δτk=expρ0LcΔτk-ρ1′i=1PIiΔτki- max1≤l≤R（ρ0LcΔτl-ρ1′i=1PIiΔτli）] （25）

由于待求參數(shù)[Δτi]具有上界和下界，這里用三角均值的概念來(lái)求解[Δτi]比較合適[9]，因此，可以表述成如下的形式：

[Δτi=12πL∠1Rk=1RF（Δτk）exp j2πΔτk（i）L] （26）

式中L是待估計(jì)時(shí)延可能的取值區(qū)間長(zhǎng)度。

3 仿真和討論

這里采用計(jì)算機(jī)的仿真來(lái)檢驗(yàn)上述的算法，仿真使用Matlab語(yǔ)言。

仿真過(guò)程中，聲源為遠(yuǎn)場(chǎng)聲源，入射角設(shè)置為與法線夾角45°，采用語(yǔ)音信號(hào)。在室溫下，設(shè)置3陣元間距為0.1 m的均勻直線麥克風(fēng)陣列，麥克風(fēng)采樣率為44.1 kHz，信號(hào)幀長(zhǎng)為4 096個(gè)采樣點(diǎn)。第一個(gè)陣元設(shè)為參考陣元，求陣元和參考陣元接收信號(hào)的時(shí)延。為了更客觀地對(duì)重要抽樣的時(shí)延估計(jì)方法進(jìn)行評(píng)估，這里將其與ML的網(wǎng)格搜索方法進(jìn)行比較。

圖1、圖2描繪出了聲源是語(yǔ)音信號(hào)時(shí)，時(shí)延估計(jì)在不同信噪比（SNR）下的均方誤差（MSE）。觀察實(shí)驗(yàn)結(jié)果，重要抽樣的方法的性能可以趨近ML網(wǎng)格搜索法，當(dāng)信噪比大于0時(shí)，可以估計(jì)出較為準(zhǔn)確的時(shí)延，在高信噪比時(shí)，性能也處于同一數(shù)量級(jí)，但是在計(jì)算復(fù)雜度上比ML網(wǎng)格搜索法低一個(gè)數(shù)量級(jí)。

圖1 估計(jì)第一個(gè)時(shí)延的均方誤差（MSE）

圖2 估計(jì)第二個(gè)時(shí)延的均方誤差（MSE）

對(duì)算法的復(fù)雜度進(jìn)行一個(gè)簡(jiǎn)要的分析。若需要估計(jì)P個(gè)時(shí)延，每個(gè)時(shí)延進(jìn)行N個(gè)點(diǎn)的搜索，對(duì)[LcΔτ]每1次搜索要進(jìn)行5次矩陣乘法運(yùn)算、2次矩陣共軛轉(zhuǎn)置、1次矩陣求逆。對(duì)[gρ1Δτ]每1次搜索要進(jìn)行2次矩陣共軛轉(zhuǎn)置、3次矩陣乘法。那么，網(wǎng)格搜索法就對(duì)[LcΔτ]進(jìn)行了[NP]次搜索，重要抽樣法對(duì)[gρ1Δτ]進(jìn)行了[P×N×R]次搜索，R為對(duì)[Δτ]的實(shí)現(xiàn)生成次數(shù)。隨著P的增加，R也要相應(yīng)增大才能保證時(shí)延估計(jì)得性能不下降。設(shè)N=100，P取1～4，R對(duì)應(yīng)取10，50，100，250。圖3描繪出了估計(jì)時(shí)延個(gè)數(shù)不同時(shí)，兩種算法的大致計(jì)算復(fù)雜度對(duì)數(shù)值。顯然隨著估計(jì)時(shí)延個(gè)數(shù)的增加，ML網(wǎng)格搜索法計(jì)算復(fù)雜度急劇上升，而重要抽樣法計(jì)算復(fù)雜度緩慢上升。因此，重要抽樣法更適合于麥克風(fēng)陣列時(shí)延的實(shí)時(shí)估計(jì)。

圖3 計(jì)算復(fù)雜度對(duì)比

4 結(jié) 語(yǔ)

將重要抽樣的思想應(yīng)用到麥克風(fēng)陣列的時(shí)延估計(jì)當(dāng)中，是對(duì)傳統(tǒng)的極大似然估計(jì)ML方法的一個(gè)改進(jìn)，有效避免了多維網(wǎng)格搜索和迭代運(yùn)算所造成的計(jì)算復(fù)雜度以及初始值的問(wèn)題。相比于傳統(tǒng)的方法，在保證時(shí)延估計(jì)準(zhǔn)確的前提下，有效降低了計(jì)算復(fù)雜度，更適合于時(shí)延的實(shí)時(shí)估計(jì)。

參考文獻(xiàn)

[1] 劉建平，張一聞，劉穎.基于麥克風(fēng)陣列的汽車笛語(yǔ)識(shí)別及笛聲定位方法[J].西安電子科技大學(xué)學(xué)報(bào)，2012，39（1）：163?167.

[2] HAHN P J， MATHEWS V J， TRAN T D. Adaptive realization of a maximum likelihood time delay estimator [C]// the 1996 IEEE International Conference on Acoustics， Speech， and Signal Processing. Atlanta， GA： IEEE， 1996： 3121?3124.

[3] 陶巍，劉建平，張一聞.基于麥克風(fēng)陣列的聲源定位系統(tǒng)[J].計(jì)算機(jī)應(yīng)用，2012，32（5）：1457?1459.

[4] NG L， BAR SHALOM Y. Multisensor multitarget time delay vector estimation [J]. IEEE Transactions on Acoustics， Speech and Signal Processing， 1986， 34（4）： 669?677.

[5] PINCUS M. A closed form solution for certain programming problems [J]. Operations Research， 1968， 16（3）： 690?694.

（上接第39頁(yè)）

[6] PALLAS M， JOURDAIN G. Active high resolution time delay estimation for large BT signals [J]. IEEE Transactions on Signal Processing， 1991， 39（4）： 781?787.

[7] SAHA S， KAY S. An exact maximum likelihood narrowband direction of arrival estimator [J]. IEEE Transactions on Signal Processing， 2008， 56（4）： 1082?1092.

[8] MCKILLIAM R G， QUINN B G. Estimating the circular mean from correlated observations [M]. [S.l.]： Daspworkshop Org， 2011.