阵列式信号处理 MUSIC算法仿真VIP

阵列式信号处理 MUSIC算法仿真渤海船舶职业学院 常 亮摘 要：波达方向(DOA)估计是阵列信号处理研究的重要内容之一多重信号分类（MUSIC）算法是一种经典有效的DOA方法。概述了阵列信号处理研究的主要内容，详细介绍了多重信号分类算法的原理，并用MATLAB进行仿真实现该算法，研究了MUSIC算法的性质。关键词：阵列信号处理 多重信号分类算法 仿真1、实验目的（1）掌握MUSIC算法进行DOA估计的基本原理（2）了解影响DOA估计性能的因素2、实验内容假设阵元个数M=16，信号方向分别为0度、3度和-35度，噪声为高斯噪声仿真：（1）快拍次数一定，分辨力与SNR的关系（2）快拍次数一定，分辨力与孔径的关系；（3）SNR一定，分辨力与快拍的关系。3、实验原理窄带远场信号的DOA数学模型为 阵列数据的协方差矩阵为 由于信号与噪声相互独立，数据协方差矩阵可以分解为与信号、噪声相关的两部分，其中是信号的协方差矩阵，是信号部分。对进行特征分解有 式中，是由大特征值对应的特征矢量张成的子空间也即信号子空间，而是由小特征值对应的特征矢量张成的子空间也即噪声子空间。根据前面所述的性质2可知，在理想的条件下数据空间中的信号子空间与噪声子空间是相互正交的，即信号空间中的导向矢量也与噪声子空间正交经典的MUSIC算法正是基于上述这个性质提出的，但考虑到实际接收数据矩阵是有限快拍数的，即数据协方差矩阵的最大似然估计为 对进行特征分解可以计算得到噪声子空间特征矢量矩阵。由于噪声的存在，与并不能完全的正交。因此，实际上求DOA是以最小优化搜索实现的，即 所以，MUSIC算法的谱估计公式为 下面给出MUSIC算法的计算步骤：（1） 由阵列的接收数据得到数据协方差矩阵；（2） 对进行特征分解；（3） 由的特征值进行信号源数的判断；（4） 确定信号子空间与噪声子空间；（5） 根据信号参数的范围进行谱峰搜索；（6） 找出极大值点对应的角度就是信号入射方向。对于上述的MUSIC算法，还应该要注意以下几点：对于非理想情况下得到协方差矩阵的特征值满足下式：所以判断信号源数需要用到有关信号源估计的方法来进行信号源数的确定。线阵的信号参数搜索范围为，而面阵的搜索范围为。另外，一个确定阵列的导向矢量由阵元的位置唯一的确定；MUSIC算法的一种归一化形式，即 在实际应用中，对于一维导向矢量有下式成立： 4、实验仿真及结果分析4.1 MUSIC随信噪比变化情况信号数共3个，分别位于2，30和31，它们功率相同，阵元数共16个。这里快拍数取500，以保证快拍数足够，可以较好得得出统计特性，以去除采样不足带来的分辨力低下。信噪比分别取为SNR=0,5,10,15,20dB，得到图2。从图中可以看到，随信噪比增大，MUSIC谱图旁瓣变低，主瓣变得更加尖锐，分辨力变得更好。 信噪比为0dB 信噪比为30dB 信噪比为50dB 信噪比为-20dB信噪比为0dB时经典music算法解出的来波角度为2.0700，31.7000，求根music算法解出的来波角度为30.5268，1.9996信噪比为30dB时经典music算法解出的来波角度为2.0700，30.9900，32.3200，求根music算法解出的来波角度为29.8490，1.9998，31.1206。信噪比为50dB时经典music算法解出的来波角度为2.0700，31，32.3100，求根music算法解出的来波角度为29.8602，2.0001，31.1104。信噪比为-20dB时经典music算法和求根music算法都无法解出来波角度。由上面的图形可以看出，信噪比越大，越有利于信号角度的分辨，以及分辨出角度差很小的两个信号，表现为峰值十分尖锐，而随着信噪比的降低，分辨力下降，当为-20dB时已经无法分辨出30和31的两个信号，更糟糕的是还出现了很多假的峰值，导致对信号源数的估计错误。4.2 MUSIC谱图随快拍数变化情况在上面的实验中，快拍数为5000，分辨效果较好，原因是当快拍足够时，能较好地得到噪声的统计特性，从而精确地分解出信号子空间和噪声子空间，从而较好的分辨信源方向，当快拍数下降时，分辨效果将会下降。 快拍数为300 快拍数为100可见当快拍数下降的时候，峰值将会没有那么尖锐，对信号的分辨力下降，甚至无法两个角度差较小的信号，原因是小样本时间平均与统计平均差距较大，噪声子空间分解不准确。5、实验结论由上面的实验可以看出，music算法可以实现对于DOA的求解，信噪比越大，快拍数越高越有利于分辨力的提高。此外不同的music算法求解时，计算时间和精度都有差别，但是由于来波来自3个方向，每个方向的求解实际上是存在相互牵制的关系，比如来自0的信号求解很准可能就会造成-20的信号求解不准。而信号的窄带信也是很关键的，当不满足窄带条件时，求解的角度会误差很大。6、实验程序clear all;clc;close all;M=16;K=3;%NN=M-K;%m=0:1:M-1;wavelength=1;%f=300MHzf=3*108/wavelength;fs=10*f;Ts=1/fs;df=f/30;%df=10MHzd=wavelength/2;%k=2*pi/wavelength;%N=500;%1%SNR=-20;%dBSCR2=0;%21/dBSCR3=0;%31/dBnoise=10(-SNR/10);%pj2=10(SCR2/10);%2pj3=10(SCR3/10);%3as2=sqrt(pj2);%2as3=sqrt(pj3);%3sigma=sqrt(noise);%fai_j1=2/180*pi;%1fai_j2=30/180*pi;%2fai_j3=31/180*pi;%3w1=f*2*pi;w2=(f+df)*2*pi;w3=(f-df)*2*pi;kj1=k*sin(fai_j1);kj2=k*w2/w1*sin(fai_j2);kj3=k*w3/w1*sin(fai_j3);akj1=exp(-1j*(kj1*m*d).;akj2=exp(-1j*(kj2*m*d).;akj3=exp(-1j*(kj3*m*d).;R=zeros(M,M);Ri=zeros(M,M);for i=1:N ri=normrnd(0,sigma,M,1); Ri=(akj1*exp(1j*w1*i*Ts)+as2*akj2*exp(1j*w2*i*Ts)+as3*akj3*exp(1j*w3*i*Ts)+ri)*(akj1*exp(1j*w1*i*Ts)+as2*akj2*exp(1j*w2*i*Ts)+as3*akj3*exp(1j*w3*i*Ts)+ri); R=R+Ri;endR=R/N;V,D=eigs(R,NN,sm);fai=(-90:0.005:90)/180*pi;ps=zeros(1,length(fai);disp(music/);for i=1:length(fai)% ki=k*w3/w1*sin(fai(i); aki=exp(-1j*(ki*m*d).; ps(i)=real(1/(aki*(V*V)*aki);endfor i=2:length(fai)-1% if ps(i)ps(i+1)&ps(i)ps(i-1)&ps(i)max(ps)/20 disp(fai(i)/pi*180);% endend figure;plot(fai/pi*180,10*log10(ps/max(ps);%syms z;Pz1=sym(zeros(1,M);Pz2=sym(zeros(1,M).; for i=1:M Pz1(i)=z(1-i); Pz2(i)=z(i-1); end tsym=solve(Pz1*(V*V)*Pz2); tdouble=double(tsym); r=abs(tdouble); rnum=length(r);%numsolve=zeros(1,K);% for i=1:K% temp1=20;% temp2=1;% for j=1:length(r)% if abs(r(j)-1)temp1% temp1=abs(r(j)-1);% temp2=j;% end% end% numsolve(i)=tdouble(temp2);% r(temp2)=10;% end tsolve=zeros(1,rnum);j=0;for i=1:rnum if(abs(r(i)-1)0.1) j=j+1; tsolve(j)=tdouble(i); endendnumsolve=tsolve(1:j);signal_position=asin(-wavelength/(2*pi*d)*angle(numsolve)/pi*180;real_position=signal_position(1);for i=2:length(signal_position)% flag=0; for j=1:i-1 if abs(signal_position