Music 算法仿真实验报告

Music算法仿真实验报告一，实验原理θ取第一个阵元为参考点，设参考点接收的入射波信号为：x1(t)=s(t)ejwt其中s(t)为信号的复振幅，w为信号的角频率阵元m接收到的信号为：xm(t)=s[t−τ(θ)]ejw[t−τm(θ)]τ(θ)为阵元m接收到的入射波相对于参考点入射波的延时，θ为入射角。对于窄带信号，相对于τ(θ)时间延迟的信号复包络变化可以忽略，则：s[t−τ(θ)]≈s(t)所以阵元m接收到的信号可简化为：xm(t)=s(t)ejw(t−τm)=x1(t)e−jwτm将M个阵元接收到的信号表示成矢量的形式，可以得到阵列接收信号矢量X为：X=[x1(t)⋮xM(t)]=x1(t)[1⋮e−jwτM(θ)]=x1(t)a(θ)a(θ)=[1,…,e−jwτM(θ)]为阵列导向矢量。如果入射到天线阵列的信号为D个远场信号，并且入射信号的入射角分别为(θ1，θ2，…θD),此时阵列接收信号矢量表示式为：X(t)=∑si(t)a(θi)=AS(t)Di=1其中，A=[a(θ1),a(θ2),…a(θD)]为阵列对信号的导向矢量矩阵。S(t)=[s1(t),s2(t),⋯,sM(t)]T,为信号矢量。考虑到各个振元存在噪声，阵列接收矢量可表示为：X(t)=AS(t)+N(t)N(t)为阵列接收的噪声矢量，表示为：N(t)=[n1(t),n2(t),⋯nM(t)]T而且ni(t)为零均值，方差为σ2的相互独立的白噪声采样对阵列输出X作相关处理，得到其协方差矩阵RxRx=E[XXH]假设信号与噪声互不相关，且噪声为零均值白噪声，因此可以得到：RX=E[(AS+N)(AS+N)H]=AE[SSH]AH+E[NNH]=ARSAH+RNRN=σ2IRN是噪声的相关矩阵，σ2是噪声功率，I是单位矩阵当所有信号互不相关是，有E[si(t)si∗(t)]={0,i≠jPi,i=j其中，Pi为第i个信号源的功率。此时信号相关矩阵为：Rs=diag(P1,P2,⋯,PD)实际应用中，通常无法直接得到RX，只能得到样本的协方差矩阵RX‘=1N∑X(i)XH(i)Ni=1对于均匀线阵，阵元间距为d,假设波长为λ，传播速度为c,在二维平面内，入射方向与阵列法向的夹角为θ，延时可表示为：τm(t)=2π(m−1)dλsinθ对于入射角为θ的信号，该均匀线阵的导向矢量表示为：a(θ)=[1,e−jφ,⋯e−j(M−1)φ]T其中φ=2πdλsinθ当有L个信号源时，接收信号导向矢量矩阵A可以表示为：A=[1⋯1⋮⋱⋮e−j(M−1)φ1⋯e−j(M−1)φL]Music算法的基本思想则是将任意阵列输出数据的协方差进行特征分解，从而得到信号分量相对应的信号子空间和信号分量相正交的噪声子空间，然后利用这两个子空间的正交性来估计信号的参数。对阵列协方差矩阵进行特征分解：RX=ARSAH+σ2I与信号有关的特征值只有D个，分别等于矩阵ARSAH的特征值与σ2之和，其余的M−D个特征值为σ2，σ2是RX的最小特征值，即有D个特征值是与信号有关的，另外M−D个是与噪声有关的。因此，可以把RX的特征向量划分为信号特征向量与噪声特征向量。设λi=σ2是矩阵RX的第i个特征值并且是最小特征值，vi是λi相对应的特征向量，则有RXvi=λivi=σ2vi带入RX=ARSAH+σ2I得出σ2vi=(ARSAH+σ2I)vi展开后得到ARSAHvi=0将上式两边同乘以RS−1(AHA)−1AH后变成RS−1(AHA)−1AHARSAHvi=0于是又AHvi=0，i=D+1,D+2,⋯,M用个噪声特征向量为列，构造一个噪声矩阵ENEN=[vD+1,vD+2,⋯vM]定义空间谱Pmu(θ)Pmu(θ)=1aH(ϑ)ENENHa(θ)=1‖ENHa(θ)‖2当a(θ)和EN的各列正交是，该分母为零，但由于噪声的存在，它实际是为最小值，因此Pmu(θ)有一尖峰。变化θ，通过变化波峰来估计到达角的方法叫做MUSIC方法。二，仿真结果波峰为与角度为20度处。三，程序代码clcclearallformatlongN=200;doa=[2040]/180*pi;w=[pi/4pi/3]';M=8;P=length(w);lamda=150;d=lamda/2;snr=15;B=zeros(P,M);fork=1;PB(k,:)=exp(-j*2*pi*d*sin(doa(k))/lamda*[0:M-1]);endB=B';xx=2*exp(j*(w*[1:N]));x=B*xx;[pp,ppp]=size(x);xx=zeros(pp,ppp);cycle=200;forkkk=1:cyclexx=xx+awgn(x,snr);endx=xx/cycle;R=x*x';J=fliplr(eye(M));R=R+J*conj(R)*J;[U,V]=eig(R);UU=U(:,1:M-P);theta=-90:0.5:90;forii=1:length(theta)AA=zeros(1,length(M));forjj=0:M-1AA(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lamda);endWW=AA*UU*UU'*AA';Pmusic(ii)=abs(1/WW);endPmusic=10*log10(Pmusic/max(Pmusic)+eps);figure(1)plot(theta,Pmusic)xlabel('Angle\theta\degree')ylabel('P(\theta)/dB')title('Spectrum')gridon