自适应波束形成仿真

1、课程设计（论文）任务书专 业电子信息工程班级0802503学生080250303指导教师刘帅、王军、周洪娟题 目自适应波束形成算法仿真子 题设计时间2011年11月28日 至 2011年12月9日 共2 周设计要求设计（论文）的任务和基本要求，包括设计任务、查阅文献、方案设 计、说明书（计算、图纸、撰写内容及规范等）、工作量等内容。设计任务：实现自适应波束形成算法，并考察算法性能。参考文献：现代信号处理，雷达原理，Matlab相关书籍等方案设计：1）由均匀线阵形式，确定阵列的导向矢量；2）由阵列导向矢量，对接收信号进行建模仿真；3）根据约束条件，实现自适应波束形成算法；4）考察算法性能与信噪比

2、，采样率，观测时间等参数的关系。说明书：设计报告内容要求包括：信号模型的建立、分析过程；建模源程序；仿真结果分析。本专业综合设计要求在两周内完成。指导教师签字:系（教研室）主任签字:年 月曰刖言近十年，自适应波束形成算法在通信雷达、声纳、生物医学工程等科技领域 中到了极为广泛的应用。在实际应用中，如果信号源、天线阵列出现误差，传统 的自适应波束形成算法性能将会下降。但是对于稳健的自适应波束形成算法。环 境及天线阵列的误差和不确定性是必须要考虑的关键问题。这里对稳健的自适应 波束形成算法的研究现状与发展动态进行了较为详细的评述。传统的算法在设计波束形成器时，都是假定在其训练数据中不含有期望信 号

3、。这时的波束形成器对于阵列响应误差和有限次快拍数据的稳健性非常好。但 是在许多情况下，干扰和噪声信号的观测数据不可避免地要被期望信号所污染。 而即使在理想情况下，这时精确的知道期望信号的方向向量，在有限次快拍训练 数据中含有的响应之间的不匹配。自适应阵列信号处理对于这类误差的反应非 常敏感。因为这时期望信号分量会被当作干扰抑制掉。阵列校正误差及观测方向 的偏差等现象在实际中经常出现。引起自适应波束形成算法性能降低的原因除了阵列响应误差外，由于信道的 不稳定、干扰和天线的移动、天线的转动等所造成的波束形成器训练数据的不稳 定性也是一个主要原因。首先，即使在训练数据 中不包含期望信号，这种不稳 定

4、性也会限制自适应波束形成算法的性能；进而，当在训练数据中含有期望信号 时，波束形成器性能会进一步严重降低；最后，如果出现快速移动的干扰。波束 形成器的性能将会崩亏。因为这时自适应权向量的收敛速度跟不上干扰的变化， 不能有效的对干扰形成零陷。 TOC o 1-5 h z 前言2 HYPERLINK l bookmark17 o Current Document 设计原理4设计过程7使用MVDR算法7EC 算法9GUI界面设计12 HYPERLINK l bookmark41 o Current Document 设计总结13 HYPERLINK l bookmark44 o Current Do

5、cument 参考文献14一.设计原理传统算法，波束形成器的输出可表示为y(k)= wh，Q)( i)，这里复向量 x( k)是在观测阵元上第k次快拍的采样值，它是Mx l维的复向量；W是波束形 成器Mx l维的权向量设阵元数目为M，( - ) n代表共轭转置，其中向量x ( k ) 可表示为(k ) = s ( k ) + i ( k ) + n ( k )。这里S ( k )，i ( k )，n ( k)分别代表期 望信号、干扰和噪声，它们都是统计独立的分量。当电波从点辐射源以球面波向 外传播，只要离辐射源足够远，在接收的局部区域，球面波就可以近似为平面波， 通信信号的传播一般都满足这一远

6、场条件。这时有 七(k)=s(k方其中七是信号 的方向向量。则最佳权向量可以通过使信十噪比(S I NR )的最大化来获得：SINRwH R wwH R w(2)g (k)S h (k)，R. AE t (k )+ nk )L(k )+ nk )H )，它们分别是 MxM维的信号协方差矩阵、干扰加噪声协方差矩阵，其中E.代表数学期望。为 了求得式(2 )中S I N R的最大值，可以通过在约束条件wHRw = 1的约束下，求 满足mi nMR+nw(3 )的权向量，也就是使输出信噪比的功率最小。当rank *= 1时，式(3 )的约束条件可简化为wHas = 1这种方法称为最小方差无畸变响应

7、(MVDR : minimum variance distortionless response )波束形成。式(3 )的结果为W”2(4 )式中，a =：*；=)在实际计算中，式(4)中的协方差矩阵Ri+n用有限次 快拍数据得到的估计值R代替RA=曷 x(n)xh (n)= XX hN n =1N(5)式中X=【x ( 1 )，x ( 2 )，x ( N )】是Mx N维矩阵，N为快拍数。稳健的自适应波束形成算法，目前最流行的算法，是使用对角加载(diagonal loading )技术，它在小快拍数据和存在阵列信号响应误差的情况下。仍然具有 较好的稳健性。它的核心思想是：通过在目标函数式(

8、3 )中增加一个二次补偿wHR w = 1项来调整权向量。可以表示为Am i nw H R w + y w H w(6 )式中y为加载因子。式(6 )的解为(A )T w = R+yI a DL ) s(7 )式中I为单位矩阵。对角加载技术能够提高自适应波束形成算法的性能，但是这种方法的 缺陷是对于参数y的选择不可靠。而且也不易获得。另外一个比较流行的稳健自适应波束形成算法是基于特征空间的波束形成算法。这种方法仅适用于点信号源。其核心思想是利用方向向量在采样信号与干扰信号子空 间上的投影来代替信号的方向向量。式(5)的特征分解为AAAA H AAA HR = EA E +GT G式中MX (

9、L + 1 )维矩阵E包含R中的L + 1个信号与干扰特征向量，A为其特征值构成的(L + I ) x ( L + 1 )维对角矩阵；MX ( M L 1 )维矩阵G包含R中的ML一1个噪声特征向量，T为其特征值构成的对角矩阵。假设干扰信号数为L ,则基于特征空间的波束形成算法的权向量可 表示为-1A HA A H式V= A空间上的正交投影。如果环境条件匹配E = E E基于特征空间的波束形成算法的性能会大幅度的中，是在信号与干扰子提高。但是这种方法的致命缺陷是要求数据必须是基于低秩的固定模型且十扰源数必须精确知道，而且这种方法只适用于高信噪比的情况下，当信噪 比较低时，由于可能会发生子空间交

10、换而导致投影矩阵失效。针对上述算法中的缺陷，已经提出了许多改进的算法。为了解决权向量的不稳定性，提出了一种改进算法 但是这种算法仅适用于慢速移动的十 扰；为了对付快速移动的干扰，要用人为的方法加宽零陷波瓣的宽度的方 法。以上的稳健的自适应波束形成算法都是一些特殊的算法，只适用于特定 的环境。只是在最近才出现一些严格意义上的稳健的波束形成算法，这种算 法的基本思想是：通过定义所谓的不确定集和最优化最差情况下的性能。这里我们假设信号和干扰协方差矩阵、噪声协方差矩阵都存在有误差， 即这些矩阵的实际值与估计值之间存在一定的失配，即：R = R + A ，R = R + A式中Rs和R.+.分别是信号和

11、干扰协方差矩阵、噪声协方差矩阵的估计值，和R1 n为 Hermitian 矩阵。令|气|气|，阵的不确定集。为了获得波束形成器的最佳权向量是它们的实际值。气和 气是未知的矩阵误差，它们均它们分别是信号和干扰协方差矩阵、噪声协方差矩 可以通过使最差情况下的输出信干噪比（S I N R ）为最大。即m n w H（R + A ）wm a x A1 7T-m a xwH（R. + a ）wv|A II1 ,|A2假设信号方向向量存在误差6 = a - a ,监 | 1, V6类似与式（8可得到稳健的MVDR波束形成问题的解为bXaHs人A -R + 人 21I 7、-1 aA -R +人 21I 7

12、-1上式表明，最差情况性能的最优化方法实质上是种新的对角加载方法。 该方法有考虑Rs的非负性，但是它并不影响该波束形成器在大多数环境下 的性能。设计过程(1)使用MVDR算法%MVDR算法M=16;fc=3e9;d=0.05;lamda=0.1;SNR = 0;JNR1 = 30;JNR2 = 30;JNR3 = 30;ang0 = 15;angl = 25;ang2 = 0;ang3 = -15;K=100;fs=18e9;% a0期望信号方向矢量% a1干扰信号1方向矢量% a2干扰信号2方向矢量% a3干扰信号3方向矢量a0=zeros(M,1);a1=zeros(M,1);a2=zer

13、os(M,1);a3=zeros(M,1);j=sqrt(-1);for i=1:Ma0(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang0*pi/180);a1(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang1*pi/180);a2(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang2*pi/180);a3(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang3*pi/180);endA=a0,a1,a2,a3;A1=a1,a2,a3;t=(1:1:K)/fs;s0=cos(2*pi*fc*t);s01=ra

14、ndn(1,K).*sin(2*pi*fc*t+12);s02=randn(1,K).*sin(2*pi*fc*t+12);s03=randn(1,K).*sin(2*pi*fc*t+12);Ps0=s0*s0/K;%计算期望信号的功率%计算干扰信号1的功率%计算干扰信号2的功率%计算干扰信号3的功率%具有输入信噪比的期望信号 %具有十噪比1的干扰信号1 %具有十噪比2的干扰信号2 %具有十噪比3的干扰信号3Ps1=s01*s01/K;Ps2=s02*s02/K;Ps3=s03*s03/K;Ps=10.A(SNR/10);Pj1=10.人(JNR1/10);Pj2=10.A(JNR2/10);

15、s=s0/Ps0A0.5*PsA0.5;Pj3=10.A(JNR3/10);s1=s01/Ps1A0.5*Pj1A0.5;s2=s02/Ps2A0.5*Pj2A0.5;s3=s03/Ps3A0.5*Pj3A0.5;S=s;s1;s2;s3;S1=s1;s2;s3;numJam=size(S1,1);Noise=randn(M,K);X=A*S+Noise;X1=A1*S1+Noise;R=X*X/K;Rjn=X1*X1/K;W1=inv(Rjn)*a0/(a0*inv(Rjn)*a0);%mvdr 算法ang=-90:0.5:90;a=zeros(M,1);for k=1:length(ang

16、)for i=1:Ma(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang(k)*pi/180);endbeam1(k)=W1*a;endbeam1=abs(beam1);beam1=beam1/max(beam1);beam1=20*log10(beam1);plot(ang,beam1)(2)M=16;fc=3e9;d=0.05;lamda=0.1;SNRJNR1JNR2JNR3=0;=30;=30;=30;ang0 angl ang2 ang3=15;=25;=0;=-15;K=100;fs=18e9;a0=zeros(M,1);a1=zeros(M,1);a2=ze

17、ros(M,1);a3=zeros(M,1);%a0期望信号方向矢量% a1干扰信号1的方向矢量% a2干扰信号2的方向矢量% a3干扰信号3的方向矢量j=sqrt(-1);for i=1:Ma0(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang0*pi/180);a1(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang1*pi/180);a2(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang2*pi/180);a3(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang3*pi/180);endA=a0,a1,a

18、2,a3;A1=a1,a2,a3;t=(1:1:K)/fs;s0=cos(2*pi*fc*t);s01=randn(1,K).*sin(2*pi*fc*t+12);s02=randn(1,K).*sin(2*pi*fc*t+12);s03=randn(1,K).*sin(2*pi*fc*t+12);%计算期望信号的功率%计算干扰信号1的功率%计算干扰信号2的功率%计素干扰信号3的功率%具有输入信噪比的期望信号%具有十噪比1的干扰信号1%具有十噪比2的干扰信号2%具有十噪比3的干扰信号3Ps0=s0*s0/K;Ps1=s01*s01/K;Ps2=s02*s02/K;Ps3=s03*s03/K;P

19、s=10.人(SNR/10);Pj1=10.人(JNR1/10);Pj2=10.人(JNR2/10);Pj3=10.人(JNR3/10);s=s0/Ps0人0.5*Ps人0.5;s1=s01/Ps1人0.5*Pj1人0.5;s2=s02/Ps2人0.5*Pj2人0.5;s3=s03/Ps3人0.5*Pj3人0.5;S=s;s1;s2;s3;S1=s1;s2;s3;numJam=size(S1,1);Noise=randn(M,K);X=A*S+Noise;X1=A1*S1+Noise;Rjn=X1*X1/K;V,D=eig(Rjn);r=M-1;%对特征值从大到小排列并对对应的特征向量进行调整

20、for rr=1:Mfor i=1:M-1if D(i,i)D(i+1,i+1)Dtemp=D(i+1,i+1);D(i+1,i+1)=D(i,i);D(i,i)=Dtemp;Vtemp=V(:,i+1);V(:,i+1)=V(:,i);V(:,i)=Vtemp;endendendUs=V(:,(M-numJam+1:M);Ds=D(M-numJam+1:M),(M-numJam+1:M);T=Us,a0;Un=V(:,(1:M-numJam);Dn=D(1:M-numJam),(1:M-numJam);dd=diag(D);dd,IX=sort(dd);dd=flipud(dd);IX=flipud(IX);U=;for i=1:rU=U,V(:,IX(i);endW1=(eye(M)-U*U)*a0;% EC 算法ang=-90:0.5:90;a=zeros(M,1);for k=