當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

基于时延法的麦克风阵列声源定位分析

發布時間：2023/12/8 编程问答 54 豆豆

生活随笔收集整理的這篇文章主要介紹了基于时延法的麦克风阵列声源定位分析小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

文章目錄

一. 關于麥克風陣列
二. 關于聲源定位
三. 基于廣義互相關(GCC)計算時延
四. 基于時延差的聲源定位法
- 1. 近場模型
- 2. 遠場模型
五. 三維空間陣列的聲源定位系統實現
- 1. 推導過程
六. 六元圓形麥克風陣列聲源定位
七. 相關鏈接

一. 關于麥克風陣列

麥克風陣列： 麥克風陣列是由一定數目的聲學傳感器（麥克風）按照一定規則排列的多麥克風系統，而基于麥克風陣列的聲源定位是指用麥克風拾取聲音信號，通過對麥克風陣列的各路輸出信號進行分析和處理，得到一個或者多個聲源的位置信息。

麥克風陣列系統的聲源定位技術研究意義在于： 輸入的信息只有兩個方向難以確定聲源的位置，人類的聽覺系統主要取決于頭和外耳氣壓差聲波實現聲源定位。假使沒有這個壓力差，只能定位在平面上聲源的位置，但就無法知道聲音是從前面，或從后面傳來的。因此，由人的聽覺系統，科技研發人員得到了靈感，使用多個麥克風系統可以實現在三維空間中的聲源位置的定位，麥克風的數量越多，所接收到的信息量也越多。聲源的聲源定位和聲源增強是實現智能處理的兩個關鍵問題，而聲源定位是實現語音增強的前提和基礎。一個麥克風的信息量較少，使得聲源定位所需的信息缺乏，而麥克風陣列克服了上述缺點，充分利用每個麥克風信號之間的數據相關性，并加以融合，可以實現聲源定位。

麥克風陣列聲源定位技術的應用： 廣泛應用于國防、智能機器人、視頻會議及語音增強等眾多領域,尤其在當下以智能辦公和智能家居為主要室內場景的遠場語音交互系統中。

二. 關于聲源定位

目前基于麥克風陣列的聲源定位方法主要有三種：基于最大輸出功率的可控波束成形的定位方法、基于高分辨譜估計的定位方法、基于到達時延差估計的定位方法(Time Difference of Arrival，TDOA)。

基于最大輸出功率的可控波束成形定位法。 波束形成法的原理是將麥克風接收到的信號進行濾波加權求和來形成波束，按照一定的規律對聲源位置進行搜索，當麥克風達到最大輸出功率時，為時搜索到的聲源位置即為真實的聲源方位。波束形成可分為常規的波束形成CBF（Conventional Beam Forming）、CBF+Adaptive Filter和自適應波束形成ABF（Adaptive Beam Forming）。

基于高分辨譜估計的定位法。 基于高分辨率譜估計的定位方法通過分解協方差矩陣估計聲源方位。適合多個聲源的情況，且聲源的分辨率與陣列尺寸無關，突破了物理限制。該方法的優點是不受采樣頻率限制，且在一定程度下可以實現任意程度的定位，但是該方法計算復雜度較高，抗噪和抗混響性能較差，因此該方法適合在一些特定的環境下使用。這類方法可以拓展到寬帶處理，但是對誤差十分敏感（如麥克風單體誤差，通道誤差），適合遠場模型，且矩陣運算量巨大。

基于到達時延差估計的定位法。 TDOA（time difference of arrival）是先后估計聲源到達不同麥克風的時延差，通過時延來計算距離差，再利用距離差和麥克風陣列的空間幾何位置來確定聲源的位置。可分為TDOA估計（估計信號到達各麥克風的時間差）和TDOA定位（運用幾何關系確定聲源位置）兩步。

三. 基于廣義互相關(GCC)計算時延

時延估計有很多種，比較經典就是廣義互相關函數 (Generalized Cross Correlation, GCC) 估計時延，這里簡單介紹基于廣義互相關函數估計時延的方法。

在噪聲存在情況下，一個由遠處聲源發出的，并且被兩個不同空間中的麥克風監聽的信號可以數學建模為：

x_1=s_1(t)+n_1(t)

?????? ???????????????????(1)

x_2=αs_1(t-D)+n_2(t)

???????? ??????(2)

其中， $s (t)$ 是聲音信號， $n_1(t)、n_2(t)$ 是兩個聲音傳感器檢測噪聲。三者是穩定的隨機過程，且互不相關。

計算 $x_1$ 與 $x_2$ 的互相關函數：

R_{x_1x_2}(τ)=E [ x_1(t)*x_2(t-τ) ]

???? ? ? ? ??? (3)

R^x1x2(τ)=1T?τ∫τTx1(t)x2(t?τ)dt\hat R_{x_1x_2}(τ)=\frac{1}{T-τ}\int_τ^T {x_1(t)x_2(t-τ)} \,{\rm d}t

???? ? (4)

其中估計的時延 $D$ 為互相關函數值達到最大值時取得的 $τ$ 值，即：

D^=argmaxτRx1x2(τ)\hat D=\underset {\ τ}{argmax} R_{x_1x_2}(τ)

???? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (5)

MATLAB 例程：

% 導入兩個麥克風的音頻數據 [y_0,Fs] = audioread('音軌-0.wav'); [y_1] = audioread('音軌-1.wav');fprintf('采樣頻率：%d\n ······\n', Fs);% 取出兩段音頻前2048個采樣點，并作互相關處理，繪制曲線。 A = y_0(1:2048); B = y_1(1:2048);[value,delay] = xcorr(A,B);subplot(1,2,1); plot(delay, value);D = zeros(1,926); % 以2048個點為一幀，計算互相關得到的時延，繪制出兩段音頻的時延變化。 for a = 1:2048:size(y_0,1)-2048A = y_0(a:a+2048);B = y_1(a:a+2048);[value,delay]=xcorr(A,B);value_max_idx = find(value==max(value));D1 = delay(value_max_idx);D((a-1)/2048+1) = D1; endsubplot(1,2,2); plot(D);

使用ReSpeaker的樹莓派六麥克風套件，拿手機放歌作為聲源圍繞麥克風陣列移動，錄制音頻，并保存為6個wav音頻文件，采用率為44100Hz。對其中兩個麥克風音頻文件做處理，如果取初始的2048個采樣點做互相關處理就能得到下面圖一，圖一的峰值的橫坐標即為估計時延。如果以2048個采樣點作為一幀，對兩段音頻的每一幀進行互相關，并找到每幀的估計時延，繪制出來就能得到時延的變化曲線，如下面圖二。

注意： 這里的時延并不是指具體的時間，而是采樣點數。想要得到具體時間還需要除以采樣頻率。以2048采樣點為一幀，就是以0.464s為一幀（ $204844100≈0.464\frac{2048}{44100}≈0.464$ s）。

四. 基于時延差的聲源定位法

在對麥克風陣列進行建模之前，我們需要分清楚什么近場與遠場。顧名思義，離麥克風近則符合近場模型 ，離得遠則符合遠場模型 。

假設 $L$ 為陣列間距， $λ$ 為聲波波長， $M$ 為聲源與麥克風的距離，我們定義 $2L2λ\frac{2L^2}{λ}$ 為遠近場臨界值。

當 $M<2L2λM<\frac{2L^2}{λ}$ 時，符合近場模型，此時聲源到達麥克風陣列的波形視為球面波。

當 $M<2L2λM<\frac{2L^2}{λ}$ 時，符合遠場模型，此時聲源到達麥克風陣列的波形視為平面波。

可聽聲的機械波頻帶為20Hz ~20000Hz，機械波波長大約在1.7cm ~ 17m（聲速取340m/s）。然而在現實生活中，過高頻率或過低頻率的聲波都是非常少量的，以人的聲音為例，人語音頻帶的范圍大概為300Hz至3400Hz，波長范圍為0.1m~ 1.12m，當陣列間隔取4cm時，遠近場臨界值范圍為 0m~3.2m。若取中間值，則可以認為1.6m內符合近場模型，1.6m外符合遠場模型。

1. 近場模型

當 $M<2L2λM<\frac{2L^2}{λ}$ 時，符合 近場模型，此時聲源到達麥克風陣列的波形視為 球面波。

近場模型至少需要3個麥克風，以最簡單的3麥克風模型為例（如圖： $y_1、y_2、y_3$ ）。假設 $τ_{12} 、τ_{13}$ 分別表示第一個麥克風與第二和第三個麥克風之間的時延，那么有：

τ12=r2?r1cτ_{12}=\frac{r_2-r_1}{c}

?????? ??(6)

τ13=r3?r1cτ_{13}=\frac{r_3-r_1}{c}

?????? ??(7)

其中， $c$ 為聲速。標準大氣壓、 $15°$ 條件下，聲速為 $340 m / s$ 。

根據麥克風陣列的的幾何關系，由余弦定理，可以得到：

r_2^2=r_1^2+d^2+2r_1dcosθ_1

?????(8)

r_3^2=r_1^2+4d^2+4r_1dcosθ_1

????(9)

其中， $τ_{12}、τ_{13}$ 可以通過互相關GCC求得（下面會細說），且 $c 、 d$ 為已知。結合(6) ~ (9) 即可求出未知量 $r_1、r_2、r_3、θ$ ，結合坐標系可以求出 $s (k)$ 的坐標。

2. 遠場模型

當 $M<2L2λM<\frac{2L^2}{λ}$ 時，符合 遠場模型，此時聲源到達麥克風陣列的波形視為 平面波。

對于兩個麥克風的情況，只能計算出方位角，無法計算出方位距離。假設 $τ$ 為聲波到達兩個麥克風之間的時延，則有：

τ=dcosθcτ=\frac{dcosθ}{c}

??????????????????????(10)

θ=arccoscτdθ=arccos\frac{cτ}ozvdkddzhkzd

??????????????(11)

五. 三維空間陣列的聲源定位系統實現

計算方法：幾何求解、最小二乘法

上面所舉的近場、遠場的定位模型例子，分析都是最簡模型下，平面聲源的定位。下面以三維麥克風陣列為例，通過幾何的方式求解三維空間下聲源的定位。

1. 推導過程

一個三維麥克風陣列，麥克風分別為 $m_i(i=0、1···n)$ ，聲源 $S$ 符合近場模型。現以麥克風 $m_0$ 為原點，如圖建立直角坐標系。推導公式之前，需要先確定以下概念：

序號概念字符

1	麥克風的坐標	$m_i$ ? $\in[0,n]$
2	聲源的估計坐標	$S$
3	聲源 $s$ 到 $m_i$ 與 $m_1$ 的估計距離差	$d^i\hat d_i$
4	麥克風 $m_i$ 到原點的的距離	$R_i$
5	聲源 $s$ 到原點的的距離	$R_s$

如上圖，根據三角形 $m_0 m_i s$ ，由余弦定理有：

R_s+d_i)^2=R_i^2-2m_i^Ts+R_s^2

??????????????????????(12)

公式 (12) 展開后得到：

R_i^2-2m_i^Ts-2R_sd_i-d_i^2=0

?????????????????????????(13)

由于 $d_i$ 是通過估計的時延得到的，與實際值之間會有一個偏差，因此公式 (10) 不為零，可寫為：

ε=(R_i^2-d_i^2)-2m_i^Ts-2R_sd_i

??????????????????????(14)

此時已經得到目標值的誤差函數，使用最小二乘法求解，問題可以轉化為：估計聲源坐標 $s (x, y, z)$ ，使得最終的誤差平方和最小，即：

argmins(x,y,z)∑i=1n[(Ri2?di2)?2miTs?2Rsdi]2\underset {\ s(x,y,z)}{argmin}\displaystyle\sum_{i=1}^{n} [(R_i^2-d_i^2)-2m_i^Ts-2R_sd_i]^2

??????(15)

將公式 (14) 寫成矩陣形式：

ε=2RsD^+2Ms?δε=2R_s\hat D+2Ms-\delta

??????????????????????????????????????????(16)

其中， $M=[m1m2m3....mn]=[x1y1z1x2y2z2x3y3z3.........xnynzn]M=\begin{bmatrix}m_1\\m_2\\m_3\\....\\m_n\end{bmatrix}=\begin{bmatrix}x_1&y_1&z_1\\x_2&y_2&z_2\\x_3&y_3&z_3\\...&...&...\\x_n&y_n&z_n\end{bmatrix}$ ??? ? $D^=[d^1d^2d^3...d^n]\hat D=\begin{bmatrix}\hat d_1\\\hat d_2\\\hat d_3\\\ ... \\\hat d_n\end{bmatrix}$ ??????? $δ=[R12?d^12R22?d^22R32?d^32..........Rn2?d^n2]\delta=\begin{bmatrix}R_1^2-\hat d_1^2\\R_2^2-\hat d_2^2\\R_3^2-\hat d_3^2\\..........\\R_n^2-\hat d_n^2\end{bmatrix}$

公式 (16) 可以簡化為：

ε=Aμ?bε=A\mu-b

???????????????????????????????????(17)

其中， $A=[MD^]A=\begin{bmatrix}M&\hat D\end{bmatrix}$ ????? $μ=[sRs]\mu=\begin{bmatrix}s\\R_s\end{bmatrix}$ ?????? $b=12δb=\frac{1}{2}\delta$

公式 (17) 的最小二乘解可以表示為：

μ^=(ATA)?1ATb\hat\mu=(A^TA)^{-1}A^Tb

???????????????????????(18)

關于矩陣形式下最小二乘解推導過程可以參考下面的鏈接（點擊跳轉），此處不詳細講解。

公式 (18) 即為計算結果，下面對結果進行進一步化簡：

定義沿 $A$ 到 $D^\hat D$ 的投影矩陣為： $PA=D^D^TD^TD^P_A=\frac{\hat D \hat D^T}{\hat D^T \hat D }$ ????????????????????????????????????????????????????(19)

沿 $D^\hat D$ 到 $A$ 的投影矩陣即為： $PD=I?PA=I?D^D^TD^TD^P_D=I-P_A=I-\frac{\hat D \hat D^T}{\hat D^T \hat D }$ （ $I$ 為單位矩陣）?????(20)

根據投影矩陣的性質，可以得到： $A=PD[M0]A=P_D\begin{bmatrix}M&0\end{bmatrix}$ ??????????????????????(21)

將公式 (21) 帶入 (17) 中得到： $ε=P_D M s-b$ ??????????????????????????????(22)

公式（22）的最小二乘解即為最終簡化結果。

最終簡化結果：

s=( M^T P_D M )^{-1}M^TP_Db

$M$ 只與麥克風陣列的個數與排列有關， $P_D$ 與 $b$ 在計算出時延估計的距離差 ( $d_i$ ) 后也可以求得。

六. 六元圓形麥克風陣列聲源定位

硬件平臺： Seed Studio 的六元麥克風陣列套件（官網)）

軟件平臺： Matlab 2019a

實際應用中，二維圓形陣列主要采用遠場模型較多，因為可以在保證計算精度的同時，大大降低復雜性，減少運算量，且不需要較大的陣列間距，此處使用近場模型純屬因為要做畢設😝。公式都是根據三維陣列來推導的，之所以使用二維平面陣列是因為：沒錢搭三維陣列，但又需要點實驗數據放進畢業論文中。

二維陣列缺少 $z$ 軸的基元，使用第六章中的公式求解出的聲源坐標值沒有實際意義，但是可以通過 $y$ 與 $x$ 的比值得到聲源的偏航角。

θyaw=arctanyxθ_{yaw}=arctan\frac{y}{x}

????????(23)

MATLAB 代碼：

#代碼寫得太爛，自行腦補

現以樹莓派Seed Studio麥克風陣列套件為硬件平臺進行環境聲音的獲取，將各個麥克風所獲取到的音頻儲存為wav格式文件，采樣率為44100Hz，將音頻文件在PC機中使用Matlab進行處理運算。麥克風陣列為六元均勻圓陣，屬于二維平面陣列，陣元之間的間距為4cm。

如上圖所示，為六個麥克風收集到的聲音曲線，其中橫坐標表示時間，縱坐標表示聲音幅值。麥克風錄制到的信號是一段長時間的音頻片段，因此需要進行分幀處理，此處以2048個采樣點為一幀，由于時間=采樣點數/采樣頻率，因此以2048個采樣點為一幀即以0.0464秒為一幀。我們假設一幀時間內，聲音信號是聯合平穩的，聲源相對于參考坐標系的相位運動足夠小可以忽略不計。

如上圖所示，左圖為該陣列對靜態聲源的估計曲線，右圖則為對動態聲源的估計曲線。

左圖： 在空曠且安靜的室內環境，使用手機播放歌曲作為聲源，在距離麥克風陣列0.5m處保持靜止，同時記錄下六個麥克風錄制的音頻信息，采樣率為44100Hz，以2048個采樣點為一幀，使用聲達時延法計算每一幀聲源的方位角。

右圖： 在空曠且安靜的室內環境，使用手機播放歌曲作為聲源，在距離麥克風陣列中心0.5m處繞麥克風陣列中心勻速旋轉，同時記錄下六個麥克風錄制的音頻信息，采樣率為44100Hz，以2048個采樣點為一幀，使用聲達時延法計算每一幀聲源的方位角。

七. 相關鏈接

相關鏈接：

1. 矩陣形式下最小二乘解推導：https://zhuanlan.zhihu.com/p/87582571

2. 關于投影矩陣與最小二乘：https://www.cnblogs.com/bigmonkey/p/9897047.html

3. 關于投影矩陣概念及其性質：https://blog.csdn.net/niu_123ming/article/details/86371308

4. 參考文獻：基于六元空間陣列的聲源定位系統實現

5. 參考文獻：基于麥克風陣列的室內三維聲源定位優化算法

總結

以上是生活随笔為你收集整理的基于时延法的麦克风阵列声源定位分析的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。

上一篇：麦克风阵列之一阶差分麦克风阵列
下一篇： android usb麦克风阵列,语音设