【学习笔记】【DOA子空间算法】2 MUSIC 算法

2 MUSIC 算法

2.1 算法原理

  根据信号模型，可以得到接收信号的协方差矩阵，理想协方差矩阵如下：
$$\begin{array}{c}\begin{array}{rl}\mathbf{R}& =\mathrm{E}[\mathbf{x}(t){\mathbf{x}}^H(t)]\\ & =\mathbf{A}\mathrm{E}[\mathbf{s}(t){\mathbf{s}}^H(t)]{\mathbf{A}}^H+\mathrm{E}[\mathbf{n}(t){\mathbf{n}}^H(t)]\\ & =\mathbf{A}{\mathbf{R}}_S{\mathbf{A}}^H+{\sigma }_n^2\mathbf{I}\end{array}\end{array}$$
其中 ${\mathbf{R}}_S$ 和 ${\sigma }_n^2\mathbf{I}$ 分别代表信号协方差矩阵和噪声协方差矩阵，${\sigma }_n^2$ 为噪声方差， $\mathbf{I}$ 为单位矩阵。需要注意到的是，由于信号源之间是独立的，${\mathbf{R}}_S$ 是对角非奇异矩阵；$\mathbf{R}$ 是非奇异的，同时 $\mathbf{R}={\mathbf{R}}^H$，因此 $\mathbf{R}$ 是 Hermitian 矩阵且是正定矩阵。
  但因为采样数 $T$ 有限，理想协方差矩阵无法直接获得，此时通常用样本的协方差矩阵 $\hat{\mathbf{R}}$ 来代替：
$$\begin{array}{c}\mathbf{R}\triangleq \hat{\mathbf{R}}=\frac{1}{T}\sum _{t=1}^T\mathbf{x}(t){\mathbf{x}}^H(t)=\frac{1}{T}\mathbf{X}{\mathbf{X}}^H\end{array}$$
  对 $\mathbf{R}$ 进行特征值分解，可得到：
$$\begin{array}{c}\begin{array}{rl}\mathbf{R}& =\mathbf{U}\Sigma {\mathbf{U}}^H=\sum _{i=1}^M{\lambda }_i{\mathbf{u}}_i{\mathbf{u}}_i^H\\ & =\left[\begin{array}{cc}{\mathbf{U}}_S& {\mathbf{U}}_N\end{array}\right]\left[\begin{array}{cc}{\Sigma }_S& \mathbf{O}\\ \mathbf{O}& {\Sigma }_N\end{array}\right]\left[\begin{array}{c}{\mathbf{U}}_S^H\\ {\mathbf{U}}_N^H\end{array}\right]\\ & ={\mathbf{U}}_S{\Sigma }_S{\mathbf{U}}_S^H+{\mathbf{U}}_N{\Sigma }_N{\mathbf{U}}_N^H\end{array}\end{array}$$
其中 ${\lambda }_i$ 和 ${\mathbf{u}}_i$ 分别为特征值和对应的特征向量，$\mathbf{O}$ 为全零矩阵，$\Sigma$ 为由 $M$ 个特征值 ${\lambda }_1,\cdots ,{\lambda }_M$ 组成的对角矩阵：
Σ S = diag ⁡ { λ 1 , ⋯   , λ M } \begin{equation*} \Sigma_S = \operatorname{diag} \{\lambda_1, \cdots, \lambda_M\} \end{equation*} ΣS​=diag{λ1​,⋯,λM​}​
且假设：
$$\begin{array}{c}{\lambda }_1>{\lambda }_2>\cdots >{\lambda }_K>{\lambda }_{K+1}>\cdots >{\lambda }_M={\sigma }_n^2\end{array}$$
${\Sigma }_S$ 为由 $K$ 个较大特征值 ${\lambda }_1,\cdots ,{\lambda }_K$ 组成的对角矩阵， ${\Sigma }_N$ 为由 $M-K$ 个较小特征值 ${\lambda }_{K+1},\cdots ,{\lambda }_M$ 组成的对角矩阵：
Σ S = diag ⁡ { λ 1 , ⋯   , λ K } , Σ N = diag ⁡ { λ K + 1 , ⋯   , λ M } \begin{equation*} \Sigma_S = \operatorname{diag} \{\lambda_1, \cdots, \lambda_K\}, \Sigma_N = \operatorname{diag} \{\lambda_{K+1}, \cdots, \lambda_M\} \end{equation*} ΣS​=diag{λ1​,⋯,λK​},ΣN​=diag{λK+1​,⋯,λM​}​
${\mathbf{U}}_S$ 表示由 $K$ 个较大特征值 ${\lambda }_1,\cdots ,{\lambda }_K$ 对应的特征矢量张成的信号子空间，${\mathbf{U}}_N$ 表示由 $M-K$ 个较小特征值 ${\lambda }_{K+1},\cdots ,{\lambda }_M$ 对应的特征矢量张成的噪声子空间。注意，$\mathbf{U}$ 为酉矩阵。
  将 $\mathbf{R}=\mathbf{A}{\mathbf{R}}_S{\mathbf{A}}^H+{\sigma }_n^2\mathbf{I}$ 左右两边乘以 ${\mathbf{U}}_N$ 可得：
$$\begin{array}{c}\mathbf{R}{\mathbf{U}}_N=\mathbf{A}{\mathbf{R}}_S{\mathbf{A}}^H{\mathbf{U}}_N+{\sigma }_n^2{\mathbf{U}}_N\end{array}$$
  同时将 $\mathbf{R}={\mathbf{U}}_S{\Sigma }_S{\mathbf{U}}_S^H+{\mathbf{U}}_N{\Sigma }_N{\mathbf{U}}_N^H$ 左右两边乘以 ${\mathbf{U}}_N$ 可得：
$$\begin{array}{c}\mathbf{R}{\mathbf{U}}_N={\mathbf{U}}_S{\Sigma }_S{\mathbf{U}}_S^H{\mathbf{U}}_N+{\mathbf{U}}_N{\Sigma }_N{\mathbf{U}}_N^H{\mathbf{U}}_N=\mathbf{O}+{\sigma }_n^2{\mathbf{U}}_{\mathbf{N}}={\sigma }_n^2{\mathbf{U}}_{\mathbf{N}}\end{array}$$
  联立上面两式子可得：
$$\begin{array}{c}\mathbf{A}{\mathbf{R}}_S{\mathbf{A}}^H{\mathbf{U}}_N=\mathbf{O}\end{array}$$
  由于 ${\mathbf{R}}_S$ 和 ${\mathbf{A}}^H\mathbf{A}$ 均满秩，所以均可逆，因此上式两边同时乘以 ${\mathbf{R}}_S^{-1}({\mathbf{A}}^H\mathbf{A}{)}^{-1}{\mathbf{A}}^H$ 可得：
$$\begin{array}{c}{\mathbf{R}}_S^{-1}({\mathbf{A}}^H\mathbf{A}{)}^{-1}{\mathbf{A}}^H\mathbf{A}{\mathbf{R}}_S{\mathbf{A}}^H{\mathbf{U}}_N={\mathbf{R}}_S^{-1}({\mathbf{A}}^H\mathbf{A}{)}^{-1}{\mathbf{A}}^H\mathbf{O}\end{array}$$
化简得：
$$\begin{array}{c}{\mathbf{A}}^H{\mathbf{U}}_N=\mathbf{O}\end{array}$$
或者可以写成如下形式：
$$\begin{array}{c}{\mathbf{A}}^H{\mathbf{u}}_i=0,i=K+1,\cdots ,M\end{array}$$
其中 $0$ 为全零向量。上式表明：噪声特征值所对应的特征向量 ${\mathbf{u}}_i$ 与方向矢量矩阵 $\mathbf{A}$ 的列向量正交，而 $\mathbf{A}$ 的各列与信号源的方向相对应的，因此可以利用噪声子空间求解信号源方向。

2.2 算法步骤

  MUSIC 算法步骤如下（输入为阵列接收矩阵 $\mathbf{X}$）：

1. 计算协方差矩阵 $\mathbf{R}=\frac{1}{T}\mathbf{X}{\mathbf{X}}^H$。
2. 对 $\mathbf{R}$ 进行特征值分解，并对特征值进行排序，然后取得 $M-K$ 个较小特征值对应的特征向量来组成噪声子空间 ${\mathbf{U}}_N$。
3. 以下式遍历 $\theta =0^{\circ },1^{\circ }\cdots ,180^{\circ }$：
   $$\begin{array}{c}P(\theta )=\frac{1}{{\mathbf{a}}^H(\theta ){\mathbf{U}}_N{\mathbf{U}}_N^H\mathbf{a}(\theta )}\end{array}$$
   此时得到一组 $P(\theta )$，$K$ 个最大值对应的 $\theta$ 就是需要返回的结果。

2.3 代码实现

% music.m
clear;
clc;
close all;

%% 参数设定
c = 3e8;                                              % 光速
fc = 500e6;                                           % 载波频率
lambda = c/fc;                                        % 波长
d = lambda/2;                                         % 阵元间距，可设 2*d = lambda
twpi = 2.0*pi;                                        % 2pi
derad = pi/180;                                       % 角度转弧度
theta = [-20, 30]*derad;                              % 待估计角度
idx = 0:1:7; idx = idx';                              % 阵元位置索引

M = length(idx);                                      % 阵元数
K = length(theta);                                    % 信源数
T = 512;                                              % 快拍数
SNR = 0;                                              % 信噪比

%% 信号模型建立
S = randn(K,T) + 1j*randn(K,T);                       % 复信号矩阵S，维度为K*T
% A = exp(-1j*twpi*d*idx*sin(theta)/lambda);           % 方向矢量矩阵A，维度为M*K
A = exp(-1j*pi*idx*sin(theta));                       % 2d = lambda，直接忽略不写
X = A*S;                                              % 接收矩阵X，维度为M*T
X = awgn(X,SNR,'measured');                           % 添加噪声

%% MUSIC 算法
% 计算协方差矩阵
R = X*X'/T;
% 特征值分解并取得噪声子空间
[U,D] = eig(R);                                       % 特征值分解
[D,I] = sort(diag(D));                                % 将特征值排序从小到大
U = fliplr(U(:, I));                                  % 对应特征矢量排序，fliplr 之后，较大特征值对应的特征矢量在前面
Un = U(:, K+1:M);                                     % 噪声子空间
Gn = Un*Un';
% 空间谱搜索
searchGrids = 0.1;                                    % 搜索间隔
ang = (-90:searchGrids:90)*derad;
Pmusic = zeros(1, length(ang));                       % 空间谱
for i = 1:length(ang)
    a = exp(-1j*pi*idx*sin(ang(i)));
    Pmusic(:, i) = 1/(a'*Gn*a);
end
% 归一化处理，单位化为 db
Pmusic = abs(Pmusic);
Pmusic = 10*log10(Pmusic/max(Pmusic));
% 作图
figure;
plot(ang/derad, Pmusic, '-', 'LineWidth', 1.5);
set(gca, 'XTick',(-90:30:90));
xlabel('\theta(\circ)', 'FontSize', 12, 'FontName', '微软雅黑');
ylabel('空间谱(dB)', 'FontSize', 12, 'FontName', '微软雅黑');
% 找出空间谱Pmusic中的峰值并得到其对应角度
[pks, locs] = findpeaks(Pmusic);                      % 找极大值
[pks, id] = sort(pks);
locs = locs(id(end-K+1:end))-1;
Theta = locs*searchGrids - 90;
Theta = sort(Theta);
disp('估计结果：');
disp(Theta);

2.4 参考内容

1. 【博客园】子空间分析方法
2. 【博客园】空间谱专题10：MUSIC算法
3. 【CSDN】子空间方法——MUSIC算法
4. 【CSDN】Traditional Comm笔记【5】：基于MUSIC Algorithm的DoA/AoA估计以及MATLAB实现
5. 【知乎】DoA 估计：多重信号分类 MUSIC 算法（附 MATLAB 代码）