MUSIC算法仿真

DOA波达方向估计

DOA（Direction Of Arrival）波达方向是指通过阵列信号处理来估计来波的方向，这里的信源可能是多个，角度也有多个。DOA技术主要有ARMA谱分析、最大似然法、熵谱分析法和特征分解法，特征分解法主要有MUSIC算法、ESPRIT算法WSF算法等。其中MUSIC是一种经典的DOA估计算法。

MUSIC算法概述

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

MUSIC算法原理

MUSIC算法是空间谱估计测向理论的重要基石。算法原理 如下：
（1） 不管测向天线阵列形状如何，也不管入射来波入射角的维数如何，如图1所示，假定阵列由$N$个阵元组成，并且符合窄带、远场模型。

图1. 均匀线性阵模型

则阵列输出模型的矩阵形式都可以表示为：$Y(t)=AS+N(t)$
其中，$Y$是观测到的阵列输出数据复向量；$X$是未知的空间信号复向量；$N(t)$是阵列输出向量中的加性噪声；$A$是阵列的方向矩阵，又称之为阵列流形矩阵；当我们以最左侧的阵元为参考位，此处，$A$矩阵表达式表示为：
A = [ α ( θ 1 ) , α ( θ 2 ) , ⋯   , α ( θ M ) ] = [ 1 1 ⋯ 1 e − j ϕ 1 e − j ϕ 2 ⋯ e − j ϕ D e j − 2 ϕ 1 e − j 2 ϕ 2 ⋯ e − j 2 ϕ D ⋮ ⋮ ⋱ ⋮ e − j ( N − 1 ) ϕ 1 e − j ( N − 1 ) ϕ 2 ⋯ e − j ( N − 1 ) ϕ D ] A=[\alpha(\theta_1),\alpha(\theta_2),\cdots,\alpha(\theta_M)]=

$$\begin{equation} \begin{bmatrix} 1 & 1 & \cdots &1 \\ e^{-j\phi_1} & e^{-j\phi_2} & \cdots & e^{-j\phi_D} \\ e^{j-2\phi_1} & e^{-j2\phi_2} & \cdots & e^{-j2\phi_D} \\ \vdots & \vdots & \ddots & \vdots \\ e^{-j(N-1)\phi_1} & e^{-j(N-1)\phi_2} & \cdots & e^{-j(N-1)\phi_D} \end{bmatrix} \end{equation}$$ A=[α(θ1​),α(θ2​),⋯,α(θM​)]= ​1e−jϕ1​ej−2ϕ1​⋮e−j(N−1)ϕ1​​1e−jϕ2​e−j2ϕ2​⋮e−j(N−1)ϕ2​​⋯⋯⋯⋱⋯​1e−jϕD​e−j2ϕD​⋮e−j(N−1)ϕD​​ ​​​
可以看到，$A$是一个范德蒙矩阵其中，${\varphi }_k=2\pi sin{\theta }_k/\lambda$，且${\theta }_k,1\le k\le D$为第$k$个信源的角度，共有$D$个角度。${\varphi }_k$在这里是指波程差引起的相移，如下图所示相邻阵元之间的波程差为$dsin\theta$，将此转换为相位，单位为弧度，则为$2\pi sin{\theta }_k/\lambda$，$\lambda$为信号波长。由下图的等相位面可知，信号到达最左侧的阵元需要走更远的距离，因此，以最左侧的阵元为参考，即左侧阵元接收到的信号为$x_1(t)$，则其他阵元接收到的信号都超前于它，此时，$x_1(t)=s(t),x_2(t)=s(t)e^{j\varphi }$。如果，我们以最右侧信号为参考，即最右侧阵元接收到的信号为$x_1(t)$，则$x_1(t)=s(t),x_2(t)=s(t)e^{-j\varphi }$。注意角度的正负号。

图2. 波程差示意图

接收信号写为向量形式：
Y ( t ) = [ y 1 ( t ) y 2 ( t ) y 3 ( t ) ⋮ y N ( t ) ] Y(t)=

$$\begin{equation} \begin{bmatrix} y_1(t) \\ y_2(t) \\ y_3(t) \\ \vdots \\ y_N(t) \end{bmatrix} \end{equation}$$ Y(t)= ​y1​(t)y2​(t)y3​(t)⋮yN​(t)​ ​​​
$y_i(t)$为第$i$个阵元接收到的信号，共有$N$个阵元。$S=[s_1(t),s_2(t),\cdots ,s_D(t)]$，$s_i(t),1\le i\le D$，表示第$i$个信源发出的信号。$N(t)=[n_1(t),n_2(t),\cdots ,n_N(t){]}^T$，表示噪声，符合独立同分布的假设。
MUSIC算法的处理任务就是设法估计出入射到阵列的空间信号的个数$D$以及空间信号源的强度及其来波方向${\theta }_k$。
（2） 在实际处理中，$Y$得到的数据是有限时间段内的有限次数的样本（也称快拍或快摄），样本的长度取决于采样率和采样时间，在这段时间内，假定来波方向不发生变化，且噪声为与信号不相关的白噪声，则定义阵列输出信号的二阶矩：$R_y$。
$$R_y=E[Y{Y}^H]$$
其中，$[\bullet {]}^H$表示共轭转置，值得注意的是，在MATLAB中我们常用的转置符号${（}^{\prime }）$其实是共轭转置，而普通转置为(.')，当然他们对于实数矩阵而言没有任何区别。$E[\bullet ]$表示统计学上的期望。
由于噪声信号符合独立同分布的假设，信号与噪声互不相关，噪声为零均值白噪声，则上式又可以写为：
$$R_y=E[Y{Y}^H]=E[(AS+N)(AS+N{)}^H]=AE[S{S}^H]A^H+E[N{N}^H]=A{R}_S{A}^H+R_N$$
其中，$R_S=E[S{S}^H]$为信号相关矩阵，$R_N=E[N{N}^H]={\sigma }^{2}I$为噪声相关矩阵，$I$为$N^{\ast }N$的单位矩阵。
（3） MUSIC算法的核心就是对$R_y$进行特征值分解，利用特征向量构建两个正交的子空间，即信号子空间和噪声子空间。对$R_y$进行特征分解。
（4） $U$是非负定的厄米特矩阵，所以特征分解得到的特征值均为非负实数，有$D$个大的特征值和$M-D$个小的特征值，大特征值对应的特征向量组成的空间Us为信号子空间，小特征值对应的特征向量组成的空间$U_n$为噪声子空间。
（5） 将噪声特征向量作为列向量，组成噪声特征矩阵 ，并张成$M-D$维的噪声子空间$U_n$，噪声子空间与信号子空间正交。而$U_s$的列空间向量恰与信号子空间重合，所以$U_s$的列向量与噪声子空间也是正交的，由此，可以构造空间谱函数。
（6） 在空间谱域求取谱函数最大值，其谱峰对应的角度即是来波方向角的估计值。
$$P_{music}=\frac{1}{{\alpha }^H(\theta )U_n{U}_n^H\alpha (\theta )}$$

MUSIC算法MATLB仿真

% MUSIC空间谱估计算法仿真
close all; clear all; clc;
J=sqrt(-1); 
source_number=4; 
source_doa=[30 40 50 70]; 
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,'measured'); 
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');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

仿真结果如下：

图3. 仿真结果

还有一种参考程序：

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;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

仿真结果如下：

图4. 仿真结果