阵列信号.docx

1、阵列信号2010 年 春 季学期研究生课程考核（读书报告、研究报告）考核科目: 阵列信号处理 报告题目: 阵列信号处理习题与实验 学生所在院（系）: 电子与信息工程学院 学生所在学科: 信息与通信工程 学 生 姓 名: 崔岩 学 号: 09S005068 学 生 类 别: 非委培 考核结果： 阅卷人： 实验一 自适应波束形成一 实验目的1、 掌握自适应波束形成的方法2、 了解影响自适应波束形成性能的因素二 实验内容假设阵元个数M=12，阵元间距为半波长，信号方向为0度，2个与信号不相关的干扰，方向分别为-30度和50度，仿真：（1）不同信噪比（SNR）下的窗函数；（2）不同干噪比（INR）下的

2、窗函数；（3）不同快拍数下的窗函数。三 实验原理3.1最小无失真方差响应最小无失真方差响应（MVDR）波束形成器阵列输出的最小均方值使波束受到限制，从而在给定信号的方向到来的波形产生无畸变的响应。采样矩阵求逆（SMI）方法。有：其中是数据的协方差矩阵的最大似然估计，s是目标的导向向量。其中，表示阵列模型在某一时刻的快拍矢量。3.2最小均方波束形成器最小均方准则条件下的权重每经过一个采样间隔都要增加一个增量。这个增量是与输出剩余功率的变化率成正比的。最小均方准则的自适应方法计算较为简单，并是慢慢收敛的。LMS算法如下图：图2.1 LMS算法方框图下面将要具体的介绍一下最小均方波束形成器的发展情况

3、。假设为一个N元阵接收到的的向量信号，波束形成器的输出是由一个复杂的权和每个阵元的输出相乘得到的，再将所有阵元相加，得到：通过权使阵列的输出功率最小，同时还要受到：的限制，其中C为的限制矩阵，f为的限制向量。为了在后面得到最小均方准则的结果，设：，需要解决的最优化问题就可以表示成：其约束条件为，这里表示协方差矩阵。最优解为：权向量w分为两部分：一个是受限子空间，另一个是与它正交的。即：维的矩阵B的列组成一个正交互补的矩阵C，（）。的向量可以改善维正交子空间的抗干扰性能。向量是固定值，并且。现在，上文提到的受限制的最优化问题可以表示为不受限制的问题：其最优解为：LMS算法步骤如下：1、设定初始静

4、态权值： 2、分别计算：、3、计算误差： 4、，其中为步长，为的协方差矩阵的最大特征值。5、重复步骤2到步骤5共N次，N为预设迭代次数。四 实验仿真及结果分析（1）不同信噪比（SNR）情况的窗函数由下述实验结果对比可知，在信干比一定，快拍数一定的条件下，当信号具有较好的信噪比条件时，三种自适应波束形成方法都能将波束指向信号方向，同时在干扰方向形成零陷。其中MVDR算法受信噪比影响很大，当信噪比很低时，旁瓣很高，对噪声的抑制效果不好；LMS算法形成的零陷最深，对干扰的抑制性能最好，并且旁瓣最接近最优波束形成方法，对噪声的抑制也很好。（a） （b）（c） （d）图4.1不同信噪比情况的窗函数SIR

5、=-10dB，快拍数L=1000（a）SNR=-20dB（b）SNR=-10dB（c）SNR=0dB（d）SNR=10dB（2） 不同干噪比（INR）情况的窗函数（a） （b）（c） （d）图4.2不同干噪比情况的窗函数SNR=0dB，快拍数L=1000（a）INR=6dB（b）INR=0dB（c）INR=-6dB（d）INR=-10dB由上述实验结果对比可知，在信噪比一定，快拍数一定的条件下，三种自适应波束形成方法都能将波束方向指向信号方向，同时在干扰方向形成零陷。其中MVDR算法受干燥比影响较大，旁瓣随着干燥比的增大而太高，降低了对噪声的音质性能；LMS算法受干燥比影响较小，旁瓣高度较为稳

6、定，对噪声和干扰的一直性能也较为稳定。（3） 不同快拍数情况的窗函数（a） （b）（c） （d）图4.3不同快拍数情况的窗函数SNR=0dB，SIR=-10dB（a）L=100（b）L=500（c）L=1000（d）L=2000由上述实验结果对比可知，在信噪比一定，信干比一定的条件下，三种自适应波束形成方法都能将波束方向指向信号方向，同时在干扰方向形成零陷。其中MVDR算法受快拍数影响较大，旁瓣随着快拍数的增大而越来越接近最优波束形成算法，这是因为MVDR需要用观测数据来估计干扰和噪声的协方差矩阵，快拍数越多，对协方差矩阵的估计就越准确；LMS算法受快拍数影响较小，这是因为LMS算法不需要估计

7、噪声和干扰的协方差矩阵。附：实验程序clear all;close all;clc% set parametersc=3e8;fc=3e6; % carrier frequencylambda=c/fc; % wavelengthfs=1e3; M=12; % number of elementsd=lambda/2; % interelement spacing% echostheta0=0; % signal directionbeta0=2*pi*d*sin(theta0*pi/180)/lambda;a0=exp(-1i*beta0*(0:M-1); % signal steering

8、 vectorL=100; %快拍次数fd=10; % signal frequencyt=0:L-1;s0=exp(1i*2*pi*fd*t/fs); % signaltheta1=-30; % interference directionbeta1=2*pi*d*sin(theta1*pi/180)/lambda;a1=exp(-1i*beta1*(0:M-1); % interference steering vectorfd1=15; % interference frequencys1=3.16*exp(1i*2*pi*fd1*t/fs); % interference 1theta

9、2=50; % interference directionbeta2=2*pi*d*sin(theta2*pi/180)/lambda;a2=exp(-1i*beta2*(0:M-1); % interference steering vectorfd2=20; % interference frequencys2=3.16*exp(1i*2*pi*fd2*t/fs); % interference 2 SIR=-10dB_3.16;SIR=-6dB_2;SIR=0dB_1;SIR=-6dB_0.56A=a0 a1 a2;S=s0;s1;s2;x=A*S;SNR=0; % SNR dBnoi

10、=zeros(1,L);x_in=a1 a2*s1;s2;noi=sqrt(0.5*10(SNR/10)*(randn(1,L);for m=1:M noi=sqrt(0.5*10(SNR/10)*(randn(1,L)+1i*randn(1,L); x(m,:)=x(m,:)+noi; x_in(m,:)=x_in(m,:)+noi;end% MVDRRx0=x*x/L;rxd=zeros(M,1);iRx0=pinv(Rx0);w_mvdr=iRx0*a0/(a0*iRx0*a0);% optR_in=x_in*x_in/L;iR_in=pinv(R_in);w_opt=(iR_in*a0

11、)/(a0*iR_in*a0);% lmsCon=a0 a1 a2; Ncon=size(Con,2); f=1;0;0; B=zeros(M,M-Ncon);D=Con;A=D(:,1:Ncon);for m=Ncon+1:M B(1:Ncon,m-Ncon)=-AD(:,m); B(m,m-Ncon)=1;endwq=Con/(Con*Con)*f;wa=zeros(M-Ncon,L);w_lms=zeros(M,L);y=zeros(1,L);yc=wq*x(:,1);z=B*x(:,1);alpha=1/M/sum(abs(x(:,1).2);wa(:,1)=alpha*z*conj(

12、yc);w_lms(:,1)=wq-B*wa(:,1);for k=2:L yc=wq*x(:,k); z=B*x(:,k); yp=yc-wa(:,k-1)*z; alpha=1/M/sum(abs(x(:,k).2); wa(:,k)=wa(:,k-1)+alpha*z*conj(yp); w_lms(:,k)=wq-B*wa(:,k);end% beamformingtheta=-90:1:90;N1=length(theta);y_mvdr=zeros(1,N1);y_opt=zeros(1,N1);y_lms=zeros(1,N1);for k=1:N1 beta=2*pi*d*si

13、n(theta(k)*pi/180)/lambda; a=exp(-1i*beta*(0:M-1); y_mvdr(k)=w_mvdr*a; y_opt(k)=w_opt*a; y_lms(k)=w_lms(:,L)*a;end% resultfigure,plot(theta,db(y_mvdr),-.)axis(-90 90 -100 0);hold onplot(theta,db(y_opt),r);hold onplot(theta,db(y_lms),g-);hold ongrid ontitle(方向图)xlabel(角度/)ylabel(归一化幅值/dB)legend(MVDR,

14、OPT,LMS);实验二 MISIC算法仿真一 实验目的1、 掌握MUSIC算法进行DOA估计的基本原理2、 了解影响DOA估计性能的因素二 实验内容假设阵元个数M=12，信号方向分别为0度、2度和-35度，噪声为高斯噪声仿真：（1）快拍次数一定，分辨力与SNR的关系（2）快拍次数一定，分辨力与孔径的关系；（3）SNR一定，分辨力与快拍的关系。三 实验原理窄带远场信号的DOA数学模型为：阵列数据的协方差矩阵为： 由于信号与噪声相互独立，数据协方差矩阵可以分解为与信号、噪声相关的两部分，其中是信号的协方差矩阵，是信号部分。对进行特征分解有：式中，是由大特征值对应的特征矢量张成的子空间也即信号子空

15、间，而是由小特征值对应的特征矢量张成的子空间也即噪声子空间。根据前面所述的性质2可知，在理想的条件下数据空间中的信号子空间与噪声子空间是相互正交的，即信号空间中的导向矢量也与噪声子空间正交经典的MUSIC算法正是基于上述这个性质提出的，但考虑到实际接收数据矩阵是有限快拍数的，即数据协方差矩阵的最大似然估计为 对进行特征分解可以计算得到噪声子空间特征矢量矩阵。由于噪声的存在，与并不能完全的正交。因此，实际上求DOA是以最小优化搜索实现的，即 所以，MUSIC算法的谱估计公式为：MUSIC算法的计算步骤如下：（1） 由阵列的接收数据得到数据协方差矩阵；（2） 对进行特征分解；（3） 由的特征值进行

16、信号源数的判断；（4） 确定信号子空间与噪声子空间；（5） 根据信号参数的范围进行谱峰搜索；（6） 找出极大值点对应的角度就是信号入射方向。对于上述的MUSIC算法，还应该要注意以下几点：对于非理想情况下得到协方差矩阵的特征值满足下式：所以判断信号源数需要用到有关信号源估计的方法来进行信号源数的确定。线阵的信号参数搜索范围为，而面阵的搜索范围为。另外，一个确定阵列的导向矢量由阵元的位置唯一的确定；MUSIC算法的一种归一化形式，即：在实际应用中，对于一维导向矢量有下式成立： 四 实验仿真及结果分析（1） 快拍次数一定，分辨力与SNR的关系快拍次数L=1000，阵元间距为半波长，信噪比SNR在-

17、10dB,30dB之间变化，得到的仿真结果如图1所示。图1 分辨率与SNR的关系总图（左）、分辨率与SNR的关系切面图（右）由图1可知，当信噪比较差，小于10dB时，主瓣内的两个信号（0度和2度）不能区分开，被误判为一个信号；当信噪比大于10dB时，可以估计出三个信号的角度，实现超分辨。（2） 快拍次数一定，分辨力与孔径的关系快拍次数L=1000，信噪比为20dB，阵元间距与半波的比在0.1，1之间变化，得到的仿真结果如图2所示。 图2 分辨率与孔径的关系总图（左）、分辨率与孔径的关系切面图（右）由图2可知，当阵元间距较小，小于0.2倍的波长时，0度和2度的信号被误判为一个信号，这是因为阵元间

18、距太小，阵列的孔径就小，主瓣很宽，角度分辨率变差；随着阵元间距的增大，角度分辨率提高，当阵元间距大于0.3倍的波长，小于0.6倍的波长时，能很好的估计出三个信号的方向；当阵元间距大于0.6倍的波长时，出现了虚警信号，这是由于阵元间距太大，产生栅瓣。（3） SNR一定，分辨力与快拍数的关系信噪比为20dB，阵元间距为半波长，快拍数在1,400之间变化，得到的仿真结果如图3所示。图3 分辨率与快拍的关系总图（左）、分辨率与快拍的关系切面图（右）由图3可知，当快拍数太少，少于50个时，会出现信号遗漏的现象，只有当快拍数足够多时，才能对三个信号做出正确的角度估计，快拍次数越多，角度分辨率越好。附：实验

19、程序实验一程序clear all;close all;clctheta=-90:0.1:90;Ltheta=length(theta);SNR=-10:30;ps=2;Lsnr=length(SNR);temp=10.(SNR/10);sigmaN=sqrt(ps./temp);DOA=zeros(Ltheta,Lsnr);N_x=1000; N=12; l=1.8;c=3e8;w0=2*pi*c/l;d=0.5*l;M=3; % Number of Signals w=pi/6 pi/10 pi/8;xx=pi/180*2;yy=0;zz=-pi/180*35; B=1 exp(-j*2*p

20、i*d*sin(xx)/l) exp(-j*2*2*pi*d*sin(xx)/l) exp(-j*2*3*pi*d*sin(xx)/l). exp(-j*2*4*pi*d*sin(xx)/l) exp(-j*2*5*pi*d*sin(xx)/l) exp(-j*2*6*pi*d*sin(xx)/l). exp(-j*2*7*pi*d*sin(xx)/l); 1 exp(-j*2*pi*d*sin(yy)/l) exp(-j*2*2*pi*d*sin(yy)/l) exp(-j*2*3*pi*d*sin(yy)/l). exp(-j*2*4*pi*d*sin(yy)/l) exp(-j*2*5*

21、pi*d*sin(yy)/l) exp(-j*2*6*pi*d*sin(yy)/l). exp(-j*2*7*pi*d*sin(yy)/l); 1 exp(-j*2*pi*d*sin(zz)/l) exp(-j*2*2*pi*d*sin(zz)/l) exp(-j*2*3*pi*d*sin(zz)/l). exp(-j*2*4*pi*d*sin(zz)/l) exp(-j*2*5*pi*d*sin(zz)/l) exp(-j*2*6*pi*d*sin(zz)/l). exp(-j*2*7*pi*d*sin(zz)/l); xxx=2*exp(j*(w+w0)*0:N_x-1); for ii=

22、1:Lsnrx=B*xxx+sigmaN(ii)*randn(8,N_x)+j*randn(8,N_x)*sigmaN(ii); R=x*x/N_x; V D=eig(R);lambda,index = sort(diag(D); UU=V(:,index(1:5); for i = 1:length(theta)AA=1 exp(-j*2*pi*d*sin(theta(i)/180*pi)/l) exp(-j*2*2*pi*d*sin(theta(i)/180*pi)/l). exp(-j*2*3*pi*d*sin(theta(i)/180*pi)/l) exp(-j*2*4*pi*d*si

23、n(theta(i)/180*pi)/l). exp(-j*2*5*pi*d*sin(theta(i)/180*pi)/l) exp(-j*2*6*pi*d*sin(theta(i)/180*pi)/l) . exp(-j*2*7*pi*d*sin(theta(i)/180*pi)/l);WW=AA*UU*UU*AA;Pmusic(i)=abs(1/WW);endPmusic = 10*log10(Pmusic);DOA(:,ii)=Pmusic;endaa=theta*ones(1,Lsnr);figure (1)subplot(2,2,1)plot(aa,DOA(:,10)xlabel(t

24、heta/) ylabel(SNR=0dB)subplot(2,2,2)plot(aa,DOA(:,20)xlabel(theta/) ylabel(SNR=10dB)subplot(2,2,3)plot(aa,DOA(:,30)xlabel(theta/) ylabel(SNR=20dB) subplot(2,2,4)plot(aa,DOA(:,40)xlabel(theta/) ylabel(SNR=30dB)aa=theta*ones(1,Lsnr);bb=ones(Ltheta,1)*SNR;figure(2)mesh(aa,bb,DOA);title(MUSIC-3D图(theta,

25、SNR);xlabel(theta/) ylabel(SNR/dB)zlabel(PMUSIC)实验二程序clear all; close all;clc;bili=0.1:0.01:1; Lbili=length(bili);theta=-90:0.1:90;Ltheta=length(theta);SNR=20ps=2;temp=10.(SNR/10);sigmaN=sqrt(ps./temp);DOA=zeros(Ltheta,Lbili);N_x=1000; N=12; l=1.8;c=3e8;w0=2*pi*c/l;d=bili*l;M=3; w=pi/6 pi/10 pi/8;xx

26、=pi/180*2;yy=0;zz=-pi/180*35; for ii=1:Lbili sss=bili(ii)*2*pi;B=1 exp(-j*sss*sin(xx) exp(-j*2*sss*sin(xx) exp(-j*3*sss*sin(xx). exp(-j*4*sss*sin(xx) exp(-j*5*sss*sin(xx) exp(-j*6*sss*sin(xx) exp(-j*7*sss*sin(xx); 1 exp(-j*sss*sin(yy) exp(-j*2*sss*sin(yy) exp(-j*3*sss*sin(yy). exp(-j*4*sss*sin(yy) e

27、xp(-j*5*sss*sin(yy) exp(-j*6*sss*sin(yy) exp(-j*7*sss*sin(yy); 1 exp(-j*sss*sin(zz) exp(-j*2*sss*sin(zz) exp(-j*3*sss*sin(zz). exp(-j*4*sss*sin(zz) exp(-j*5*sss*sin(zz) exp(-j*6*sss*sin(zz) exp(-j*7*sss*sin(zz); xxx=2*exp(j*(w+w0)*0:N_x-1);x=B*xxx+sigmaN*randn(8,N_x)+j*randn(8,N_x)*sigmaN; R=x*x/N_x

28、; V D=eig(R);lambda,index = sort(diag(D);UU=V(:,index(1:5); for i = 1:length(theta)AA=1 exp(-j*sss*sin(theta(i)/180*pi) exp(-j*2*sss*sin(theta(i)/180*pi). exp(-j*3*sss*sin(theta(i)/180*pi) exp(-j*4*sss*sin(theta(i)/180*pi). exp(-j*5*sss*sin(theta(i)/180*pi) exp(-j*6*sss*sin(theta(i)/180*pi). exp(-j*

29、7*sss*sin(theta(i)/180*pi);WW=AA*UU*UU*AA;Pmusic(i)=abs(1/WW);endPmusic = 10*log10(Pmusic);DOA(:,ii)=Pmusic;endaa=theta*ones(1,Lbili);bb=ones(Ltheta,1)*bili;subplot(2,2,1)plot(aa,DOA(:,15)xlabel(theta/) ylabel(d/lambda=0.15)subplot(2,2,2)plot(aa,DOA(:,40)xlabel(theta/) ylabel(d/lambda=0.4)subplot(2,2,3)plot(aa,DOA(:,50)xlabel(theta/) ylabel(d/lambda=0.5)subplot(2,2,4)plot(aa,DOA(:,70)xlabel(theta/) ylabel(d/lambda=0.7)figure(2)mesh(aa,bb,DOA)