DOA算法1：MUSIC算法（一）

注：本博文为本人阅读论文、文章后的原创笔记，未经授权不允许任何转载或商用行为，否则一经发现本人保留追责权利。有问题可留言联系，欢迎指摘批评，共同进步！！！

一、 经典MUSIC算法

MUSIC算法全称为：MUltiple Signal Classification.是波达角(DOA)估计问题中使用最广泛的算法之一。

1.1 算法介绍(第一阶段)

模型选用加性噪声模型：
$$\begin{array}{rl}& \mathbf{x}=\mathbf{S}\mathbf{\alpha }+\mathbf{n}.\\ & \mathbf{S}=[\mathbf{s}({\varphi }_1)\mathbf{s}({\varphi }_2)\cdots ,\mathbf{s}({\varphi }_{\mathbf{M}})],\\ & \mathbf{\alpha }=[{\mathbf{\alpha }}_1,{\mathbf{\alpha }}_2\cdots {\mathbf{\alpha }}_{\mathbf{M}}{]}^{\mathbf{T}}.\end{array}$$

1. $\mathbf{s}({\varphi }_{\mathbf{i}})$是$N\times 1$大小的矩阵：$\mathbf{s}({\varphi }_{\mathbf{i}})=[s({\varphi }_{i1})s({\varphi }_{i2})\cdots s({\varphi }_{iN}){]}^T$，代表第i个信号的转向矢量(steering vector)，也可理解为第i个信号在接收阵列N个元素上产生的效果。因此矩阵$\mathbf{S}$是$N\times M$大小的矩阵，表示入射信号；
2. ${\mathbf{\alpha }}_{\mathbf{i}}$是$1\times N$大小的矩阵：${\mathbf{\alpha }}_{\mathbf{i}}=[{\alpha }_{i1}{\alpha }_{i2}\cdots {\alpha }_{iN}]$，代表信道特征。因此矩阵$\mathbf{\alpha }$是$M\times N$大小的矩阵.

假设不同信号之间互不相关。计算矩阵$\mathbf{x}$的相关矩阵$\mathbf{R}$：
$$\begin{array}{rl}\mathbf{R}& =E[{\mathbf{x}\mathbf{x}}^H],\\ & =E[\mathbf{S}{\mathbf{\alpha }\mathbf{\alpha }}^H{\mathbf{S}}^H]+E[{\mathbf{n}\mathbf{n}}^H],\\ & ={\mathbf{S}\mathbf{A}\mathbf{S}}^H+{\sigma }^2\mathbf{I},\\ & ={\mathbf{R}}_{\mathbf{s}}+{\sigma }^2\mathbf{I},\end{array}$$
其中：
$$\begin{array}{rl}{\mathbf{R}}_{\mathbf{s}}& ={\mathbf{S}\mathbf{A}\mathbf{S}}^H\\ \mathbf{A}& ={\mathbf{\alpha }\mathbf{\alpha }}^H\\ & =\left[\begin{array}{cccc}E[|{\alpha }_1{|}^2]& 0& \cdots & 0\\ 0& E[|{\alpha }_2{|}^2]& \cdots & 0\\ 0& 0& \cdots & E[|{\alpha }_M{|}^2]\end{array}\right].\end{array}$$

• 矩阵$\mathbf{A}$是一个对角阵，其中的0元素是因为$E[{\mathbf{\alpha }}_{\mathbf{i}}\cdot {\mathbf{\alpha }}_{\mathbf{j}}]=0.(i\ne j)$表示不同信号之间互不相关。
• $\mathbf{S}$是$N\times M$大小的矩阵。$\mathbf{A}$是$M\times M$大小的方阵。由此可以得到，${\mathbf{R}}_{\mathbf{s}}$是$N\times N$ 大小的方阵，且秩为M。

因为${\mathbf{R}}_{\mathbf{s}}$不是满秩的，所以它有$(N-M)$个特征向量对应于0特征值。我们就设${\mathbf{q}}_{\mathbf{m}}$为这样的特征向量之一，${\mathbf{q}}_{\mathbf{m}}$是$N\times 1$大小的列矩阵。由此可以得到如下结果：
$$\begin{array}{rl}{\mathbf{R}}_{\mathbf{s}}{\mathbf{q}}_{\mathbf{m}}& ={\mathbf{S}\mathbf{A}\mathbf{S}}^H{\mathbf{q}}_{\mathbf{m}}=0,\\ \Rightarrow {{\mathbf{q}}_{\mathbf{m}}}^H{\mathbf{S}\mathbf{A}\mathbf{S}}^H{\mathbf{q}}_{\mathbf{m}}& =0,\\ \Rightarrow {\mathbf{S}}^H{\mathbf{q}}_{\mathbf{m}}& =0\end{array}$$

MUSIC算法的基础：
由 ${\mathbf{S}}^H{\mathbf{q}}_{\mathbf{m}}=0$这一结果可得：矩阵${\mathbf{R}}_{\mathbf{s}}$的所有对应于0特征值的$(N-M)$个特征向量(${\mathbf{q}}_{\mathbf{m}}$)与所有$M$个信号转向矢量(signal steering vector)都正交。

设定：矩阵${\mathbf{Q}}_{\mathbf{n}}$由这些特征向量${\mathbf{q}}_{\mathbf{m}}$组成，${\mathbf{Q}}_{\mathbf{n}}$就是$N\times (N-M)$大小的矩阵。MUSIC算法所需计算的功率谱表达如下：
$$\begin{array}{rl}{P}_{MUSIC}(\varphi )& =\frac{1}{{\sum }_{m=1}^{N-M}|{\mathbf{s}}^H(\varphi ){\mathbf{q}}_{\mathbf{m}}{|}^2}\\ & =\frac{1}{{\mathbf{s}}^H(\varphi ){{\mathbf{Q}}_{\mathbf{n}}{\mathbf{Q}}_{\mathbf{n}}}^H\mathbf{s}(\varphi )}\\ & =\frac{1}{||{{\mathbf{Q}}_{\mathbf{n}}}^{H}s(\varphi )|{|}^2}\end{array}$$

注意到，${\mathbf{Q}}_{\mathbf{n}}$是由${\mathbf{q}}_{\mathbf{m}}$组成的，因此矩阵${\mathbf{Q}}_{\mathbf{n}}$与信号旋转矢量$\mathbf{s}(\varphi )$也是正交的，即当$\varphi$是一个波达角时，${{\mathbf{Q}}_{\mathbf{n}}}^{H}s(\varphi )=0$，这表现在$P_{MUSIC}$谱线上就是一个波峰（因为分母为0）。因此波达角的估计体现在功率谱上就是M个最大的峰值对应的$\varphi$.

1.2 算法介绍(第二阶段)

在上述MUSIC算法中有一缺陷：在现实情况下，信号的协方差矩阵${\mathbf{R}}_{\mathbf{s}}$不可知，我们最多可以估计矩阵$\mathbf{R}$(含噪声)，此二者的关键联系点就是：特征向量矩阵${\mathbf{Q}}_{\mathbf{n}}$是相同的，下面就证明这一点。

对于任意的${\mathbf{q}}_{\mathbf{m}}\in \mathbf{Q}$ ($\mathbf{Q}$是矩阵${\mathbf{R}}_{\mathbf{s}}$所有特征向量组成的矩阵，形状是$N\times N$大小的方阵)我们有：
$$\begin{array}{rl}{\mathbf{R}}_{\mathbf{s}}{\mathbf{q}}_{\mathbf{m}}& =\lambda {\mathbf{q}}_{\mathbf{m}}\\ \Rightarrow \mathbf{R}{\mathbf{q}}_{\mathbf{m}}& ={\mathbf{R}}_{\mathbf{s}}{\mathbf{q}}_{\mathbf{m}}+{\sigma }^2\mathbf{I}{\mathbf{q}}_{\mathbf{m}},\\ & =({\lambda }_m+{\sigma }^2){\mathbf{q}}_{\mathbf{m}},\end{array}$$

由此可见，矩阵${\mathbf{R}}_{\mathbf{s}}$的特征向量也同样是矩阵$\mathbf{R}$的特征向量，只不过特征值从原来的$\lambda$变成了$(\lambda +{\sigma }^2)$。注意，如前所述，矩阵${\mathbf{R}}_{\mathbf{s}}$有$(N-M)$个对应于0特征值的特征向量，它们的$\lambda =0$，对应到矩阵$\mathbf{R}$中，特征值$\lambda +{\sigma }^2={\sigma }^2$，也一共有$(N-M)$个

用矩阵${\mathbf{R}}_{\mathbf{s}}$的$N$个特征向量组成一个$N\times N$大小的对角阵$\mathbf{\Lambda }$,可得：
$$\begin{array}{rl}{\mathbf{R}}_{\mathbf{s}}\mathbf{Q}& =\mathbf{Q}\mathbf{\Lambda }\\ \Rightarrow {\mathbf{R}}_{\mathbf{s}}& ={\mathbf{Q}\mathbf{\Lambda }\mathbf{Q}}^H.\end{array}$$
$$\begin{array}{rl}\mathbf{\Lambda }& =\left[\begin{array}{ccccccc}{\lambda }_1& 0& \cdots & 0& 0& \cdots & 0\\ 0& {\lambda }_2& \cdots & 0& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮\\ 0& 0& \cdots & {\lambda }_M& 0& \cdots & 0\\ 0& 0& \cdots & 0& 0& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 0& 0& \cdots & 0\end{array}\right]\end{array}$$

同理，对于矩阵$\mathbf{R}$可得：
$$\begin{array}{rl}\mathbf{R}& =\mathbf{Q}[\mathbf{\Lambda }+{\sigma }^2{\mathbf{I}]\mathbf{Q}}^H\\ & =\mathbf{Q}\left[\begin{array}{ccccccc}{\lambda }_1+{\sigma }^2& 0& \cdots & 0& 0& \cdots & 0\\ 0& {\lambda }_2+{\sigma }^2& \cdots & 0& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮\\ 0& 0& \cdots & {\lambda }_M+{\sigma }^2& 0& \cdots & 0\\ 0& 0& \cdots & 0& {\sigma }^2& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 0& 0& \cdots & {\sigma }^2\end{array}\right]{\mathbf{Q}}^H\end{array}$$

依据特征分解，我们可以将矩阵$\mathbf{Q}$分解成信号矩阵${\mathbf{Q}}_{\mathbf{s}}$(矩阵$\mathbf{Q}$的前M列，对应M个信号特征值)和噪声特征向量矩阵${\mathbf{Q}}_{\mathbf{n}}$(矩阵$\mathbf{Q}$的后$(N-M)$列，对应$(N-M)$个噪声特征值${\sigma }^2$)

1. 矩阵${\mathbf{Q}}_{\mathbf{s}}$定义了信号子空间
2. 矩阵${\mathbf{Q}}_{\mathbf{n}}$定义了噪声子空间

注意到：矩阵${\mathbf{Q}}_{\mathbf{n}}$与矩阵${\mathbf{R}}_{\mathbf{s}}$对应于0特征值的特征向量所组成的矩阵完全相同（就是这个矩阵）。表明，矩阵$\mathbf{R}$和矩阵${\mathbf{R}}_{\mathbf{s}}$有完全相同的特征向量矩阵。

1.3 对MUSIC算法的几点观察

1. 在以上理论中，矩阵$\mathbf{R}$的最小的几个特征值都是噪声特征值，等于${\sigma }^2$,这样一来就可以区分信号和噪声特征值了：找出几个小的、相等的特征值，这就是噪声特征值，其余的是信号特征值；
2. ${\mathbf{Q}}_{\mathbf{s}}\perp {\mathbf{Q}}_{\mathbf{n}}$

依据上述两点，我们可以得出MUSIC算法的前提/基础：
所有的噪声特征值与信号转向矢量都正交
$$\begin{array}{rl}{P}_{MUSIC}(\varphi )& =\frac{1}{{\sum }_{m=M+1}^N|{\mathbf{q}}_{\mathbf{m}}\mathbf{s}(\varphi ){|}^2}\\ & =\frac{1}{{\mathbf{s}}^H(\varphi ){{\mathbf{Q}}_{\mathbf{n}}{\mathbf{Q}}_{\mathbf{n}}}^H\mathbf{s}(\varphi )},\end{array}$$

式中${\mathbf{q}}_{\mathbf{m}}$是$(N-M)$个噪声特征矩阵中的一个。如果$\varphi$就是信号的DOA之一，那么$\mathbf{s}(\varphi )\perp {\mathbf{q}}_{\mathbf{m}}$，所以$P_{MUSIC}$的分母就恒等于0。因此，DOA就是$P_{MUSIC}$的几个峰值。

二、现实情况下的MUSIC算法

2.1 现实与理论的差距

• 难点1：相关矩阵$\mathbf{R}$是未知的，它必须从接收信号中估计得来。而这项估计采用的是在多个数据快照(snapshots of data)中取平均值的方法：
  $$\mathbf{R}=\frac{1}{K}\sum _{k=1}^K{\mathbf{x}}_{\mathbf{k}}{{\mathbf{x}}_{\mathbf{k}}}^H,$$
  在参考文献[2]中，作者表明，快照数量$K>2N$时，信噪比最佳能达到3dB.

• 难点2：使用估计出来的协方差矩阵$R$会使得噪声的特征值不再相等了（原本理想情况下噪声特征值就是${\sigma }^2$）

• 难点3：在现实情况下，特征值更趋向于一个连续体(continuum)，信号和噪声特征值的区分不再明显了。

• 难点4：如果信号的数量($M$)未知，那么就不容易确定哪些特征值是相等的（至少在理论上是对应相等的，它们代表噪声空间）

2.2 在现实情况下使用MUSIC算法的步骤

1. 依据式$\mathbf{R}=\frac{1}{K}{\sum }_{k=1}^K{\mathbf{x}}_{\mathbf{k}}{{\mathbf{x}}_{\mathbf{k}}}^H$来估计出相关矩阵$\mathbf{R}$，并找到其特征分解：${\mathbf{R}=\mathbf{Q}\mathbf{\Lambda }\mathbf{Q}}^H$

2. 将矩阵$\mathbf{Q}$进行分解，得到矩阵${\mathbf{Q}}_{\mathbf{n}}$，它对应于矩阵$\mathbf{Q}$的$(N-M)$个最小（不一定相等）的特征值，表示噪声空间。

3. 得到$P_{MUSIC}(\varphi )$并绘制图像

4. $M$个信号角度就是$P_{MUSIC}(\varphi )$函数的$M$个峰值的横坐标值。

参考文献（仅写文章的标题，以做记录）

[1]. Direction of Arrival Estimation
[2]. Rapid convergence rate in adaptive arrays