较为详细的MUSIC算法原理及MATLAB实现

DOA（Direction Of Arrival）波达方向定位技术主要有ARMA谱分析、最大似然法、熵谱分析法和特征分解法，特征分解法主要有MUSIC算法、ESPRIT算法WSF算法等。 MUSIC (Multiple Signal Classification)算法，即多信号分类算法，由Schmidt等人于1979年提出。MUSIC算法是一种基于子空间分解的算法，它利用信号子空间和噪声子空间的正交性，构建空间谱函数，通过谱峰搜索，估计信号的参数。对于声源定位来说，需要估计信号的DOA。MUSIC算法对DOA的估计有很高的分辨率，且对麦克风阵列的形状没有特殊要求，因此应用十分广泛。

运用矩阵的定义，可得到更为简洁的表达式：

X = A S + N X=AS+N X=AS+N

式中

X = [ x 1 ( t ) , x 2 ( t ) , . . . x M ( t ) ] T X=[x_1(t),x_2(t),...x_M(t)]^T X=[x1​(t),x2​(t),...xM​(t)]T， S = [ S 1 ( t ) , S 2 ( t ) , . . . S D ( t ) ] T S=[S_1(t),S_2(t),...S_D(t)]^T S=[S1​(t),S2​(t),...SD​(t)]T， A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T A=[a(\theta_1),a(\theta_2),...a(\theta_D)]^T A=[a(θ1​),a(θ2​),...a(θD​)]T， N = [ n 1 ( t ) , n 2 ( t ) , . . . n M ( t ) ] T N=[n_1(t),n_2(t),...n_M(t)]^T N=[n1​(t),n2​(t),...nM​(t)]T。

X X X为阵元的输出， A A A为方向响应向量， S S S是入射信号， N N N表示阵列噪声。

其中 φ k = 2 π d λ s i n θ k \varphi_k=\frac{2\pi d}{\lambda}sin\theta_k φk​=λ2πd​sinθk​有

A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T = [ 1 1 ⋯ 1 e − j φ 1 e − j φ 2 ⋯ e − j φ D ⋮ ⋮ ⋱ ⋮ e − j ( M − 1 ) φ 1 e − j ( M − 1 ) φ 2 ⋯ e − j ( M − 1 ) φ D ] A=[a(\theta_1),a(\theta_2),...a(\theta_D)]^T=\left[ \begin{matrix} 1 & 1 & \cdots & 1 \\ {e^{-j\varphi_1}} & {e^{-j\varphi_2}} & \cdots & {e^{-j\varphi_D}} \\ \vdots & \vdots & \ddots & \vdots \\ {e^{-j(M-1)\varphi_1}} & {e^{-j(M-1)\varphi_2}} & \cdots & {e^{-j(M-1)\varphi_D}} \\ \end{matrix} \right] A=[a(θ1​),a(θ2​),...a(θD​)]T=⎣⎢⎢⎢⎡​1e−jφ1​⋮e−j(M−1)φ1​​1e−jφ2​⋮e−j(M−1)φ2​​⋯⋯⋱⋯​1e−jφD​⋮e−j(M−1)φD​​⎦⎥⎥⎥⎤​ 对 x m ( t ) x_m(t) xm​(t)进行N点采样，要处理的问题就变成了通过输出信号 x m ( t ) x_m(t) xm​(t)的采样 { x m ( i ) = 1 , 2 , . . . , M } \{ x_m (i)=1,2,...,M\} {xm​(i)=1,2,...,M}估计信号源的波达方向角 θ 1 , θ 2 . . . θ D \theta_1,\theta_2...\theta_D θ1​,θ2​...θD​,由此可以很自然的将阵列信号看作是噪声干扰的若干空间谐波的叠加，从而将波达方向估计问题与谱估计联系起来。 对阵列输出X做相关处理，得到其协方差矩阵

R x = E [ X X H ] R_x=E[XX^H] Rx​=E[XXH]

其中 H H H表示矩阵的共轭转置。 根据已假设信号与噪声互不相关、噪声为零均值白噪声，因此可得到：

R x = E [ ( A S + N ) ( A S + N ) H ] = A E [ S S H ] A H + E [ N N H ] = A R S A H + R N R_x=E[(AS+N)(AS+N)^H] =AE[SS^H]A^H+E[NN^H]=AR_SA^H+R_N Rx​=E[(AS+N)(AS+N)H]=AE[SSH]AH+E[NNH]=ARS​AH+RN​ 其中 R s = E [ S S H ] R_s=E[SS^H] Rs​=E[SSH]称为信号相关矩阵

R N = σ 2 I R_N=\sigma^2I RN​=σ2I是噪声相关阵

σ 2 \sigma^2 σ2是噪声功率

I I I是 M × M M\times M M×M阶的单位矩阵

在实际应用中通常无法直接得到 R x R_x Rx​，能使用的只有样本的协方差矩阵： R x ^ = 1 N ∑ i = 1 N X ( i ) X H ( i ) \hat{R_x}=\frac{1}{N} \sum_{i=1}^{N}X(i)X^H (i) Rx​^​=N1​∑i=1N​X(i)XH(i)， R x ^ \hat{R_x} Rx​^​是 R x R_x Rx​的最大似然估计。

当采样数 N → ∞ N\to\infty N→∞，他们是一致的，但实际情况将由于样本数有限而造成误差。根据矩阵特征分解的理论，可对阵列协方差矩阵进行特征分解，首先考虑理想情况，即无噪声的情况： R x = A R s A H R_x=AR_sA^H Rx​=ARs​AH，对均匀线阵，矩阵A由 A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T = [ 1 1 ⋯ 1 e − j φ 1 e − j φ 2 ⋯ e − j φ D ⋮ ⋮ ⋱ ⋮ e − j ( M − 1 ) φ 1 e − j ( M − 1 ) φ 2 ⋯ e − j ( M − 1 ) φ D ] A=[a(\theta_1),a(\theta_2),...a(\theta_D)]^T=\left[ \begin{matrix} 1 & 1 & \cdots & 1 \\ {e^{-j\varphi_1}} & {e^{-j\varphi_2}} & \cdots & {e^{-j\varphi_D}} \\ \vdots & \vdots & \ddots & \vdots \\ {e^{-j(M-1)\varphi_1}} & {e^{-j(M-1)\varphi_2}} & \cdots & {e^{-j(M-1)\varphi_D}} \\ \end{matrix} \right] A=[a(θ1​),a(θ2​),...a(θD​)]T=⎣⎢⎢⎢⎡​1e−jφ1​⋮e−j(M−1)φ1​​1e−jφ2​⋮e−j(M−1)φ2​​⋯⋯⋱⋯​1e−jφD​⋮e−j(M−1)φD​​⎦⎥⎥⎥⎤​ 所定义的范德蒙德矩阵，只要满足 θ i ≠ θ j , i ≠ j \theta_i\neq \theta_j,i\neq j θi​​=θj​,i​=j，则他的各列相互独立。

若 R s R_s Rs​为非奇异矩阵 R a n k ( R s ) = D Rank(R_s)=D Rank(Rs​)=D，各信号源两两不相干，且 M > D M>D M>D，则 r a n d ( A R s A H ) = D rand(AR_sA^H)=D rand(ARs​AH)=D，

由于 R x = E [ X X H ] R_x=E[XX^H] Rx​=E[XXH]，有：

R s H = R x R_s^H=R_x RsH​=Rx​

即 R s R_s Rs​为Hermite矩阵，它的特性是都是实数，又由于 R s R_s Rs​为正定的，因此 A R s A … … H AR_sA……H ARs​A……H为半正定的，它有D个正特征值和 M − D M-D M−D个零特征值。

再考虑有噪声存在的情况

R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2I Rx​=ARs​AH+σ2I

由于 σ 2 > 0 \sigma^2>0 σ2>0， R x R_x Rx​为满秩阵，所以 R x R_x Rx​有M个正实特征值 λ 1 , λ 2 . . . λ M \lambda_1,\lambda_2...\lambda_M λ1​,λ2​...λM​

分别对应于M个特征向量 v 1 , v 2 . . . v M v_1,v_2...v_M v1​,v2​...vM​。又由于 R x R_x Rx​为Hermite矩阵，所以各特征向量是正交的，即： v i H v j = 0 , i ≠ j v_i^Hv_j=0,i\neq j viH​vj​=0,i​=j与信号有关的特征值只有D个,分别等于矩阵 A R s A H AR_sA^H ARs​AH的各特征值与 σ 2 \sigma^2 σ2之和，其余 M − D M-D M−D个特征值为 σ 2 \sigma^2 σ2，即 σ 2 \sigma^2 σ2为 R R R的最小特征值，它是 M − D M-D M−D维的，对应的特征向量 v i , i = 1 , 2 , . . . , M v_i,i=1,2,...,M vi​,i=1,2,...,M中，也有D个是与信号有关的，另外 M − D M-D M−D个是与噪声有关的，可利用特征分解的性质求出信号源的波达方向 θ k \theta_k θk​。

通过对协方差矩阵的特征值分解，可得到如下结论：

将矩阵 R x R_x Rx​的特征值进行从小到大的排序，即 λ 1 ≥ λ 2 ≥ . . . ≥ λ M > 0 \lambda_1 \geq \lambda_2\geq...\geq\lambda_M>0 λ1​≥λ2​≥...≥λM​>0，其中D个较大的特征值对应于信号， M − D M-D M−D个较小的特征值对应于噪声。

矩阵 R x R_x Rx​的属于这些特征值的特征向量也分别对应于各个信号和噪声，因此可把 R x R_x Rx​的特征值（特征向量）划分为信号特征（特征向量）与噪声特征（特征向量）。

设 λ i \lambda_i λi​为 R x R_x Rx​的第 i i i个特征值， v i v_i vi​是与 λ i \lambda_i λi​个相对应的特征向量，有：

R x v i = λ i v i R_xv_i=\lambda_iv_i Rx​vi​=λi​vi​

再设 λ i = σ 2 \lambda_i=\sigma^2 λi​=σ2是 R x R_x Rx​的最小特征值 R x v i = σ 2 v i i = D + 1 , D + 2... M R_xv_i=\sigma^2v_i i=D+1,D+2...M Rx​vi​=σ2vi​i=D+1,D+2...M，

将 R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2I Rx​=ARs​AH+σ2I代入可得 σ 2 v i = ( A R s A H + σ 2 I ) v i \sigma^2v_i=(AR_sA^H+\sigma^2I)v_i σ2vi​=(ARs​AH+σ2I)vi​，

将其右边展开与左边比较得：

A R s A H v i = 0 AR_sA^Hv_i=0 ARs​AHvi​=0

因 A H A A^HA AHA是 D ∗ D D*D D∗D维的满秩矩阵， ( A H A ) − 1 (A^HA)^{-1} (AHA)−1存在；

而 R s − 1 R_s^{-1} Rs−1​同样存在，则上式两边同乘以 R s − 1 ( A H A ) − 1 A H R_s^{-1}(A^HA)^{-1}A^H Rs−1​(AHA)−1AH，

有：

R s − 1 ( A H A ) − 1 A H A R s A H v i = 0 R_s^{-1}(A^HA)^{-1}A^HAR_sA^Hv_i=0 Rs−1​(AHA)−1AHARs​AHvi​=0

于是有

A H v i = 0 ， i = D + 1 , D + 2 , . . . , M A^Hv_i=0，i=D+1,D+2,...,M AHvi​=0，i=D+1,D+2,...,M

上式表明：噪声特征值所对应的特征向量（称为噪声特征向量） v i v_i vi​，与矩阵 A A A的列向量正交，而 A A A的各列是与信号源的方向相对应的，这就是利用噪声特征向量求解信号源方向的出发点。

用各噪声特征向量为例，构造一个噪声矩阵 E n E_n En​:

E n = [ v D + 1 , v D + 2 , . . . v M ] E_n=[v_{D+1},v_{D+2},...v_{M}] En​=[vD+1​,vD+2​,...vM​]

定义空间谱 P m u ( θ ) P_{mu}(\theta) Pmu​(θ):

P m u ( θ ) = 1 a H ( θ ) E n E n H a ( θ ) = 1 ∥ E n H a ( θ ) ∥ 2 P_{mu}(\theta)=\frac{1}{{a^H}(\theta)}E_nE_n^Ha(\theta)=\frac{1}{\Vert E_n^Ha(\theta)\Vert^2} Pmu​(θ)=aH(θ)1​En​EnH​a(θ)=∥EnH​a(θ)∥21​

该式中分母是信号向量和噪声矩阵的内积，当 a ( θ ) a(\theta) a(θ)和 E n E_n En​的各列正交时，该分母为零，但由于噪声的存在，它实际上为一最小值，因此 P m u ( θ ) P_{mu}(\theta) Pmu​(θ)有一尖峰值，由该式，使 θ \theta θ变化，通过寻找波峰来估计到达角。

1.根据N个接收信号矢量得到下面协方差矩阵的估计值：

R x = 1 N ∑ i = 1 N X ( i ) X H ( i ) R_x=\frac{1}{N} \sum_{i=1}^NX(i)X^H(i) Rx​=N1​∑i=1N​X(i)XH(i)

对上面得到的协方差矩阵进行特征分解

R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2I Rx​=ARs​AH+σ2I

2.按特征值的大小排序将与信号个数 D D D相等的特征值和对应的特征向量看做信号部分空间，将剩下的 M − D M-D M−D个特征值和特征向量看做噪声部分空间，得到噪声矩阵 E n E_n En​:

A H v i = 0 ， i = D + 1 , D + 2 , . . . , M A^Hv_i=0，i=D+1,D+2,...,M AHvi​=0，i=D+1,D+2,...,M

E n = [ v D + 1 , v D + 2 , . . . v M ] E_n=[v_{D+1},v_{D+2},...v_{M}] En​=[vD+1​,vD+2​,...vM​]

3.使 θ \theta θ变化，按照式

P m u ( θ ) = 1 a H ( θ ) E n E n H a ( θ ) P_{mu}(\theta)=\frac{1}{{a^H}(\theta)E_nE_n^Ha(\theta)} Pmu​(θ)=aH(θ)En​EnH​a(θ)1​

来计算谱函数，通过寻求峰值来得到波达方向的估计值。