MUSIC算法原理及MATLAB代码阵列信号处理

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

MUSIC算法是空间谱估计测向理论的重要基石。算法原理 [2] 如下：
（1）不管测向天线阵列形状如何，也不管入射来波入射角的维数如何，假定阵列由M个阵元组成，则阵列输出模型的矩阵形式都可以表示为：Y(t)=AX(t)+N(t)
其中，Y是观测到的阵列输出数据复向量；X是未知的空间信号复向量；N是阵列输出向量中的加性噪声；A是阵列的方向矩阵；此处，A矩阵表达式由图册表示。
MUSIC算法的处理任务就是设法估计出入射到阵列的空间信号的个数D以及空间信号源的强度及其来波方向。
（2）在实际处理中，Y得到的数据是有限时间段内的有限次数的样本（也称快拍或快摄），在这段时间内，假定来波方向不发生变化，且噪声为与信号不相关的白噪声，则定义阵列输出信号的二阶矩：Ry。

（3） MUSIC算法的核心就是对Ry进行特征值分解，利用特征向量构建两个正交的子空间，即信号子空间和噪声子空间。对Ry进行特征分解，即是使得图册中的公式成立。
（4） U是非负定的厄米特矩阵，所以特征分解得到的特征值均为非负实数，有D个大的特征值和M-D个小的特征值，大特征值对应的特征向量组成的空间Us为信号子空间，小特征值对应的特征向量组成的空间Un为噪声子空间。
（5）将噪声特征向量作为列向量，组成噪声特征矩阵，并张成M-D维的噪声子空间Un，噪声子空间与信号子空间正交。而Us的列空间向量恰与信号子空间重合，所以Us的列向量与噪声子空间也是正交的，由此，可以构造空间谱函数。

（6）在空间谱域求取谱函数最大值，其谱峰对应的角度即是来波方向角的估计值。

MUSIC（Multiple Signal Classification多信号分类）算法是1979年由美国人R.O.Schmidt提出的，它标志着空间谱估计测向进入了繁荣发展的阶段。它将“向量空间”的概念引入了空间谱估计领域，经过三十年的发展，可以说其理论已经比较成熟。
自80年代以来，人们对基于特征分解的超分辨率空间谱估计算法进行了广泛深入的研究，并提出了一系列高效的处理方法，其中最经典的是多信号分类（MUSIC）算法，这种算法要经过一维搜索才能求出信源的来向，而相对最大似然（ML）和加权子空间拟合（WSF）等多维搜索算法的运算量已经减少了很多。以MUSIC为代表的算法存在一个缺点，即对相干信号处理的不理想。在针对相干信号源的一系列处理方案中，比较经典的是空间平滑技术，如空间平滑（SS）和修正的空间平滑（MSS）算法。然而，空间平滑技术是以损失阵列有效孔径为代价的，而且只适用于等距均匀线阵（ULA）。
事实上空间谱估计算法都是在已知信号源数目下计算的，而在实际应用中这是不可能的，只能根据观测数据对源数目进行估计。R.O.Schmidt在他的经典之作中提出了依据阵列协方差矩阵特征值的分布来估计信号源的方法。这种方法在理论上是完美的，至少对独立源和部分相关源是正确的，但实际上由于数据长度有限，很大程度上只能依靠主观判断来确定源数。

MUSIC算法的分辨率要优于传统波束形成算法（空间傅里叶变换算法），有时候比capon算法分辨率也要高

J=sqrt(-1);
source_number=4;
source_doa=[30 45 60 135];
sensor_number=7;
snapshot_number=2000;
snr=10; A=exp(-J*(0:sensor_number-1)'*pi*sin(source_doa*pi/180));
s=(randn(source_number,snapshot_number)+J*randn(source_number,snapshot_number))/sqrt(2);
x=A*s;
y=awgn(x,snr);
R=y*y'/snapshot_number; [V,D]=eig(R);
Un=V(:,1:sensor_number-source_number);
Gn=Un*Un'; searching_doa=0:0.1:90;
for i=1:length(searching_doa) a_theta=exp(-J*(0:sensor_number-1)'*pi*sin(pi*searching_doa(i)/180))P_con(i)=abs(a_theta'*R*a_theta);P_BF(i)=abs((a_theta'*R*a_theta)./(a_theta'*a_theta)); P_capon(i)=1./abs((a_theta'*inv(R)*a_theta)); P_music(i)=1./abs((a_theta'*Gn*a_theta));
end
plot(searching_doa,P_con/max(P_con),'k');hold on;
plot(searching_doa,P_BF/max(P_BF),'r'); hold on;
plot(searching_doa,P_capon/max(P_capon),'g'); hold on;
plot(searching_doa,P_music/max(P_music),'b'); hold off;grid on;
xlabel('ang');
ylabel('功率谱估计');
legend('conditional spectrum','Bartlett spectrum','Capon spectrum','Music spectrum');

人人为我，我为人人
记得点个赞啊