DOA估计中的Capon、MUSIC算法

最近也在上阵列信号处理的课程，目前学习的主要内容包括阵列的信号模型、方向图、空时等效、Capon波束形成器、MMSE波束形成器、LCMV波束形成器、GSC波束形成器、传统测向法(比相法、CBF/Bartlett估计器)、MUSIC算法、ESPRIT算法以及空间平滑法。下面主要介绍DOA估计算法，主要内容有：MUSIC算法及改进、MUSIC及Capon法的比较、LS-ESPRIT以及TLS-ESPRIT。本文将介绍前两部分。

DOA估计

时域频谱表示信号在各个频率上的能量分布；空间谱表示信号在空间各个方向上的能量分布。所以如果能够得到信号的空间谱，就能够得到信号的波达方向 (Direction of arrival, DOA) 所以空间谱估计也被称为DOA估计。

1. DOA估计的分类主要有多信号源DOA、相干源DOA、宽带信号源DOA以及复杂环境下的DOA…
2. DOA估计的方法主要包括传统法、子空间法、最大似然法和综合法。传统法存在瑞利限；最大似然法高分辨最优但计算量大；子空间法属于次最优算法；综合法也会用到子空间法。过多的就不再展开赘述了，感兴趣者可自行查阅文献。

瑞利限: 一般的，在一定阵列长度下可达到的最小分辨率称为瑞利限，超过瑞利限的方法被称为 超分辨方法。

• 1979年，Schmidt提出MUSIC算法，1986年重新发表
• 1986年，Roy提出ESPRIT算法

以上两种方法属于早期经典的超分辨方法，今天大多数子空间法都是在这两种方法基础上改进的。

MUSIC算法

算法假设：

1.均匀线阵且阵元间距不大于处理最高频率信号波长的二分之一
2.处理器噪声为加性高斯分布，不同阵元间距噪声均为平稳随机过程，独立同分布，各阵元噪声方差相同
3.信号为零均值平稳随机过程且与噪声独立
4.阵元数大于信源数，信号源为窄带信号且入射信源之间弱相关或者不相关。
##窄带：信号通过天线阵列的时间远远小于信号带宽的倒数

MUSIC算法步骤

1.设置信源、阵列模型、接收信号
2.计算接收信号的协方差矩阵
3.协方差矩阵特征分解，得到特征值
4.特征值由大到小排序，求对应的信号子空间和噪声子空间
5.利用计算MUSIC谱，即噪声功率谱平方和的倒数
6.搜索MUSIC谱的峰值，找到最大的K个谱峰，对应的K个DOA估计的角度

• 计算MUSIC谱的公式如下

$$P_{music}=\frac{1}{\sum _{i=P+1}^N{\left|v_i^{H}a(\theta )\right|}^2}$$，其中N为阵元数，P为信号子空间数，$v_i$为特征向量。

• Capon方法的计算公式如下
  $$P_{capon}=\frac{1}{a^H(\theta )R_x^{-1}a(\theta )}$$
  其中$R_x^{-1}$为接收信号协方差的逆矩阵。

下面通过仿真来对比一下MUSIC与Capon的区别，参考的文章中给出了不同信噪比下以及不同角度间隔下两种方法的估计误差。
1. 信噪比
设置阵元数为10，阵元间隔为半波长，信源数为3(-10度,0度,20度)，快拍数为1024，下图为估计得到的信号谱，左侧信噪比设置为-8dB，右侧信噪比设置为10dB。

从上图可以看出，在信噪比较低时两种估计算法性能都急剧下降，但是MUSIC算法略优于Capon，而当信噪比较大时，两种算法基本一致。

关于谱峰强度的解释：

 MUSIC谱峰只反映阵列流形矢量与噪声子空间的正交性，与信噪比无关；
 Capon谱峰是真正的输出功率，与信噪比有关。

2.角度间隔/分辨率
设置阵元数为10，阵元间隔为半波长，信源数为3(-0.5度,0度,0.5度)，快拍数为1024，信噪比设置为10dB，下图为估计得到的信号谱，为方便观察进行了归一化。
可以看到这种情况下，MUSIC的分辨率是优于Capon法的。
上述仿真分析感觉也不是特别严谨，至于两种方法谁好谁坏，在处理具体问题的时候都试一下哪个好就用哪个吧😄😄
下面给出仿真代码

clear all;
close all
clc;
source_number=3;%信元数
sensor_number=10;%阵元数
N_x=1024; %信号长度
snapshot_number=N_x;%快拍数
w=[pi/4 pi/6 pi/3].';%信号频率  三个非相干元
l=sum(2*pi*3e8./w)/3;%信号波长  
d=0.5*l;%阵元间距
snr=10;%信噪比
source_doa=[-0.5 0 0.5];%两个信号的入射角度
A=[exp(-1j*(0:sensor_number-1)*d*2*pi*sin(source_doa(1)*pi/180)/l);exp(-1j*(0:sensor_number-1)*d*2*pi*sin(source_doa(2)*pi/180)/l);exp(-1j*(0:sensor_number-1)*d*2*pi*sin(source_doa(3)*pi/180)/l)].';%阵列流型
s=sqrt(10.^(snr/10))*exp(1j*w*[0:N_x-1]);%仿真信号
x=awgn(A*s,snr);
R=x*x'/snapshot_number;
inv_R=inv(R);
[V,D]=eig(R);
Un=V(:,1:sensor_number-source_number);%信源已知，直接划分
Gn=Un*Un';

searching_doa=-90:0.1:90;%线阵的搜索范围为-90~90度
 for i=1:length(searching_doa)
   a_theta=exp(-1j*(0:sensor_number-1)'*2*pi*d*sin(pi*searching_doa(i)/180)/l);
   Pmusic(i)=a_theta'*a_theta./abs((a_theta)'*Gn*a_theta);
   Pcapon(i)=1./abs((a_theta)'*inv_R*a_theta);
 end
plot(searching_doa,10*log10(Pmusic/max(Pmusic)),'k-',searching_doa,10*log10(Pcapon/max(Pcapon)),'b--');
xlabel('DOAs/degree');
xlim([-5 5])
ylabel('Normalized Spectrum/dB');
legend('Music Spectrum','Capon Spectrum');
title('Comparation of MUSIC and Capon for DOA Estimation');
grid on;

从代码上看还是很简单的，可以通过修改线阵搜索的步进值实现分辨率一定程度上地提高。

MUSIC算法的改进

1. 信源数未知——>一维噪声子空间法(加权子空间法的特例)
2. 色噪声——>基于高阶统计量的阵列处理
3. 快拍数有限——>求根MUSIC
4. 信噪比低——>加权MUSIC
5. 相干源
...

相干源下的MUSIC算法
对于相干源，采取的措施是将信号取共轭重排之后计算自相关后与原信号自相关相加，再进行MUSIC算法。
仿真参数与上文一致，只是信号频率变为[pi/4, pi/6, pi/4],入射角度为[-10, 0, 10],结果如下图所示

从上图可以看到，MUSIC算无法估计出正负10度的两个相干信号，而相干的MUSIC算法实现了精确的估计。
下面给出仿真代码

clc
clear all
close all

source_number=3;%信元数
sensor_number=10;%阵元数
N_x=1024; %信号长度
snapshot_number=N_x;%快拍数
w=[pi/4 pi/6 pi/4].';%信号频率  三个非相干元
l=sum(2*pi*3e8./w)/3;%信号波长  
d=0.5*l;%阵元间距
snr=10;%信噪比

source_doa=[-10 0 10];%两个信号的入射角度
A=[exp(-1j*(0:sensor_number-1)*d*2*pi*sin(source_doa(1)*pi/180)/l);exp(-1j*(0:sensor_number-1)*d*2*pi*sin(source_doa(2)*pi/180)/l);exp(-1j*(0:sensor_number-1)*d*2*pi*sin(source_doa(3)*pi/180)/l)].';%阵列流型

s=sqrt(10.^(snr/10))*exp(1j*w*[0:N_x-1]);%仿真信号
x=awgn(A*s,snr);
% x=A*s+(1/sqrt(2))*(randn(sensor_number,N_x)+1j*randn(sensor_number,N_x));%加了高斯白噪声后的阵列接收信号

R=x*x'/snapshot_number;
B=eye(sensor_number);
J=fliplr(B);
Y=J*conj(x);
R1=Y*Y'/snapshot_number;
R2=R1+R;
%相干music
[V,D]=eig(R2);
[lamda,index]=sort((diag(D)));
UU=V(:,index(1:source_number)); %噪声子空间
Gn_music=UU*UU';
%music
[V_m,D_m]=eig(R);
[lamda_m,index_m]=sort((diag(D_m)));
Un=V_m(:,index(1:source_number));
Gn_mmusic=Un*Un';

searching_doa=-90:0.1:90;%线阵的搜索范围为-90~90度
 for i=1:length(searching_doa)
   a_theta=exp(-1j*(0:sensor_number-1)'*2*pi*d*sin(pi*searching_doa(i)/180)/l);
   Pmusic(i)=a_theta'*a_theta./abs((a_theta)'*Gn_music*a_theta);
   Pmmusic(i)=a_theta'*a_theta./abs((a_theta)'*Gn_mmusic*a_theta);
 end

plot(searching_doa,10*log10(Pmusic/max(Pmusic)),'k-',searching_doa,10*log10(Pmmusic/max(Pmmusic)),'b--');
xlabel('DOAs/degree');
% xlim([-5 5])
ylabel('Normalized Spectrum/dB');
legend('相干Music','Music');
title('Comparation of MUSIC and Capon for DOA Estimation');
grid on;

#-----------------------------------------------------------------------------总结---------------------------------------------------------------------------------#
本文主要进行了较为基础的MUSIC算法的实现以及与Capon方法的对比，涉及到的相关代码均已给出，供大家参考。
#-------------------------------------------------------------------------下一步计划----------------------------------------------------------------------------#

• ESPRIT算法的仿真与实现
• LS_ESPRIT算法
• TLS_ESPRIT算法

#--------------------------------------------------------------------------参考文献------------------------------------------------------------------------------#
《现代数字信号处理》——王展，李双勋等
《最优阵列处理技术》——Harry L. Van Trees…
《The Difference Between Capon and MUSIC Algorithms》 Huiping Huang