基于SRP-PHAT 的實(shí)時聲源定位算法設(shè)計與實(shí)現(xiàn)

劉 生

（南京熊貓電子裝備有限公司，江蘇南京 210046）

0 引言

近年來，隨著人工智能及交互技術(shù)的發(fā)展，越來越多智能設(shè)備可以實(shí)現(xiàn)人機(jī)對話。為了提高自動語音識別率，業(yè)內(nèi)通常使用波束形成技術(shù)進(jìn)行定向收音，獲得了較好的收音效果［1］。針對一些如智能音箱、機(jī)器人等智能對話設(shè)備，為了實(shí)現(xiàn)全向的高效收音，通常會在對話開始時利用聲源定位技術(shù)獲得聲音的方向，并將其作為參數(shù)傳給后續(xù)的波束形成部分進(jìn)行處理。

常用的聲源定位算法主要有基于最大可控響應(yīng)功率的波束形成方法、基于高分辨率譜估計的方法與基于時延估計的定位方法［2-4］?；谙辔蛔儞Q加權(quán)的可控波束響應(yīng)功率（Steered Response Power-Phase Transform，SRP-PHAT）聲源定位算法由于具有較強(qiáng)的魯棒性，可實(shí)現(xiàn)真實(shí)環(huán)境中的聲源定位，因此也得到了廣泛應(yīng)用［5-6］。由于算法運(yùn)算量較大，在嵌入式系統(tǒng)中，為提高算法的實(shí)時性，一般使用FPGA 或?qū)Ｓ肅MOS 芯片進(jìn)行算法實(shí)現(xiàn)［7-8］。在智能語義交互產(chǎn)品中，目前的主流技術(shù)是使用Linux 或Android 系統(tǒng)進(jìn)行網(wǎng)絡(luò)通信，連接語音對話服務(wù)器，實(shí)現(xiàn)語義識別、語義處理及語音合成。這部分功能需要使用通用ARM 處理器進(jìn)行實(shí)現(xiàn)，如果使用FPGA 或?qū)Ｓ肅MOS 芯片實(shí)現(xiàn)語音定位及前端語音處理部分功能，勢必會提高產(chǎn)品復(fù)雜度，成本也會隨之升高。本文基于ARM 平臺處理器，使用Python 及C++語言對算法進(jìn)行設(shè)計與實(shí)現(xiàn)。將原來利用專用芯片實(shí)現(xiàn)的前端語音功能在嵌入式ARM 平臺上重新進(jìn)行設(shè)計與實(shí)現(xiàn)，對產(chǎn)品集成度進(jìn)行了提升。

1 SPR-PHAT 算法原理

對空間中特定位置或方向信號進(jìn)行定向增強(qiáng)的技術(shù)是陣列信號處理中的一項關(guān)鍵技術(shù)。這些技術(shù)大多基于波束形成理論，主要通過麥克風(fēng)陣列捕捉聲音并進(jìn)行相關(guān)信號處理。傳統(tǒng)意義上的波束形成可通過濾波求和處理，在對麥克風(fēng)信號求和以產(chǎn)生定向信號聚焦之前，會在時域?qū)溈孙L(fēng)信號進(jìn)行濾波。濾波器通常會在波束形成過程中進(jìn)行動態(tài)調(diào)整，以增強(qiáng)所需的信號源，同時衰減其他信號源。一種最簡單的對信號進(jìn)行濾波處理的方法就是對源信號進(jìn)行時間偏移，該偏移值為各個陣列間的信號傳播延遲，這種方法被稱為延遲求和（Delay-and-Sum）波束形成［9-10］。通過延遲麥克風(fēng)信號，將各通道的源信號在時間上對準(zhǔn)，之后對所有通道的源信號進(jìn)行求和。應(yīng)用更復(fù)雜濾波求和技術(shù)的濾波器以及其他信號增強(qiáng)處理算法通常都是基于這種時間對準(zhǔn)方法。

波束形成技術(shù)在麥克風(fēng)陣列中得到了廣泛應(yīng)用［11-13］。對于這種應(yīng)用，濾波求和技術(shù)不僅需要抑制背景噪聲以及不需要的信號，同時還必須保證想要的信號通過濾波器后不能失真。當(dāng)將波束形成技術(shù)應(yīng)用于聲源定位時，這些濾波器只需在完成波束形成運(yùn)算后，對特定方向的信號進(jìn)行功率增強(qiáng)即可。本文使用的PHAT 算法正是利用了該特性，使用可控波束響應(yīng)功率算法對所有可能的相位變換進(jìn)行求和運(yùn)算。SRP-PHAT 可直接對多路麥克風(fēng)信號進(jìn)行變換處理，利用多個麥克風(fēng)提高位置估計的準(zhǔn)確性。

SRP 可以使用塊處理方案加以實(shí)現(xiàn)，該方案使用短時數(shù)字傅里葉變換作為麥克風(fēng)信號頻譜的估計［14-15］。將陣列信號在時域上分割成多個塊，并計算每個塊的轉(zhuǎn)向響應(yīng)。信號塊的數(shù)字傅里葉用Xm,b[k]表示，其中m 是麥克風(fēng)索引，b 是塊索引，Gm,b[k]是麥克風(fēng)m 離散時間濾波器的傅里葉變換，會在每個塊b 中單獨(dú)進(jìn)行計算。塊b 的轉(zhuǎn)向響應(yīng)可定義如下：

[k,Δ1...ΔM)是索引為k 的離散頻率函數(shù)和連續(xù)的轉(zhuǎn)向延遲（Δ1...ΔM表示陣列所有連續(xù)的轉(zhuǎn)向延遲）。理論上需要對信號中所有頻率段的數(shù)據(jù)進(jìn)行處理，但在實(shí)際使用過程中，一般會選擇一個或多個頻率段的數(shù)據(jù)進(jìn)行處理。同時，盡管M 個轉(zhuǎn)向延遲是連續(xù)的，但在實(shí)際使用過程中，會在預(yù)定義的一組空間位置或方向進(jìn)行采樣，并通過對K 個離散頻率進(jìn)行求和，從而獲得轉(zhuǎn)向響應(yīng)功率：

離散濾波器如下：

將式（3）代入式（1），相位變換加權(quán)的可控響應(yīng)表示為：

將式（4）代入式（2），可獲得相位變換加權(quán)的可控響應(yīng)功率SRP-PHAT 為：

2 算法實(shí)現(xiàn)

為方便計算，首先對多個通道的麥克風(fēng)數(shù)據(jù)進(jìn)行DFT變換，根據(jù)麥克風(fēng)的物理參數(shù)計算在某特定方向上的可控時延；然后基于音頻數(shù)據(jù)的DFT 以及可控時延，結(jié)合公式（5），可求得該方向所有頻率的SRP-PHAT；最后重復(fù)上述運(yùn)算，直到獲得所有方向的SRP-PHAT 值，并選取最大值對應(yīng)的方向作為音源方向。具體流程如圖1 所示。

Fig.1 Flow of SRP-PHAT algorithm implementation圖1 SRP-PHAT 算法實(shí)現(xiàn)流程

上述程序可利用Python 的NumPy 庫方便地實(shí)現(xiàn)如矩陣內(nèi)積、點(diǎn)積以及共軛轉(zhuǎn)置等諸多矩陣運(yùn)算，并可使用NumPy 庫中的fft 函數(shù)實(shí)現(xiàn)傅里葉變換［16］。但如果要在嵌入式平臺上運(yùn)行，考慮到部署方便的問題，實(shí)驗一般使用C++實(shí)現(xiàn)。由于在算法中需要使用傅里葉變換、線性代數(shù)等運(yùn)算，必須使用相應(yīng)的數(shù)學(xué)庫。本文使用了Xtensor、Xtl、OpenBLAS 等數(shù)學(xué)運(yùn)算庫［17-19］。Xtensor 是一個C++庫，用于多維數(shù)組表達(dá)式的數(shù)值計算，可提供支持延遲廣播的可擴(kuò)展表達(dá)式系統(tǒng)，Xtensor 提供的API 遵循C++標(biāo)準(zhǔn)。

基于Xtensor 庫有很多擴(kuò)展包，其中Xtensor-fftw 主要實(shí)現(xiàn)了xt：：fftw 函數(shù)以進(jìn)行fft 運(yùn)算。fftw 是一個用來計算一維或多維離散傅里葉變換（DFT）的C 語言函數(shù)庫，其支持輸入任意大小的實(shí)數(shù)及復(fù)數(shù)運(yùn)算。Xtensor-blas 是對Xtensor 庫的擴(kuò)展，基于FLENS 項目的CXXBLAS 和cxxla?pack 將BLAS 與LAPPACK 庫進(jìn)行綁定，具有計算逆矩陣、求特征值、解線性方程組及求解行列式等功能。其提供的函數(shù)接口基本可直接與NumPy 庫兼容，因而方便程序移植。上述幾個庫還依賴Xtl（提供了基本的科學(xué)計算函數(shù)庫）和OpenBLAS（Basic Linear Algebra Subprograms 基礎(chǔ)線性代數(shù)子程序庫），這些庫之間的關(guān)系如圖2 所示。

Fig.2 Algorithm library structure圖2 算法庫結(jié)構(gòu)

3 實(shí)驗分析

3.1 實(shí)驗環(huán)境

在性能對比實(shí)驗中，主要對算法的準(zhǔn)確性及實(shí)時性進(jìn)行分析。通過選取編程語言、運(yùn)行環(huán)境、頻率范圍、采樣時間等幾種不同變量進(jìn)行多組組合實(shí)驗，以獲得這些參數(shù)對算法準(zhǔn)確性及性能的影響。實(shí)驗采用的錄音設(shè)備為Re?Speaker Core v2.0［20］，該麥克風(fēng)陣列排列以及實(shí)驗時收音方向如圖3 所示。

Fig.3 Microphone array arrangement and sound reception diagram圖3 麥克風(fēng)陣列排列及收音示意圖

3.2 準(zhǔn)確性實(shí)驗

準(zhǔn)確性實(shí)驗主要對算法定位的準(zhǔn)確性進(jìn)行測試。定義中心到MIC 的方向為0°，中心到MIC4 的方向為180°，錄制在180°方向、距離約50cm 位置處人聲閱讀文章的音頻。經(jīng)6 麥克風(fēng)錄制后的音頻數(shù)據(jù)格式為6 通道，采樣深度為16bit，采樣頻率為16K。具體音頻如圖4 所示，音頻時間長度為1min。

Fig.4 6-channel audio of microphone array圖4 麥克風(fēng)陣列6 聲道收音音頻

在實(shí)際應(yīng)用中，頻段范圍的選擇比較重要，由于人聲的頻率范圍集中在300～3 000Hz，選取該頻率范圍可在盡量保證準(zhǔn)確性的前提下，減小算法運(yùn)算量。為考察頻率范圍對算法準(zhǔn)確性的影響，本實(shí)驗選取0～8 000Hz 與300～3 000Hz 范圍進(jìn)行對比實(shí)驗。具體方法為將1min 時間以50ms 進(jìn)行切片，分別進(jìn)行處理，共獲得1 200 個樣本。對處理結(jié)果進(jìn)行統(tǒng)計，獲得不同頻段范圍的準(zhǔn)確率如圖5 所示。

從實(shí)驗結(jié)果可以看出，在0～8 000Hz 及300～3 000Hz頻段范圍進(jìn)行處理時，定位偏差在±5°范圍內(nèi)的概率分別為92%及81%。0～8 000Hz 頻段的誤差范圍可控制在±14°以內(nèi)，而300～3 000Hz 頻段的誤差范圍為±18°以內(nèi)。可以看出，在采用了較小頻率范圍后，整體精度還是降低了10%～20%。

為了得到不同時間片對處理準(zhǔn)確性的影響，分別以25ms、50ms 及500ms 對數(shù)據(jù)進(jìn)行切片，然后對切片數(shù)據(jù)分別進(jìn)行處理，并對處理結(jié)果的偏差值進(jìn)行統(tǒng)計。結(jié)果如圖6 所示。

Fig.6 Accuracy of different time slices圖6 不同時間片的準(zhǔn)確率

從實(shí)驗結(jié)果來看，25ms 的數(shù)據(jù)切片在0°方向的準(zhǔn)確率最高，但在±4°以后的準(zhǔn)確率小于其他兩種切片方式，整體準(zhǔn)確率基本在控制±17°以內(nèi)。從整體來看，不同時間切片對聲源定位的準(zhǔn)確性并無明顯影響。

3.3 實(shí)時性實(shí)驗

很多嵌入式設(shè)備對算法的實(shí)時性都有一定要求，因此本節(jié)對該算法的實(shí)時性能進(jìn)行測試與分析。選取1 000ms的6 麥克風(fēng)數(shù)據(jù)進(jìn)行處理，相位變換加權(quán)角度的精度分別為1°、2°與4°。選取X86 平臺以及嵌入式ARM 平臺，X86處理器為i7-4790，嵌入式ARM 處理器為RK3399，運(yùn)行的操作系統(tǒng)均為Linux Ubuntu 發(fā)行版。0～8 000Hz 及300～3 000Hz 頻段在各平臺運(yùn)行耗時如圖7 所示。

Fig.7 Time consuming of each platform圖7 各平臺運(yùn)行耗時

從實(shí)驗結(jié)果來看，相對于0～8 000Hz 頻段范圍，采用300～3 000Hz 頻段范圍能夠減少約40%的處理時間。

業(yè)內(nèi)認(rèn)為一般情況下，C++語言運(yùn)行效率高于Python。為了對比使用Python 和C++的算法性能區(qū)別，對應(yīng)用兩種語言在X86 平臺的處理速度進(jìn)行測試?？煽吹酵ㄟ^兩種不同語言實(shí)現(xiàn)的程序運(yùn)行時的性能差距較小，這是因為在使用Python 語言實(shí)現(xiàn)算法時，已利用cPython 庫加速矩陣預(yù)算過程，因此在使用C++語言后，性能并未得到較大提升。

最后，使用C++語言實(shí)現(xiàn)該算法，并在嵌入式ARM 處理器RK3399 芯片上進(jìn)行性能測試。ARM 平臺由于處理能力有限，運(yùn)行耗時約為X86 平臺的2.8 倍。

在一些產(chǎn)品交互應(yīng)用，特別是語音機(jī)器人中，主要是在喚醒操作時對聲源進(jìn)行定位，提供聲源方向，為后續(xù)波束形成提供音源方向參數(shù)。根據(jù)實(shí)驗結(jié)果，在嵌入式平臺通過本算法對1 000ms 的6 麥克風(fēng)數(shù)據(jù)進(jìn)行處理需要110ms 的時間，在人機(jī)對話交互過程中，該時間在可接受范圍內(nèi)。如果希望進(jìn)一步提高算法響應(yīng)速度，可根據(jù)需求減小相位變換角度分辨率。繼續(xù)在ARM 處理器上進(jìn)行實(shí)驗，分別選取角度分辨率為1°、2°與4°，運(yùn)行多次取均值，實(shí)驗結(jié)果如圖8 所示。

根據(jù)實(shí)驗結(jié)果，可看出通過降低角度分辨率，可以線性地降低算法復(fù)雜度，減少運(yùn)行時間，但需要在實(shí)現(xiàn)過程中通過綜合考慮進(jìn)行效率及精度取舍。

對于音源持續(xù)時間與處理耗時的關(guān)系也進(jìn)行對比實(shí)驗，選取音頻時間為1 000ms、500ms 和50ms，在ARM 平臺上運(yùn)行多次取均值，結(jié)果如圖9 所示。

Fig.9 Time consuming of different audio duration圖9 不同音頻時長運(yùn)行耗時

可以看到，對于不同時長的音頻進(jìn)行處理，其耗時并未隨音頻時長的減少而大幅下降，減少的時間大部分來自對源數(shù)據(jù)進(jìn)行預(yù)處理及傅里葉變換導(dǎo)致的時間變化。因為傅里葉變換后將時域數(shù)據(jù)轉(zhuǎn)為頻域信息，對頻域數(shù)據(jù)進(jìn)行相位變換加權(quán)后的數(shù)據(jù)矩陣大小與時長無關(guān)，比較耗時的部分是在角度搜索中進(jìn)行互功率最大值運(yùn)算。

4 結(jié)語

本文簡要分析了SRP-PHAT 算法原理，并提出使用Py?thon 及C++語言實(shí)現(xiàn)該聲源定位算法。針對在工程實(shí)踐中的角度分辨率及音源時長等相關(guān)參數(shù)對算法準(zhǔn)確性與實(shí)時性的影響進(jìn)行實(shí)驗與分析，為算法部署時的參數(shù)選取提供了依據(jù)。實(shí)驗結(jié)果表明，該算法在嵌入式平臺處理耗時在200ms 以內(nèi)，能基本滿足嵌入式平臺的應(yīng)用要求，也可根據(jù)需求選取不同的音頻頻率段，調(diào)整搜索精度以進(jìn)一步提高算法的實(shí)時性能。根據(jù)實(shí)驗結(jié)果，該方案可應(yīng)用于智能語音交互產(chǎn)品中，使用單個ARM 芯片即可實(shí)現(xiàn)聲源定位，并通過與服務(wù)器的通信實(shí)現(xiàn)語音識別、語義分析等關(guān)鍵功能。下一步在產(chǎn)品具體實(shí)施過程中，需要結(jié)合噪聲抑制以及回聲消除等其他音頻處理技術(shù)，進(jìn)一步提高聲源定位的準(zhǔn)確性。