MUSIC算法及MATLAB实现

MUSIC (Multiple Signal Classification)算法，即多信号分类算法，由Schmidt等人于1979年提出。MUSIC算法是一种基于矩阵特征空间分解的方法。从几何角度讲，信号处理的观测空间可以分解为信号子空间和噪声子空间，显然这两个空间是正交的。信号子空间由阵列接收到的数据协方差矩阵中与信号对应的特征向量组成，噪声子空间则由协方差矩阵中所有最小特征值（噪声方差）对应的特征向量组成.

MUSIC算法利用信号子空间和噪声子空间的正交性，构建空间谱函数，通过谱峰搜索，估计信号的参数。对于声源定位来说，需要估计信号的DOA。MUSIC算法对DOA的估计有很高的分辨率.但是该算法对于入射信号的要求非常高，应用该算法的前提是入射信号必须是互不相干的。

假设由 N 个阵元组成的均匀线性阵列，阵列之间的间距为半波长 d = λ / 2 ，假设有M个入射远场信号源，它们的功率相同且互不相干，分别从不同的方向入射到接收阵列，其中入射信号和噪声互不相关。以最左侧的阵元为参考点，各个阵元的位置相比前一个增加距离d。其中一个信号 $S_1(t)$ 的入射方向角为 θ，表示入射信号与线阵法线的夹角。一般假设信号源数M 小于阵元数目 N 。

在第 t次快拍（也就是第t个时刻），得到的接收数据向量为：

$X(t)=A(\theta )S(t)+N(t)$

窄带远场信号的DOA数学模型：

其中：

$X(t)=[X_1(t),X_2(t)...X_N(t) ] ^T$ ， X是一个N*L的矩阵，代表第t时刻下，阵列输出矢量， $X_1(t)$ 是1*L的信号矢量，即该信号有L个采样点

$N=[N_1(t),N_2(t),...N_N(t) ] ^T$ ，噪声

$S=[S_1(t),S_2(t),...S_M(t) ] ^T$ , S是一个M*L的输入信号矩阵，空间中一共有M个入射信号, $S_1(t)$ 是1*L的信号矢量，即该信号有L个采样点

$A=[\alpha(\theta _1),\alpha(\theta _2),...\alpha (\theta _M) ]$ , N*M的空间阵列的导向矢量阵

我们要求解的问题是：通过对采样的输出信号 X(k) ，估计出M个信号源各自的波达方向角 $\theta_i(i=1,2...M)$ 。

在实际处理中，输出信号数据是有限时间段内的有限次数的样本（也称快拍或快摄），在这段时间内，假定来波方向不发生变化，且噪声为与信号不相关的白噪声，则定义阵列输出信号的协方差矩阵： $R_{XX}$ 。

其中 $R_S$ 是输入信号的协方差矩阵。 $R_N$ 是噪声相关阵

MUSIC算法的核心就是对 $R_{XX}$ 进行特征值分解( $R_{XX}$ 是对称矩阵)，利用特征向量构建两个正交的子空间，即信号子空间和噪声子空间(信号和噪声相互独立，数据协方差可分解为信号子空间和噪声子空间)。

其中 $U=\left[ U_S,U_N\right], \Sigma = \text{diag} \left( \sigma_{1}^2,\cdots,\sigma_{N}^2 \right)$ ， $U_S$ 是由大的特征值对应的特征向量张成的子空间，即信号子空间。 $U_N$ 是由小的特征值对应的特征矢量张成的子空间，即为噪声子空间。

由于噪声子空间 $U_N$ 与信号子空间 $U_S$ 正交，则：

则：

上式的充分必要条件： $A^HU_N=\mathbf{0}$

由于实际的干扰因素， $\mathbf{a}^H(\theta)U_N=\mathbf{0}$ 不完全成立，即不完全满足正交性。所以可采用最小优化搜索来估计波达方向：

MUSIC算法实现的步骤

1.根据N个接收信号矢量得到输出信号的协方差矩阵的估计值：

2.对上面得到的协方差矩阵进行特征分解

3.按特征值的大小排序将与信号个数M相等的特征值和对应的特征向量看做信号部分空间，将剩下的特征值和特征向量看做噪声部分空间，得到噪声矩阵

4.由于阵列方向向量与噪声子空间正交的缘故， $\theta _{MUSIC}$ 在信号入射方向上会出现极小值。整理得到空间谱函数，进行极大值的谱峰搜索，使 $\theta$ 变化：

clear; close all;
%%%%%%%% MUSIC for Uniform Linear Array%%%%%%%%
derad = pi/180;      %角度->弧度
N = 8;               % 阵元个数        
M = 3;               % 信源数目
theta = [-30 0 60];  % 待估计角度
snr = 10;            % 信噪比
K = 512;             % 快拍数
 
dd = 0.5;            % 阵元间距 
d=0:dd:(N-1)*dd;
A=exp(-1i*2*pi*d.'*sin(theta*derad));  %方向矢量

%%%%构建信号模型%%%%%
S=randn(M,K);             %信源信号，入射信号
X=A*S;                    %构造接收信号
X1=awgn(X,snr,'measured'); %将白色高斯噪声添加到信号中
% 计算协方差矩阵
Rxx=X1*X1'/K;
% 特征值分解
[EV,D]=eig(Rxx);                   %特征值分解
EVA=diag(D)';                      %将特征值矩阵对角线提取并转为一行
[EVA,I]=sort(EVA);                 %将特征值排序 从小到大
EV=fliplr(EV(:,I));                % 对应特征矢量排序
                 
 
% 遍历每个角度，计算空间谱
for iang = 1:361
    angle(iang)=(iang-181)/2;
    phim=derad*angle(iang);
    a=exp(-1i*2*pi*d*sin(phim)).'; 
    En=EV(:,M+1:N);                   % 取矩阵的第M+1到N列组成噪声子空间
    Pmusic(iang)=1/(a'*En*En'*a);
end
Pmusic=abs(Pmusic);
Pmmax=max(Pmusic)
Pmusic=10*log10(Pmusic/Pmmax);            % 归一化处理
h=plot(angle,Pmusic);
set(h,'Linewidth',2);
xlabel('入射角/(degree)');
ylabel('空间谱/(dB)');
set(gca, 'XTick',[-90:30:90]);
grid on;

from：https://blog.csdn.net/qq_23947237/article/details/82318222

原文链接：https://blog.csdn.net/qq_30815237/article/details/109692896