unitary_MUSIC_算法

2 MUSIC 算法

2.1 算法原理

  根据信号模型，可以得到接收信号的协方差矩阵，理想协方差矩阵如下：

$$\begin{aligned} \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{aligned}$$

其中 $\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{equation*} \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{equation*}$$

  对 $\mathbf{R}$ 进行特征值分解，可得到：

$$\begin{equation*} \begin{aligned} \mathbf{R} &= \mathbf{U} \Sigma \mathbf{U}^H = \sum\limits_{i=1}^M \lambda_i \mathbf{u}_i \mathbf{u}_i^H \\ &= \begin{bmatrix} \mathbf{U}_s & \mathbf{U}_n \end{bmatrix} \begin{bmatrix} \Sigma_s & \mathbf{O} \\ \mathbf{O} & \Sigma_n \end{bmatrix} \begin{bmatrix} \mathbf{U}_s^H \\ \mathbf{U}_n^H \end{bmatrix} \\ &= \mathbf{U}_s \Sigma_s \mathbf{U}_s^H + \mathbf{U}_n \Sigma_n \mathbf{U}_n^H \end{aligned} \end{equation*}$$

其中 $\lambda_i$ 和 $\mathbf{u}_i$ 分别为特征值和对应的特征向量，$\mathbf{O}$ 为全零矩阵，$\Sigma$ 为由 $M$ 个特征值 $\lambda_1, \cdots, \lambda_M$ 组成的对角矩阵：

$$\begin{equation*} \Sigma = \operatorname{diag} \{\lambda_1, \cdots, \lambda_M\} \end{equation*}$$

且假设：

$$\begin{equation*} \lambda_1 > \cdots > \lambda_K > \lambda_{K+1} = \cdots = \lambda_M = \sigma_n^2 \end{equation*}$$

$\Sigma_s$ 为由 $K$ 个较大特征值 $\lambda_1, \cdots, \lambda_K$ 组成的对角矩阵， $\Sigma_n$ 为由 $M-K$ 个较小特征值 $\lambda_{K+1}, \cdots, \lambda_M$ 组成的对角矩阵：

$$\begin{equation*} \Sigma_s = \operatorname{diag} \{\lambda_1, \cdots, \lambda_K\}, \Sigma_n = \operatorname{diag} \{\lambda_{K+1}, \cdots, \lambda_M\} \end{equation*}$$

$\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{equation*} \mathbf{R}\mathbf{U}_n = \mathbf{A}\mathbf{R}_s\mathbf{A}^H \mathbf{U}_n+ \sigma_n^2\mathbf{U}_n \end{equation*}$$

  同时将 $\mathbf{R} =\mathbf{U}_s \Sigma_s \mathbf{U}_s^H + \mathbf{U}_n \Sigma_n \mathbf{U}_n^H$ 左右两边乘以 $\mathbf{U}_n$ 可得：

$$\begin{equation*} \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}_n = \sigma_n^2 \mathbf{U}_n \end{equation*}$$

  联立上面两式子可得：

$$\begin{equation*} \mathbf{A}\mathbf{R}_s\mathbf{A}^H \mathbf{U}_n = \mathbf{O} \end{equation*}$$

  由于 $\mathbf{R}_s$ 和 $\mathbf{A}^H \mathbf{A}$ 均满秩，所以均可逆，因此上式两边同时乘以 $\mathbf{R}_s^{-1} (\mathbf{A}^H\mathbf{A})^{-1} \mathbf{A}^H$ 可得：

$$\begin{equation*} \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{equation*}$$

化简得：

$$\begin{equation*} \mathbf{A}^H \mathbf{U}_n = \mathbf{O} \end{equation*}$$

或者可以写成如下形式：

$$\begin{equation*} \mathbf{A}^H\mathbf{u}_i = \mathbf{0}, i = K+1, \cdots, M \end{equation*}$$

其中 $\mathbf{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 \in [-90^{\circ}, 90^{\circ}]$：

   $$\begin{equation*} P(\theta) = \frac{1}{\mathbf{a}^H(\theta)\mathbf{U}_n\mathbf{U}_n^H\mathbf{a}(\theta)} \end{equation*}$$

   此时得到一组 $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);
   lang-matlab复制代码

   2.4 参考内容

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

---

转载请注明来源，欢迎对文章中的引用来源进行考证，欢迎指出任何有错误或不够清晰的表达。可以在下面评论区评论，也可以邮件至 jaytp@qq.com