MUSIC(Multiple Signal Classification)算法

MUSIC算法是一种基于子空间分解的算法，它利用信号子空间和噪声子空间的正交性，构建空间谱函数，通过谱峰搜索，估计信号的参数。MUSIC算法对DOA的估计有很高的分辨率，且对阵列的形状没有要求，应用十分广泛。

System model

对于muti-user的情况，假设有$M$根天线，$K$个用户，系统建模为

$\boldsymbol{y}(t)=\boldsymbol{A}s(t)+n(t)$

其中$\boldsymbol{A}=\left[a(\theta_1),a(\theta_2),\cdots，a(\theta_K)\right]\in\mathbb{C}^{M\times K}$是K个用户的导向向量，要求是LOS信道，这里导向向量是远场的形式，其实是可以拓展到近场的，在极域进行二维的MUSIC。

$s(t)\in\mathbb{C}^{K\times L}$是导频数据，L是snapshot的数量。这里我们要求不同用户的导频相互独立，那么导频的协方差是一个满秩的矩阵，这一点至关重要。

$n(t)\sim\mathcal{CN}(0,\sigma_z^2\boldsymbol{I}_r)$是接收端的噪声矢量，服从circularly symmetric complex Gaussian分布。

Method

The covariance matrix of the received signal is given by

$\mathbf{R}=\mathbb{E}\left\{\boldsymbol{y}(t)\boldsymbol{y}^\mathrm{H}(t)\right\}=\mathbf{A}\mathbf{S}\mathbf{A}^\mathrm{H}+\sigma^2\mathbf{I}_M$

在实际应用中通常无法直接得到$\mathbf{R}$，能使用的只有样本的协方差矩阵：

$\mathbf{\hat{R}}=\frac{1}{L} \sum_{t=1}^{L} \boldsymbol{y}(t)\boldsymbol{y}^\mathrm{H}(t)$

snapshot数量足够大时，$\mathbf{\hat{R}}$近似等于$\mathbf{R}$。
$\mathbf{A}\mathbf{S}\mathbf{A}^\mathrm{H}$是半正定的Hermit矩阵，因为用户角度不同的话$\mathbf{A}$是满秩的，$\mathbf{S}$也是满秩的。对它做特征值分解：

$\mathbf{A}\mathbf{S}\mathbf{A}^\mathrm{H}=\mathbf{U}\dot{\Sigma}\mathbf{U}^\mathrm{H}$

这里$\mathbf{S}$的秩为$K$，$\mathbf{A}\mathbf{S}\mathbf{A}^\mathrm{H}$的大小为$M\times M$，秩为$K$，是非满秩的，$\dot{\Sigma}=\text{diag}(\lambda_1,\lambda_2,\cdots,\lambda_K,0,0,\cdots,0)$。那么$\mathbf{R}$的特征值分解可以写成

$\mathbf{R}=\mathbf{ASA}^{\mathrm{H}}+\sigma^{2}\mathbf{I}_{M}=\mathbf{U}\underbrace{\left(\dot{\boldsymbol{\Sigma}}+\sigma^{2}\mathbf{I}_{M}\right)}_{=\Sigma}\mathbf{U}^{\mathrm{H}}=\mathbf{U}_{s}\mathbf{\Sigma}_{s}\mathbf{U}_{s}^{\mathrm{H}}+\mathbf{U}_{n}\mathbf{\Sigma}_{n}\mathbf{U}_{n}^{\mathrm{H}}$

其中$\mathbf{U}_s$和$\mathbf{U}_n$分别是噪声空间和信号空间的特征向量，是酉矩阵。为了方便说明，给出$\mathbf{\Sigma}$的具体表达

$\mathbf{\Sigma}=\begin{bmatrix}&\lambda_1+\sigma^2&0&\dots&0&0&\dots&0\\&0&\lambda_2+\sigma^2&\dots&0&0&\dots&0\\&\dots&\dots&\dots&\dots&\dots&\dots\\&0&0&\dots&\lambda_K+\sigma^2&0&\dots&0\\&0&0&\dots&0&\sigma^2&\dots&0\\&\dots&\dots&\dots&\dots&\dots&\dots\\&0&0&\dots&0&0&\dots&\sigma^2\end{bmatrix}$

前$K$次是被噪声污染的能量，后$M-K$次是是神能量，但注意高SNR下，信号能量远大于噪声能量。
下面证明导向矢量与噪声空间的正交性，让$\mathbf{R}$右乘噪声空间

$\begin{aligned} \mathbf{RU}_n=\mathbf{ASA}^\mathrm{H}\mathbf{U}_n+\sigma^2\mathbf{U}_n& =\mathbf{U}_s\boldsymbol{\Sigma}_s\mathbf{U}_s^\mathrm{H}\mathbf{U}_n+\mathbf{U}_n\boldsymbol{\Sigma}_n\mathbf{U}_n^\mathrm{H}\mathbf{U}_n \\ &=\mathbf{U}_n\boldsymbol{\Sigma}_n=\sigma^2\mathbf{U}_n \end{aligned}$

那么可以得到

$\mathbf{ASA}^\mathrm{H}\mathbf{U}_n=\mathrm{0}$

两边同乘$\mathbf{S}^\mathrm{H}$

$\mathbf{S}^\mathrm{H}\mathbf{ASA}^\mathrm{H}\mathbf{U}_n=\mathrm{0}$

由于$\mathbf{Q}x=0$和$\mathbf{Q}^\mathrm{H}\mathbf{Q}x=0$同解，那么可以得到

$\mathbf{SA}^\mathrm{H}\mathbf{U}_n=\mathrm{0}$

由于$\mathbf{S}$是满秩的方阵，可逆，最终可以得到

$\mathbf{A}^\mathrm{H}\mathbf{U}_n=\mathrm{0}$

为了便于谱峰搜素，把0向量标量化放分母上去

$P_{\text{MUSIC}}^{\text{far}}(\theta)=\frac{1}{\mathbf{a}^{\text{H}}(\theta)\hat{\mathbf{U}}_{n}\hat{\mathbf{U}}_{n}^{\text{H}}\mathbf{a}(\theta)}$

当角度遍历到用户时，出现谱峰，实现DOA估计。注意这算法对于阵列结构和是否远近场没有要求。算法要求导频数据协方差矩阵满秩，并且能估计最大的用户数量是天线数量减一，因为要分出噪声空间。