Traditional Comm笔记【5】：基于MUSIC Algorithm的DoA/AoA估计以及MATLAB实现

1. DoA/AoA

"AoA"：Angle of Arrival，波达角。
"DoA"：Direction of Arrival，波达角。
波达方向是指空间信号的到达方向（各个信号到达阵列参考阵元的方向角，简称波达方向），“波达角（AoA）与波达方向（DoA）是同一回事”。

“波达方向（Direction Of Arrival）估计”，又称为“角谱估计（Angle spectral estimation）”、“波达角（Angle Of Arrival）估计”。一个信源有很多可能的传播路径和到达角。如果几个发射机同时工作，每个信源在接收机处形成潜在的多径分量。因此，接收天线能估计出这些到达角就显得很重要，目的是估计出哪个发射机在工作以及发射机所处的方向，简单的说就是利用己方雷达接收来自目标发射机的来波方向进行估计；其物理原理是利用电磁波的直线传播原理。

对于 $3D$ 空间来说，$DoA/AoA$包括 $(\theta ,\varphi )$。对于 $2D$ 空间来说，$DoA$ 即为方位角 $\theta$。

figure1 DoA/AoA示意图

2. 导向向量steering vector

1. 导向矢量的本质是描述空间相位差的；
2. 导向向量的结构和阵元之间的相对位置有关系。$ULA$ 的特点是 $Vandermonde$ 结构。对于其它集合形状的阵列，导向向量的结构将会有不同。
3. 导向向量的值，是来波方向的函数，若 $\theta$ 不同，则导向矢量的值会有所不同。
4. 对于同一方向 $(e.g.30^{°})$ ，若选取的参考点不同，那么导向矢量的值也会不同。但是阵元之间的相对相位差不会变化。

3. MUSIC Algorithm

(1). MUSIC算法概述：

“MUSIC（Multiple Signal Classification，多重信号分类）”，是一类“空间谱估计算法”。其思想是利用接收数据的协方差矩阵进行特征分解，分离出信号子空间和噪声子空间，利用信号“方向向量/导向向量（steering vector）”与噪声子空间的正交性来构成空间扫描谱，进行全域搜索谱峰，从而实现信号的参数估计。

在众多性能优良的高分辨率 $DoA$ 估计算法中， $MUSIC$ 算法最经典，它在空域内进行谱峰搜索求出信源来向。与最大似然、加权子空间拟合等多维搜索算法相比，$MUSIC$ 算法运算量要小很多。在$MUSIC$ 算法的基础上，发展出了加权$MUSIC$ 和改进$MUSIC$ 算法等。

$MUSIC$ 算法建立在以下基础之上：

1. 阵列形式为线性均匀阵，阵元间距不大于处理最高频率信号波长的二分之一。
2. 信号源个数小于阵元数目，以确保阵列流形矩阵的各个列线性独立。若传感器的数量比信源的个数多，则阵列数据的信号分量一定位于一个低秩的子空间，在一定条件下，这个子空间的将唯一确定信号的波达方向，并且可以使用数值稳定的奇异值分解精确的确定波达方向。
3. 处理器的噪声为加性高斯分布，不同阵元间距噪声均为平稳随机过程，各阵元间噪声相互独立，空间平稳（各噪声方差相等）。
4. 空间信号为零均值平稳随机过程，信号与阵元噪声相互独立。
5. 信号源通常为窄带远场信号。

正是由于$MUSIC$ 算法在特定的条件下具有很高的分辨力、估计精度及稳定性，从而吸引了大量的学者对其进行深入的研究和分析。

(2). MUSIC算法原理：

将任意阵列输出数据的协方差矩阵进行特征值分解，对应不同特征值的特征向量构成相互正交的信号子空间和噪声子空间，大特征值对应的特征向量构成的是信号子空间，小特征值对应的特征向量构成的是噪声子空间，然后利用两个子空间之间的正交性来估计信号的方向。

假设 $N$ 元等距线阵，阵元间距为 $d$，信号的波长为 $\lambda$，空间中有 $N$ 个信源，那么接收到的观测信号为：
$$X(t)=\sum _{i=1}^N{S}_i(t)A({\theta }_i)+N(t)$$
其中，$S_i(t)$ 是第 $i$ 个信号，${\theta }_i$ 是空间中第 $i$ 个信号的入射角度。这里假设线性阵列总共接收到 $N$ 个同源且不相干的信号。

figure2 均匀线阵模型/等距线阵模型

将最左边的天线信道响应作为接受基准，假设为 $1$。
所以有：
$$\begin{gathered} X(t)=A\cdot S(t)+N(t) \\ A=[A({\theta }_1),A({\theta }_2),...,A({\theta }_N)] \\ S(t)=[S_1(t),S_2(t),...,S_N(t){]}^T \end{gathered}$$

在基于天线阵列协方差矩阵的特征分解类 $DoA$ 估计算法中，$MUSIC$ 算法具有普遍的适用性，只要已知天线阵的布阵形式，无论是线阵还是圆阵，不管阵元是否等间隔分布，都可以得到高分辨率的估计结果。阵列协方差矩阵 $R$ 可以划分为两个空间，即 $R=U_S{\Sigma }_S{U}_S^H+U_N{\Sigma }_N{U}_N^H$ 。可得
$$R{U}_N=A(\theta )R_S{A}^H(\theta )U_N+{\sigma }_n^2{U}_N={\sigma }_n^2{U}_N(1)$$
根据式 $(1)$可得
$$A(\theta )R_S{A}^H(\theta )U_N=0(2)$$
矩阵 $R_S$ 为满秩阵，非奇异，所以有逆存在。于是，上式可变为 $A^H(\theta )U_N=0$ ，这说明矩阵 $A(\theta )$ 中的各个列向量与噪声子空间正交，故有
$$U_N^{H}a({\theta }_i)=0,i=1,2,...,K(3)$$
由噪声特征向量和信号向量的正交关系，得到阵列空间谱函数
$$P_{MUSIC}(\theta )=\frac{1}{a^H(\theta )U_N{U}_N^{H}a(\theta )}(4)$$
由上式，使 $\theta$ 变化，通过寻找波峰来估计到达角。则 $(4)$ 式还可表示为
θ i = a r g θ   m i n   a H ( θ ) U N U N H a ( θ ) = a r g θ   m i n   t r { P a U N U N H }       ( 5 ) \theta_i=arg_{\theta}\ min\ a^H(\theta)U_NU_N^Ha(\theta)=arg_{\theta}\ min\ tr\{P_aU_NU_N^H\}\ \ \ \ \ (5) θi​=argθ​ min aH(θ)UN​UNH​a(θ)=argθ​ min tr{Pa​UN​UNH​}     (5)
其中， $P_a=a(\theta )[a^H(\theta )a(\theta ){]}^{-1}{a}^H(\theta )$ 为 $a(\theta )$ 的投影矩阵。

(3). 基本MUSIC算法步骤：

根据以上讨论，将 $MUSIC$ 算法的步骤归纳如下：

• $step1.$ 根据 $N$ 个接收信号向量得到协方差矩阵的估计值
  $$R=\frac{1}{N}{\Sigma }_{n=1}^{N}x(n)x^H(n)$$
• $step2.$ 对式 $(6)$ 得到的协方差矩阵进行特征值分解 $R=U\Sigma U^H$。
• $step3.$ 按特征值的大小顺序，把与信号个数 $K$ 相等的最大特征值对应的特征向量看成信号子空间，把剩下的 $M-K$ 个特征值对应的特征向量看成噪声子空间，则 $R=U_S{\Sigma }_S{U}_S^H+U_N{\Sigma }_N{U}_N^H$。
• $step4.$ 使 $\theta$ 变化，按照 $P_{MUSIC}(\theta )=1/[a^H(\theta )U_N{U}_N^{H}a(\theta )]$ 来计算谱函数，通过寻求峰值来得到波达方向的估计值。对信号 $AoA$ 进行遍历搜索时，当搜索到的角度与 $N$ 个信号中某个信号的 $AoA$ 相同时，空间谱函数就会产生一个极值点，因此遍历完所有的角度之后，就能产生 $N$ 个极值点。

(4). 补充：

1. 在信噪比一定的情况下，随着阵元数量的减少算法的分辨能力以及精度都在下降，当阵元数量小于或者等于信源个数时，算法性能急剧下降，甚至无法工作。这是由于当阵元数量小于信源个数时，噪声子空间与信号子空间已经不再正交，因为利用 $MUSIC$ 算法估计多个信源角度时，阵元数一定要大于信源个数，而且阵元数越多角度的分辨率越高。
2. $MUSIC$ 算法基本核心思想是子空间正交性，即接收数据协方差可以分解为信号子空间和噪声子空间，这两个空间正交，又因为信号子空间和阵列流形空间相同，所以把协方差矩阵特征分解得到噪声子空间的特征向量，与阵列流形一起构成了空间频谱。

4. MATLAB实现

(1). 通过steering vector构造接收信号：

按照基本 $MUSIC$ 算法的步骤得到MATLAB代码：$MUSIC.m$

clear all
close all
derad = pi/180;   %角度->弧度
radeg = 180/pi;   %弧度->角度
twpi = 2*pi;
kelm = 8;         %阵元个数     
dd = 0.5;         %阵元间距       
d=0:dd:(kelm-1)*dd;     
iwave = 3;        %信源数             
 
%波达方向/待估计角度为15 28 60
%1.给出待估计角度，接收的数据直接为X(t)，公式中是将接收数据分解为A*S+N，所以才需要生成导向矩阵。
%2.仿真中需要生成导向矢量是因为我们需要模拟天线接收到的信号，所以会给出角度。
%3.实际应用中是不知道这个角度的，所以不需要生成导向矩阵A，导向矢量是接收天线本身的属性。
theta = [15 28 60]; 

snr = 10;                                %信噪比
n = 500;                                 %快拍数或者称为采样数
A=exp(-1i*twpi*d.'*sin(theta*derad));    %方向向量,.'转置，'共轭转置
S=randn(iwave,n);                        %信源信号，正态分布随机矩阵
X=A*S;                                   %接收信号
X1=awgn(X,snr,'measured');               %添加高斯白噪声
Rxx=X1*X1'/n;                            %计算协方差矩阵
[EV,D]=eig(Rxx);     %计算Rxx的特征值对应的对角阵D和特征向量构成的矩阵EV
EVA=diag(D)';        %diag抽取矩阵对角线元素
[EVA,I]=sort(EVA);   %特征值从小到大排序，sort只能从小到大排列
EVA=fliplr(EVA);     %特征值左右翻转，从大到小序
EV=fliplr(EV(:,I));  %对应特征向量排序

%构造MUSIC谱函数
for iang = 1:361
        angle(iang)=(iang-181)/2;
        phim=derad*angle(iang);
        a=exp(-1i*twpi*d*sin(phim)).';
        L=iwave;    
        En=EV(:,L+1:kelm);               %得到噪声子空间
        SP(iang)=(a'*a)/(a'*En*En'*a);
end

%作图
SP=abs(SP);
SPmax=max(SP);
SP=10*log10(SP/SPmax);
h=plot(angle,SP);
set(h,'Linewidth',2)
xlabel('angle (degree)')
ylabel('magnitude (dB)')
axis([-90 90 -60 0])
set(gca, 'XTick',(-90:10:90))
grid on  

figure3 基于MUSIC的DoA/AoA估计结果

(2). 直接读取接收信号：

所用数据为：data1.txt，提取码：pyoo；读取数据之后仍然按照基本 $MUSIC$ 算法的步骤得到MATLAB代码：$MusicData1.m$，该资源上传到：基于MUSIC算法的DoA/AoA估计MATLAB实现.

5.Reference

[1] 阵列信号处理及MATLAB实现（第2版）@张小飞 李建峰 徐大专 等 著.
[2] 百度百科：DOA（波达方向定位技术）.
[3] 【算法】 music算法的总结.