第一章 緒論
1.2 文獻回顧
估測訊號在空間中的方位，或者是入射方向與感測裝置之間的關係，此 類研究一般來說稱為Direction of arrival (DOA) estimation，DOA 的技術 層面其中一類為基於聲音到達各陣列元件的時間差(time difference of arrival, TDOA)的定位技術，TDOA的估測最常使用的方法就是取兩顆麥 克風間的GCC (Generalized Cross-Correlation)，藉由找出這個關係式的最 大值就會等效於找到了聲音傳遞到這兩顆麥克風的時間差，而GCC的定 義卻延伸出相當多的變化，也是推動此類方法持續演進的關鍵技術，這 類方法主要都是轉換到頻率域，並且定義出各自的Weighting function。
Roth[1]所提出的Roth Processor，好處在於可以壓低第一顆麥克風電子雜 訊和空間雜訊的所在頻率所產生的影響，估計的正確性。Carter等人[2]
提出的PHAT，如果兩顆麥克風彼此的雜訊的分佈情況是無關的此法可以 有相當好的效果，這也是本文中使用的方法。Carter等人[3]又提出的 SCOT，此法是同時去考慮到兩顆麥克風的雜訊所造成的影響，而想較 Roth只考慮到第一顆麥克風的雜訊所造成的影響，可說是SCOT的一個特 例。Knapp以及Carter提出的ML (Maximum Likelihood)方法[4]，放鬆了高 SNR 這樣的限制，且假設雜訊以及聲源的統計分佈都是高斯分佈
(Gaussian)，最後推演出新的Weighting function、Brandstein等人的[5]，此 方法並不需要求取channel的響應，而是利用新的統計函數(Tukey’s
Biweight)來重新定義Weighting function，並且推演出對於反射情況下更穩 健的估測方法，以上都是在探討這方面的問題。
以上方法都是用於1-D情況，藉由時間差的資訊，便可以推測出聲源相 對於這兩顆麥克風的方位，而一旦麥克風數量增加時，麥克風的位置就不 一定單單只是一個一維的情況。在多麥克風的情況下，兩兩一對的麥克風 對也就增加了，因此可以獲得的資訊也就變多，除了分別利用外，另一方
面也可以整合這些資訊來提升估測的效果，但原本的關係式為一個非線性 的方程式，所以最常用的方法是固定其中一顆麥克風當做參考點，於是問 題就可以轉換成一個線性方程式，Yao[6]中提出利用TDOA和最小平方法 (Least square)，在假設聲速不知道的情況下去估測聲源位置與聲速，更加 在[7]中提出了一套有限制的最小平方演算法，其中針對估測值之間的關係 加以限制，在最小平方法的基礎和利用已知聲速的理想範圍，加以改良估 測求得較合理的估測值，有利於提升估測的準確度。而後[8][9]也在聲速已 知的條件下，分別提出限制的最小平方演算法的封閉解。[8]中也是利用估 測值之間有個關係式，利用其當成一條限制方程式，使用Lagrange multiplier 來達成限制的最小平方演算法，但因為要求得Lagrange multiplier的係數沒 有一個封閉解，對於計算上比較麻煩，所以它提中另外一套近似的最小平 方演算法，利用泰勒展開式線性，則可找到對估測的封閉疊代解。[9]中則 是深入探討Lagrange multiplier的方法，且針對Range measurement和
Range-difference measurement分別加以討論。
探討陣列形狀的文章，從最簡單的直線分布，[10]探討了在直線不同位置 擺放，根據Power spectral dispersion，結論出要利用Stochastic Region Contraction (SRC) algorithm來找到最佳的解， [11]主要對線段上的sensor 做最佳化，不限定在直線上。[11]中，限制麥克風只能分佈在已定的線段 上，基於Cramér-Rao Lower Bound(CRBL)去最小化變異量的下界，並分 別去探討三種不同的情況，角度估測、距離估測以及位置估測，都有其 不同的擺放方式。在[12]中基於Cramér-Rao Lower Bound (CRLB)對聲源 位置的變異量做最佳化，在不限制麥克風和聲源位置的情況下，可以得 到當2維陣列以聲源為中心呈正多邊形分佈或3維陣列以聲源為中心呈正 多面體分布時，可以得到最小的位置變異量，這跟一般的直覺是一樣的。
[13]則延續[12]的工作，給予更多除了[12]提到之外的排法，更特別的探
討maximum Fisher information matrix, minimum Cramér-Rao Bound,
spherical codes, uniform angular arrays之間的關係。[14]則是基於線性方程 式係數矩陣的條件數(condition number)會影響估測誤差的上界，利用使條 件數最小化來決定感應器(sensor)的位置，但只探討最簡單的情況，也就 是4顆感應器依序分布於原點、以及3軸上，結果可看出來跟其他非條件 數最佳化的陣列相比可以有較小的誤差。[15]則是針對遠聲場的聲源，考 慮訊號及雜訊的功率譜密度(signal and noise power spectral densities)對時 間延遲(time delay)的影響，找出方位角(azimuth)和仰角(elevation)的 CRLB，進而可以推論出麥克風位置對其變異量的影響，為了最小化CRLB 而提出對麥克風位置的一些條件，於這些條件下可以決定麥克風的位 置。[16]則是也是針對遠聲場的聲源，加上在聲速已知的條件下，基於線 性方程式對azimuth和elevation的估測，分別探討azimuth和elevation的 CRLB，結論出線性方程式係數，也就是麥克風位置對CRLB的影響。
1.3 論文貢獻
本文假設聲源方位以及聲速未知的情況下，以最小方差無偏估計 (minimum variance unbiased estimator)為基準，結合奇異數(singular value)的 觀念，對於麥克風位置加以探討。以遠聲場聲源而言，對於全方位估測，
結論出一套可任意擺放和即時調整的方法，該陣列對於所有方向的聲源都 有相同的估測準度，並給予其滿足條件的方程式，更佳也證明正多邊形和 正多面體都滿足條件；對於已知的小範圍方位估測，可近一步利用縮小角 度變異量，可精確其估測準度。對近聲場聲源，對於誤差分布有進一步的 分析，分析在方程式係數也有誤差的情況下，對近聲場該如何估測，在估 測位置方面，一樣對全方位估測和區域方位估測有所分別探討，並提出一 套演算法供區域方位估測設計出精確較佳的麥克風擺放結構。
1.4 論文架構
本論文從聲源估測方法講起，再依序介紹近聲場與遠聲場的差異及分 析，最後透過模擬來比較，主要內容如下：
第二章：主要GCC演算法與聲源估測演算法及其改良方法。
第三章：從聲源估測演算法中探討奇異數與麥克風位置的關聯性，更加深 入分析如何找出最好的麥克風位置，分開遠聲場和近聲場的情況，最後為 模擬與討論。
第四張：結論與未來展望
第二章 GCC 和聲源方位估測演算法
uncorrelated。其中D就是真正的延遲(Delay)，而 代表的是改變大小的Scale value，並且其D和的改變相較於s t1( )是緩慢變化的，而此時麥克風之間的 cross correlation可以表示為：
1,2( ) 1( ) 2( )
Rx x E x t x t (2.1.3)
其中E代表期望值，而使得(2.1.3)達到最大值的 即為兩個麥克風之間的時 間延遲量，因為實際觀測的時間是有限的，所以cross correlation的估算可 以表示為：
其中T代表觀測的時間間隔，而cross correlation 和cross power spectrum 之 間的關係可以用Fourier Transform 表示成下面這個形式：
1 2 1 2 狀態，因此實際麥克風之間的cross power spectrum 可以表示為:
1 2 1 2 空間轉換函數，因此定義麥克風對之間的Generalized correlation 為:
( ) 2
, ( ) ( ) , ( )
g j f
x x g x x
R
f G f e df (2.1.7)其中g( )f H1( )f H2*( )f
實際上，在有限的時間觀察下只能以估算的Gˆx x1,2( )f 取代Gx x1,2( )f ，因此改寫 為:
( )
1 2 1 2
2
, ˆ ,
ˆx xg ( ) g( ) x x ( ) j f R f G f e df
(2.1.8)由(2.1.8)就可以估算麥克風對之間的時間延遲量，而g( )f 的選取對於時間 延遲的估計也有一定的影響，而其中一種方法是由Carter 等人提出的名為 PHAT (The Phase Transform)的方法，也就是:
1,2
( ) 1
g ( )
x x
f
G f
此方法在兩個麥克風彼此的雜訊分佈是無關的時候，有相當好的效果。
2.2 聲源方位估測演算法
因此，(2.2.3)可以改寫為:
當聲場接近遠聲場模型的時候(i 1)，由(2.2.5)式可以得到新的方程式為:
2.2.2 有限制的最小平方演算法
(2.2.20)帶入(2.2.21)式
2 3 2 2 2
( )v ( )v ( )v 0
(2.2.22) 其中
2 2 2
1 2 3
2
1 1 2 2 3 3 4
2 2 2
1 2 3 4 4
2 4
2( )
2
q q q
p q p q p q q
p p p p q
p
對於聲速v來說，這個方程式最多存在3個解，利用聲速的數值是有範圍 的，可以選取一個最適合的聲速ˆv，所以對於聲源估測可以寫成
ˆ2
y p v q
但此方法必須確定所估測的TDOA都是正確的，否則在錯誤的TDOA情況 下，他也可能算到合理的聲速，導致估測錯誤。但在TDOA估測無誤的情 況下，相較一般的LS可以改善誤差的大小。
第三章 麥克風位置擺放最佳化
V和U都必須是正交矩陣(orthogonal matrix)。
3.1.2 奇異數與麥克風位置的關聯
首先探討觀測(observation) b的誤差對於估測(estimator) x的影響 假設觀測到的ˆb有一個誤差b
利用Cramer-Rao Lower Bound(CRLB)，可以求得我們的最小方差無偏估計 (minimum variance unbiased estimator)和誤差下界(error lower bound)
ˆ N( ,2 )
利用去計算CRLB，(3.1.5)對估測值x做偏微分
2 2
跟觀測值ˆb的誤差有關，假設此刻的影響固定，想要讓估測值的變異量和 根據柯西不等式(Cauchy-Schwarz inequality)
3.2 麥克風最佳化排放位置
根據最小平方(Least square)方法
ˆs ( Ts s)1 Ts 根據(3.1.13)，如果要讓估測的變異量最少，所有的singular value要相同。
先從簡單的2維開始分析
= M是奇數跟偶數有點差別所以把(3.2.10b)分開表示
2 1 1 當M顆麥克風的時候我們的eigenvalue為底下所表示
2
( 3) / 2
由(3.2.14b)式可看出滿足最小的情況，也就是兩個eigenvalue相同。
可以看出如果是使用硬算eigenvalue的情況下，算式會很複雜，一旦擴展到 3維的空間，會變的更複雜，所以這邊考慮另一種分析方式：
L. Yu. Kolotilina的[19][20]中提到，如果有一個n n Hermitian positive definite matrix A，把他表示成 11* 12 12
12 22
2 2
如果對Ri 沒有任何限制，也就是說第1到M 1顆麥克風任意擺放後，可以
接著拓展到3維分析
以L. Yu. Kolotilina的理論來說，(3.2.25)式要滿足condition number=1，需要 以下的式子成立：
假設ai Risin
i cos
i , bi Risin
i sin
i and ci Ricos
i ，如圖二示意可稍微合併和整理(3.2.26a) ,(3.2.26b) ,(3.2.26c)
1
1-3cos 0
M
所以當然當2顆到第M顆是以正多面體的分佈時，是其中的一個特殊解，如 圖三所示。
tetrahed. cube octahed. dodecahed. icosahed.
-10
cos sin sin sin
cos
其中 2 1 所以CRLB分析的Fish imformation要改寫成
1
sin cos cos cos
1 0 ( ˆ ( ))
cos sin sin sin cos
T T
但藉由ML(maximum-likelihood)
sin cos cos cos
0
( ˆ ( ))
cos sin sin sin cos
T T
接著來分析Cramer-Rao Lower Bound，也就是誤差下界
sin cos cos sin cos sin
0 cos
sin sin
= v ( ) sin cos
cos cos sin cos sin sin
0 sin
cos cos var( )ˆ v cos cos sin cos sin ( ) sin cos
sin sin
M
如果A AT cI
cos cos var( )ˆ v cos cos sin cos sin ( ) sin cos eigenvector，由此eigenvector可反推出跟。
(A AT )1最小eigenvalue也就是A AT 的最大eigenvalue。
先分析矩陣裡面元素，目前關心的元素是11
如果考慮對聲速的估測，聲速跟x的關係為
cos sin
cos sin sin sin cos sin sin
= ( )
cos sin
= cos sin sin sin cos ( ) sin sin
cos sin ˆ
但如果是一個2維的陣列在3維空間
sin cos sin cos
sin sin sin cos
= v ( )
cos sin cos sin
cos cos sin cos
T
sin sin
M
cos cos
M
一樣也可以由已知的跟去分析陣列位置
因為這個陣列是針對90度方向設計的，可以知道對其他角度效果一定不好
cos sin
sin
可以看到跟角度的估測其實是呈現相反的情況，也就是說當估測角度好的