Music算法介绍

MUSIC算法，全称为Multiple Signal Classification（多重信号分类），是用于波达方向估计的一种高分辨率谱估计算法。它最初由Schmidt在1979年提出，广泛应用于雷达、声纳、无线通信等领域中，以确定多个信号源的方向。该算法利用了信号子空间和噪声子空间正交性的特点，构造噪声空间然后通过谱峰搜索来检测信号的波达方向。

1. MUSIC 算法步骤

输入：多快拍接收数据 $Y\in {\mathbb{C}}^{M\times L}$
1）计算采样协方差矩阵 $R$；
2）对 $R$ 进行特征分解，求出对应特征值及特征矢量；
3）按照特征值大小排序，根据特征值判断信号数 $K$ 并确定信号子空间 $U_S$ 与噪声子空间 $U_N$；
4）通过 $P(\theta )=\frac{1}{a(\theta {)}^H{U}_N{U}_N^{H}a(\theta )}$ 得到空间谱估计；
5）对空间谱进行谱峰搜索，前 $K$ 个极大值点对应的角度就是目标入射方向。
输出：目标入射角度

2. MUSIC 算法原理证明

首先，将观测数据向量的协方差矩阵特征分解
$$\begin{array}{r}{R}_{xx}=A{R}_{ss}{A}^H+{\sigma }^{2}I=U\Sigma U^H=\left[\begin{array}{cc}{U}_S& U_N\end{array}\right]\left[\begin{array}{cc}{\Sigma }_S& O\\ O& {\Sigma }_N\end{array}\right]\left[\begin{array}{c}{U}_S^H\\ U_N^H\end{array}\right]\end{array}$$

右乘 $U_N$ 有

$$\begin{array}{r}{R}_{xx}{U}_N=\left[\begin{array}{c}{U}_S,U_N\end{array}\right]\left[\begin{array}{cc}{\Sigma }_S& O\\ O& {\Sigma }_N\end{array}\right]\left[\begin{array}{c}{U}_S^H\\ U_N^H\end{array}\right]U_N=\left[\begin{array}{c}{U}_S,U_N\end{array}\right]\left[\begin{array}{cc}{\Sigma }_S& O\\ O& {\Sigma }_N\end{array}\right]\left[\begin{array}{c}O\\ I\end{array}\right]={\sigma }^2{U}_N\end{array}$$

以及

$$\begin{array}{r}{R}_{xx}{U}_N=A{R}_{ss}{A}^H{U}_N+{\sigma }^2{U}_N\end{array}$$

联立两式得到

$$A{R}_{ss}{A}^H{U}_N=O$$

左乘 $U_N^H$ 有

$$U_N^{H}A{R}_{ss}{A}^H{U}_N=O$$

若 Q 非奇异时 $t^{H}Qt=0$，则 $t=0$ 。故当$R_{SS}$非奇异，也就是入射信号不相干时，有

$$A^H{U}_N=O$$
也就是
$${\left[\begin{array}{c}a({\omega }_1),a({\omega }_2),\dots ,a({\omega }_N)\end{array}\right]}^H{U}_N=O$$
即
$$\begin{array}{r}a(\omega {)}^H{U}_N=0^T,\omega ={\omega }_1,{\omega }_2,\dots ,{\omega }_N\\ a(\omega {)}^H{U}_N\ne 0^T,\omega \ne {\omega }_1,{\omega }_2,\dots ,{\omega }_N\end{array}$$
写成标量形式，定义一种类似功率谱的函数
$$P(\omega )=\frac{1}{a(\omega {)}^H{U}_N{U}_N^{H}a(\omega )}$$
上式的N个峰值对应N个信号的来波方向

3. MATLAB仿真

clear all
close all

theta = [10 30 60];            % 信源入射角度
source_number = length(theta); % 信元数
sensor_number = 8;             % 阵元数
w = pi/4;                      % 信号频率
lambda = (2*pi*3e8)/w;         % 信号波长
d = 0.5*lambda;                % 阵元间距
snr = 10;                      % 信噪比为10
n = 500;                       % 快拍数500

A = exp(-1j*2*pi*(0:d:(sensor_number-1)*d).'*sind(theta)/lambda); % 阵列流形
S = randn(source_number,n);                                       % 随机数矩阵
X = A*S;                                                          % 信号矩阵
X1 = awgn(X,snr,'measured');                                      % 添加噪声
Rxx = X1*X1'/n;                                                   % 协方差矩阵

InvS = inv(Rxx);
[EV,D] = eig(Rxx);             % 协方差矩阵对角化 EV特征向量矩阵 D对角阵
EVA = diag(D)';                % 提取对角元素
[EVA,I] = sort(EVA);           % EVA对角元素从小到大排序 I位置索引
EVA = fliplr(EVA);             % 左右翻转
EV = fliplr(EV(:,I));          % 使用索引I对EV列重新排列 左右翻转

for k = 1:361
angle(k)=(k-181)/2;
a = exp(-1j*2*pi*(0:d:(sensor_number-1)*d)*sind((k-181)/2)/lambda).';       % 方向向量
En=EV(:,source_number+1:sensor_number);                                     % 提取source_number+1列到sensor_number列
SP(k)=(a'*a)/(a'*En*En'*a);                                                 % 存储计算值
end
SP=abs(SP);                                                                 % SP的绝对值
SPmax=max(SP);    
SP=10*log10(SP/SPmax);                                                      % 以dB为单位

h=plot(angle,SP);                                                           % 绘制angle与SP的关系图
set(h,'Linewidth',2)                                                        % 绘制线条宽度为 2
xlabel('angle(degree)')                                                     % 为 x 轴和 y 轴添加标签
ylabel('magnituded(dB)')
title('Music Algorithm For Doa', 'fontsize', 16);
set(gca,'XTick',[-90:30:90])                                                % x 轴刻度
grid on

运行结果：

这里可以看到，music算法能够识别多个方向，而且具有很高的分辨率。需要注意的是，当入射信号相干时，music算法会失效，在推导过程中也是以入射信号不相干为前提的。对于相干信号，可采用平滑music算法进行处理，这种方法在之后的文章中进行介绍。

4. 不同信噪比下的RMSE

在信源数为1，阵元数为8，间距 $\frac{\lambda }{2}$，蒙托卡罗次数为1000 时的 RMSE

clear all
close all

source_number = 1;         % 信元数
sensor_number = 8;         % 阵元数
w = pi/4;                  % 信号频率
lambda = (2*pi*3e8)/w;     % 信号波长
d = 0.5*lambda;            % 阵元间距

kps = 100;                 % 快拍数
theta = 30;                % 信源入射角度
mengtimes = 1000;          % 蒙特卡洛次数


A = exp(-1j*2*pi*(0:d:(sensor_number-1)*d).'*sind(theta)/lambda); % 阵列流形
S = randn(source_number,kps);                                     % 随机数矩阵
X = A*S;                                                          % 信号矩阵

snr = 0:1:25;                                                     % 信噪比区间
eror = zeros(length(snr),1);                                      % 估计误差

for n = 1:size(snr,2)
    for m = 1:mengtimes
        X1 = awgn(X,snr(n),'measured'); % 添加噪声
        Rxx = X1*X1'/kps;               % 协方差矩阵
        InvS = inv(Rxx);
        [EV,D] = eig(Rxx);              % 协方差矩阵对角化 EV特征向量矩阵 D对角阵
        EVA = diag(D)';                 % 提取对角元素
        [EVA,I] = sort(EVA);            % EVA对角元素从小到大排序 I位置索引
        EVA = fliplr(EVA);              % 左右翻转
        EV = fliplr(EV(:,I));           % 使用索引I对EV列重新排列 左右翻转
        for k = 1:3601
            a = exp(-1j*2*pi*(0:d:(sensor_number-1)*d)*sind((k-1801)/20)/lambda).';       % 方向向量
            En = EV(:,source_number+1:sensor_number);                                     % 提取source_number+1列到sensor_number列
            SP(k) = (a'*a)/(a'*En*En'*a);                                                 % 存储计算值
        end
        SP = abs(SP);                                                                     % SP的绝对值
        [maxsp,index] = max(SP);                                                          % 搜索谱峰最大值坐标
        estmt(m) = (index-1801)/20;
        eror(n) = (estmt(m)-theta)^2 + eror(n);                                           % 误差平方和
    end
rmse(n) = sqrt(eror(n)/mengtimes);
end

figure
plot(snr,rmse,'Linewidth',2)
grid on
xlabel('SNR(dB)','Fontname','Times New Roman','FontSize',17)
ylabel('\fontname{Times New Roman}\fontsize{17}RMSE°')
legend('\fontname{Times New Roman}\fontsize{17}MUSIC','Location','Best')

运行结果：

参考：
Music算法原理部分参考于：《矩阵分析与应用》张贤达