DOA估计算法
DOA(Direction Of Arrival)波达方向定位技术主要有ARMA谱分析、最大似然法、熵谱分析法和特征分解法,特征分解法主要有MUSIC算法、ESPRIT算法WSF算法等。
MUSIC (Multiple Signal Classification)算法,即多信号分类算法,由Schmidt等人于1979年提出。MUSIC算法是一种基于子空间分解的算法,它利用信号子空间和噪声子空间的正交性,构建空间谱函数,通过谱峰搜索,估计信号的参数。对于声源定位来说,需要估计信号的DOA。MUSIC算法对DOA的估计有很高的分辨率,且对麦克风阵列的形状没有特殊要求,因此应用十分广泛。
运用矩阵的定义,可得到更为简洁的表达式:
X=AS+NX=AS+NX=AS+N
式中
X=[x1(t),x2(t),...xM(t)]TX=[x_1(t),x_2(t),...x_M(t)]^TX=[x1(t),x2(t),...xM(t)]T,
S=[S1(t),S2(t),...SD(t)]TS=[S_1(t),S_2(t),...S_D(t)]^TS=[S1(t),S2(t),...SD(t)]T,
A=[a(θ1),a(θ2),...a(θD)]TA=[a(\theta_1),a(\theta_2),...a(\theta_D)]^TA=[a(θ1),a(θ2),...a(θD)]T,
N=[n1(t),n2(t),...nM(t)]TN=[n_1(t),n_2(t),...n_M(t)]^TN=[n1(t),n2(t),...nM(t)]T。
XXX为阵元的输出,AAA为方向响应向量,SSS是入射信号,NNN表示阵列噪声。
其中 φk=2πdλsinθk\varphi_k=\frac{2\pi d}{\lambda}sin\theta_kφk=λ2πdsinθk有
A=[a(θ1),a(θ2),...a(θD)]T=[11⋯1e−jφ1e−jφ2⋯e−jφD⋮⋮⋱⋮e−j(M−1)φ1e−j(M−1)φ2⋯e−j(M−1)φD]
A=[a(\theta_1),a(\theta_2),...a(\theta_D)]^T=\left[
\begin{matrix}
1 & 1 & \cdots & 1 \\
{e^{-j\varphi_1}} & {e^{-j\varphi_2}} & \cdots & {e^{-j\varphi_D}} \\
\vdots & \vdots & \ddots & \vdots \\
{e^{-j(M-1)\varphi_1}} & {e^{-j(M-1)\varphi_2}} & \cdots & {e^{-j(M-1)\varphi_D}} \\
\end{matrix}
\right]
A=[a(θ1),a(θ2),...a(θD)]T=⎣⎢⎢⎢⎡1e−jφ1⋮e−j(M−1)φ11e−jφ2⋮e−j(M−1)φ2⋯⋯⋱⋯1e−jφD⋮e−j(M−1)φD⎦⎥⎥⎥⎤
对xm(t)x_m(t)xm(t)进行N点采样,要处理的问题就变成了通过输出信号xm(t)x_m(t)xm(t)的采样{xm(i)=1,2,...,M}\{ x_m (i)=1,2,...,M\}{xm(i)=1,2,...,M}估计信号源的波达方向角θ1,θ2...θD\theta_1,\theta_2...\theta_Dθ1,θ2...θD,由此可以很自然的将阵列信号看作是噪声干扰的若干空间谐波的叠加,从而将波达方向估计问题与谱估计联系起来。
对阵列输出X做相关处理,得到其协方差矩阵
Rx=E[XXH]R_x=E[XX^H]Rx=E[XXH]
其中HHH表示矩阵的共轭转置。
根据已假设信号与噪声互不相关、噪声为零均值白噪声,因此可得到:
Rx=E[(AS+N)(AS+N)H]=AE[SSH]AH+E[NNH]=ARSAH+RNR_x=E[(AS+N)(AS+N)^H] =AE[SS^H]A^H+E[NN^H]=AR_SA^H+R_NRx=E[(AS+N)(AS+N)H]=AE[SSH]AH+E[NNH]=ARSAH+RN
其中Rs=E[SSH]R_s=E[SS^H]Rs=E[SSH]称为信号相关矩阵
RN=σ2IR_N=\sigma^2IRN=σ2I是噪声相关阵
σ2\sigma^2σ2是噪声功率
III是M×MM\times MM×M阶的单位矩阵
在实际应用中通常无法直接得到RxR_xRx,能使用的只有样本的协方差矩阵:
Rx^=1N∑i=1NX(i)XH(i)\hat{R_x}=\frac{1}{N} \sum_{i=1}^{N}X(i)X^H (i)Rx^=N1∑i=1NX(i)XH(i),Rx^\hat{R_x}Rx^是RxR_xRx的最大似然估计。
当采样数N→∞N\to\inftyN→∞,他们是一致的,但实际情况将由于样本数有限而造成误差。根据矩阵特征分解的理论,可对阵列协方差矩阵进行特征分解,首先考虑理想情况,即无噪声的情况:Rx=ARsAHR_x=AR_sA^HRx=ARsAH,对均匀线阵,矩阵A由
A=[a(θ1),a(θ2),...a(θD)]T=[11⋯1e−jφ1e−jφ2⋯e−jφD⋮⋮⋱⋮e−j(M−1)φ1e−j(M−1)φ2⋯e−j(M−1)φD]
A=[a(\theta_1),a(\theta_2),...a(\theta_D)]^T=\left[
\begin{matrix}
1 & 1 & \cdots & 1 \\
{e^{-j\varphi_1}} & {e^{-j\varphi_2}} & \cdots & {e^{-j\varphi_D}} \\
\vdots & \vdots & \ddots & \vdots \\
{e^{-j(M-1)\varphi_1}} & {e^{-j(M-1)\varphi_2}} & \cdots & {e^{-j(M-1)\varphi_D}} \\
\end{matrix}
\right]
A=[a(θ1),a(θ2),...a(θD)]T=⎣⎢⎢⎢⎡1e−jφ1⋮e−j(M−1)φ11e−jφ2⋮e−j(M−1)φ2⋯⋯⋱⋯1e−jφD⋮e−j(M−1)φD⎦⎥⎥⎥⎤
所定义的范德蒙德矩阵,只要满足θi≠θj,i≠j\theta_i\neq \theta_j,i\neq jθi=θj,i=j,则他的各列相互独立。
若RsR_sRs为非奇异矩阵Rank(Rs)=DRank(R_s)=DRank(Rs)=D,各信号源两两不相干,且M>DM>DM>D,则rand(ARsAH)=Drand(AR_sA^H)=Drand(ARsAH)=D,
由于Rx=E[XXH]R_x=E[XX^H]Rx=E[XXH],有:
RsH=RxR_s^H=R_xRsH=Rx
即RsR_sRs为Hermite矩阵,它的特性是都是实数,又由于RsR_sRs为正定的,因此ARsA……HAR_sA……HARsA……H为半正定的,它有D个正特征值和M−DM-DM−D个零特征值。
再考虑有噪声存在的情况
Rx=ARsAH+σ2IR_x=AR_sA^H+\sigma^2IRx=ARsAH+σ2I
由于σ2>0\sigma^2>0σ2>0,RxR_xRx为满秩阵,所以RxR_xRx有M个正实特征值λ1,λ2...λM\lambda_1,\lambda_2...\lambda_Mλ1,λ2...λM
分别对应于M个特征向量v1,v2...vMv_1,v_2...v_Mv1,v2...vM。又由于RxR_xRx为Hermite矩阵,所以各特征向量是正交的,即:viHvj=0,i≠jv_i^Hv_j=0,i\neq jviHvj=0,i=j与信号有关的特征值只有D个,分别等于矩阵ARsAHAR_sA^HARsAH的各特征值与σ2\sigma^2σ2之和,其余M−DM-DM−D个特征值为σ2\sigma^2σ2,即σ2\sigma^2σ2为RRR的最小特征值,它是M−DM-DM−D维的,对应的特征向量vi,i=1,2,...,Mv_i,i=1,2,...,Mvi,i=1,2,...,M中,也有D个是与信号有关的,另外M−DM-DM−D个是与噪声有关的,可利用特征分解的性质求出信号源的波达方向θk\theta_kθk。
MUSIC算法的原理及实现
通过对协方差矩阵的特征值分解,可得到如下结论:
将矩阵RxR_xRx的特征值进行从小到大的排序,即λ1≥λ2≥...≥λM>0\lambda_1 \geq \lambda_2\geq...\geq\lambda_M>0λ1≥λ2≥...≥λM>0,其中D个较大的特征值对应于信号,M−DM-DM−D个较小的特征值对应于噪声。
矩阵RxR_xRx的属于这些特征值的特征向量也分别对应于各个信号和噪声,因此可把RxR_xRx的特征值(特征向量)划分为信号特征(特征向量)与噪声特征(特征向量)。
设λi\lambda_iλi为RxR_xRx的第iii个特征值,viv_ivi是与λi\lambda_iλi个相对应的特征向量,有:
Rxvi=λiviR_xv_i=\lambda_iv_iRxvi=λivi
再设λi=σ2\lambda_i=\sigma^2λi=σ2是RxR_xRx的最小特征值Rxvi=σ2vii=D+1,D+2...MR_xv_i=\sigma^2v_i i=D+1,D+2...MRxvi=σ2vii=D+1,D+2...M,
将Rx=ARsAH+σ2IR_x=AR_sA^H+\sigma^2IRx=ARsAH+σ2I代入可得σ2vi=(ARsAH+σ2I)vi\sigma^2v_i=(AR_sA^H+\sigma^2I)v_iσ2vi=(ARsAH+σ2I)vi,
将其右边展开与左边比较得:
ARsAHvi=0AR_sA^Hv_i=0ARsAHvi=0
因AHAA^HAAHA是D∗DD*DD∗D维的满秩矩阵,(AHA)−1(A^HA)^{-1}(AHA)−1存在;
而Rs−1R_s^{-1}Rs−1同样存在,则上式两边同乘以Rs−1(AHA)−1AHR_s^{-1}(A^HA)^{-1}A^HRs−1(AHA)−1AH,
有:
Rs−1(AHA)−1AHARsAHvi=0R_s^{-1}(A^HA)^{-1}A^HAR_sA^Hv_i=0Rs−1(AHA)−1AHARsAHvi=0
于是有
AHvi=0,i=D+1,D+2,...,MA^Hv_i=0,i=D+1,D+2,...,MAHvi=0,i=D+1,D+2,...,M
上式表明:噪声特征值所对应的特征向量(称为噪声特征向量)viv_ivi,与矩阵AAA的列向量正交,而AAA的各列是与信号源的方向相对应的,这就是利用噪声特征向量求解信号源方向的出发点。
用各噪声特征向量为例,构造一个噪声矩阵EnE_nEn:
En=[vD+1,vD+2,...vM]E_n=[v_{D+1},v_{D+2},...v_{M}]En=[vD+1,vD+2,...vM]
定义空间谱Pmu(θ)P_{mu}(\theta)Pmu(θ):
Pmu(θ)=1aH(θ)EnEnHa(θ)=1∥EnHa(θ)∥2P_{mu}(\theta)=\frac{1}{{a^H}(\theta)}E_nE_n^Ha(\theta)=\frac{1}{\Vert E_n^Ha(\theta)\Vert^2}Pmu(θ)=aH(θ)1EnEnHa(θ)=∥EnHa(θ)∥21
该式中分母是信号向量和噪声矩阵的内积,当a(θ)a(\theta)a(θ)和EnE_nEn的各列正交时,该分母为零,但由于噪声的存在,它实际上为一最小值,因此Pmu(θ)P_{mu}(\theta)Pmu(θ)有一尖峰值,由该式,使θ\thetaθ变化,通过寻找波峰来估计到达角。
MUSIC算法实现的步骤
1.根据N个接收信号矢量得到下面协方差矩阵的估计值:
Rx=1N∑i=1NX(i)XH(i)R_x=\frac{1}{N} \sum_{i=1}^NX(i)X^H(i)Rx=N1∑i=1NX(i)XH(i)
对上面得到的协方差矩阵进行特征分解
Rx=ARsAH+σ2IR_x=AR_sA^H+\sigma^2IRx=ARsAH+σ2I
2.按特征值的大小排序将与信号个数DDD相等的特征值和对应的特征向量看做信号部分空间,将剩下的M−DM-DM−D个特征值和特征向量看做噪声部分空间,得到噪声矩阵EnE_nEn:
AHvi=0,i=D+1,D+2,...,MA^Hv_i=0,i=D+1,D+2,...,MAHvi=0,i=D+1,D+2,...,M
En=[vD+1,vD+2,...vM]E_n=[v_{D+1},v_{D+2},...v_{M}]En=[vD+1,vD+2,...vM]
3.使θ\thetaθ变化,按照式
Pmu(θ)=1aH(θ)EnEnHa(θ)P_{mu}(\theta)=\frac{1}{{a^H}(\theta)E_nE_n^Ha(\theta)}Pmu(θ)=aH(θ)EnEnHa(θ)1
来计算谱函数,通过寻求峰值来得到波达方向的估计值。
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;