阵列信号.docx
《阵列信号.docx》由会员分享，可在线阅读，更多相关《阵列信号.docx（22页珍藏版）》请在冰豆网上搜索。
阵列信号
2010年春季学期研究生课程考核
（读书报告、研究报告）
考核科目
:
阵列信号处理
报告题目
:
阵列信号处理习题与实验
学生所在院（系）
:
电子与信息工程学院
学生所在学科
:
信息与通信工程
学生姓名
:
崔岩
学号
:
09S005068
学生类别
:
非委培
考核结果
：
阅卷人：
实验一自适应波束形成
一实验目的
1、掌握自适应波束形成的方法
2、了解影响自适应波束形成性能的因素
二实验内容
假设阵元个数M=12，阵元间距为半波长，信号方向为0度，2个与信号不相关的干扰，方向分别为-30度和50度，仿真：
（1）不同信噪比（SNR）下的窗函数；
（2）不同干噪比（INR）下的窗函数；（3）不同快拍数下的窗函数。
三实验原理
3.1最小无失真方差响应
最小无失真方差响应（MVDR）波束形成器阵列输出的最小均方值使波束受到限制，从而在给定信号的方向到来的波形产生无畸变的响应。
采样矩阵求逆（SMI）方法。
有：
其中是数据的协方差矩阵的最大似然估计，s是目标的导向向量。
其中，表示阵列模型在某一时刻的快拍矢量。
3.2最小均方波束形成器
最小均方准则条件下的权重每经过一个采样间隔都要增加一个增量。
这个增量是与输出剩余功率的变化率成正比的。
最小均方准则的自适应方法计算较为简单，并是慢慢收敛的。
LMS算法如下图：
图2.1LMS算法方框图
下面将要具体的介绍一下最小均方波束形成器的发展情况。
假设为一个N元阵接收到的的向量信号，波束形成器的输出是由一个复杂的权和每个阵元的输出相乘得到的，再将所有阵元相加，得到：
通过权使阵列的输出功率最小，同时还要受到：
的限制，其中C为的限制矩阵，f为的限制向量。
为了在后面得到最小均方准则的结果，设：
，需要解决的最优化问题就可以表示成：
其约束条件为，这里表示协方差矩阵。
最优解为：
权向量w分为两部分：
一个是受限子空间，另一个是与它正交的。
即：
维的矩阵B的列组成一个正交互补的矩阵C，（）。
的向量可以改善维正交子空间的抗干扰性能。
向量是固定值，并且。
现在，上文提到的受限制的最优化问题可以表示为不受限制的问题：
其最优解为：
LMS算法步骤如下：
1、设定初始静态权值：
2、分别计算：
、
3、计算误差：
4、
，
，其中
为步长，
，
为
的协方差矩阵
的最大特征值。
5、重复步骤2到步骤5共N次，N为预设迭代次数。
四实验仿真及结果分析
（1）不同信噪比（SNR）情况的窗函数
由下述实验结果对比可知，在信干比一定，快拍数一定的条件下，当信号具有较好的信噪比条件时，三种自适应波束形成方法都能将波束指向信号方向，同时在干扰方向形成零陷。
其中MVDR算法受信噪比影响很大，当信噪比很低时，旁瓣很高，对噪声的抑制效果不好；LMS算法形成的零陷最深，对干扰的抑制性能最好，并且旁瓣最接近最优波束形成方法，对噪声的抑制也很好。
（a）（b）
（c）（d）
图4.1不同信噪比情况的窗函数SIR=-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算法受干燥比影响较小，旁瓣高度较为稳定，对噪声和干扰的一直性能也较为稳定。
（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算法不需要估计噪声和干扰的协方差矩阵。
附：
实验程序
clearall;
closeall;
clc
%%setparameters
c=3e8;
fc=3e6;%carrierfrequency
lambda=c/fc;%wavelength
fs=1e3;
M=12;%numberofelements
d=lambda/2;%interelementspacing
%%echos
theta0=0;%signaldirection
beta0=2*pi*d*sin（theta0*pi/180）/lambda;
a0=exp（-1i*beta0*（0:
M-1）'）;%signalsteeringvector
L=100;%快拍次数
fd=10;%signalfrequency
t=0:
L-1;
s0=exp（1i*2*pi*fd*t/fs）;%signal
theta1=-30;%interferencedirection
beta1=2*pi*d*sin（theta1*pi/180）/lambda;
a1=exp（-1i*beta1*（0:
M-1）'）;%interferencesteeringvector
fd1=15;%interferencefrequency
s1=3.16*exp（1i*2*pi*fd1*t/fs）;%interference1
theta2=50;%interferencedirection
beta2=2*pi*d*sin（theta2*pi/180）/lambda;
a2=exp（-1i*beta2*（0:
M-1）'）;%interferencesteeringvector
fd2=20;%interferencefrequency
s2=3.16*exp（1i*2*pi*fd2*t/fs）;%interference2SIR=-10dB_3.16;SIR=-6dB_2;SIR=0dB_1;SIR=-6dB_0.56
A=[a0a1a2];
S=[s0;s1;s2];
x=A*S;
SNR=0;%SNRdB
noi=zeros（1,L）;
x_in=[a1a2]*[s1;s2];
noi=sqrt（0.5*10^（SNR/10））*（randn（1,L））;
form=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
%%MVDR
Rx0=x*x'/L;
rxd=zeros（M,1）;
iRx0=pinv（Rx0）;
w_mvdr=iRx0*a0/（a0'*iRx0*a0）;
%%opt
R_in=x_in*x_in'/L;
iR_in=pinv（R_in）;
w_opt=（iR_in*a0）/（a0'*iR_in*a0）;
%%lms
Con=[a0a1a2];
Ncon=size（Con,2）;
f=[1;0;0];
B=zeros（M,M-Ncon）;
D=Con';
A=D（:
1:
Ncon）;
form=Ncon+1:
M
B（1:
Ncon,m-Ncon）=-A\D（:
m）;
B（m,m-Ncon）=1;
end
wq=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（yc）;
w_lms（:
1）=wq-B*wa（:
1）;
fork=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
%%beamforming
theta=-90:
1:
90;
N1=length（theta）;
y_mvdr=zeros（1,N1）;
y_opt=zeros（1,N1）;
y_lms=zeros（1,N1）;
fork=1:
N1
beta=2*pi*d*sin（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
%%result
figure,plot（theta,db（y_mvdr）,'-.'）
axis（[-9090-1000]）;
holdon
plot（theta,db（y_opt）,'r'）;
holdon
plot（theta,db（y_lms）,'g--'）;
holdon
gridon
title（'方向图'）
xlabel（'角度/°'）
ylabel（'归一化幅值/dB'）
legend（'MVDR','OPT','LMS'）;
实验二MISIC算法仿真
一实验目的
1、掌握MUSIC算法进行DOA估计的基本原理
2、了解影响DOA估计性能的因素
二实验内容
假设阵元个数M=12，信号方向分别为0度、2度和-35度，噪声为高斯噪声仿真：
（1）快拍次数一定，分辨力与SNR的关系
（2）快拍次数一定，分辨力与孔径的关系；（3）SNR一定，分辨力与快拍的关系。
三实验原理
窄带远场信号的DOA数学模型为：
阵列数据的协方差矩阵为：
由于信号与噪声相互独立，数据协方差矩阵可以分解为与信号、噪声相关的两部分，其中
是信号的协方差矩阵，
是信号部分。
对
进行特征分解有：
式中，
是由大特征值对应的特征矢量张成的子空间也即信号子空间，
而是由小特征值对应的特征矢量张成的子空间也即噪声子空间。
根据前面所述的性质2可知，在理想的条件下数据空间中的信号子空间与噪声子空间是相互正交的，即信号空间中的导向矢量也与噪声子空间正交
经典的MUSIC算法正是基于上述这个性质提出的，但考虑到实际接收数据矩阵是有限快拍数的，即数据协方差矩阵的最大似然估计为
对
进行特征分解可以计算得到噪声子空间特征矢量矩阵
。
由于噪声的存在，
与
并不能完全的正交。
因此，实际上求DOA是以最小优化搜索实现的，即
所以，MUSIC算法的谱估计公式为：
MUSIC算法的计算步骤如下：
（1）由阵列的接收数据得到数据协方差矩阵
；
（2）对
进行特征分解；
（3）由
的特征值进行信号源数的判断；
（4）确定信号子空间
与噪声子空间
；
（5）根据信号参数的范围进行谱峰搜索；
（6）找出极大值点对应的角度就是信号入射方向。
对于上述的MUSIC算法，还应该要注意以下几点：
对于非理想情况下得到协方差矩阵的特征值满足下式：
所以判断信号源数需要用到有关信号源估计的方法来进行信号源数的确定。
线阵的信号参数搜索范围为
，而面阵的搜索范围为
。
另外，一个确定阵列的导向矢量由阵元的位置唯一的确定；
MUSIC算法的一种归一化形式，即：
在实际应用中，对于一维导向矢量有下式成立：
四实验仿真及结果分析
（1）快拍次数一定，分辨力与SNR的关系
快拍次数L=1000，阵元间距为半波长，信噪比SNR在[-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度的信号被误判为一个信号，这是因为阵元间距太小，阵列的孔径就小，主瓣很宽，角度分辨率变差；随着阵元间距的增大，角度分辨率提高，当阵元间距大于0.3倍的波长，小于0.6倍的波长时，能很好的估计出三个信号的方向；当阵元间距大于0.6倍的波长时，出现了虚警信号，这是由于阵元间距太大，产生栅瓣。
（3）SNR一定，分辨力与快拍数的关系
信噪比为20dB，阵元间距为半波长，快拍数在[1,400]之间变化，得到的仿真结果如图3所示。
图3分辨率与快拍的关系总图（左）、分辨率与快拍的关系切面图（右）
由图3可知，当快拍数太少，少于50个时，会出现信号遗漏的现象，只有当快拍数足够多时，才能对三个信号做出正确的角度估计，快拍次数越多，角度分辨率越好。
附：
实验程序
实验一程序
clearall;
closeall;
clc
theta=-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;%NumberofSignals
w=[pi/6pi/10pi/8]';
xx=pi/180*2;yy=0;zz=-pi/180*35;
B=[1exp（-j*2*pi*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）;
1exp（-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*pi*d*sin（yy）/l）exp（-j*2*6*pi*d*sin（yy）/l）...
exp（-j*2*7*pi*d*sin（yy）/l）;
1exp（-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]）;
forii=1:
Lsnr
x=B*xxx+sigmaN（ii）*randn（8,N_x）+j*randn（8,N_x）*sigmaN（ii）;
R=x*x'/N_x;
[VD]=eig（R）;
[lambda,index]=sort（（diag（D）））;
UU=V（:
index（1:
5））;
fori=1:
length（theta）
AA=[1exp（-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*sin（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）;
end
Pmusic=10*log10（Pmusic）;
DOA（:
ii）=Pmusic';
end
aa=theta'*ones（1,Lsnr）;
figure
（1）
subplot（2,2,1）
plot（aa,DOA（:
10））
xlabel（'\theta/°'）
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,SNR）'）;
xlabel（'\theta/°'）
ylabel（'SNR/dB'）
zlabel（'PMUSIC'）
实验二程序
clearall;
closeall;
clc;
bili=0.1:
0.01:
1;
Lbili=length（bili）;
theta=-90:
0.1:
90;
Ltheta=length（theta）;
SNR=20
ps=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/6pi/10pi/8]';
xx=pi/180*2;yy=0;zz=-pi/180*35;
forii=1:
Lbili
sss=bili（ii）*2*pi;
B=[1exp（-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））;
1exp（-j*sss*sin（yy））exp（-j*2*sss*sin（yy））exp（-j*3*sss*sin（yy））...
exp（-j*4*sss*sin（yy））exp（-j*5*sss*sin（yy））exp（-j*6*sss*sin（yy））exp（-j*7*sss*sin（yy））;
1exp（-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;
[VD]=eig（R）;
[lambda,index]=sort（（diag（D）））;
UU=V（:
index（1:
5））;
fori=1:
length（theta）
AA=[1exp（-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*7*sss*sin（theta（i）/180*pi））];
WW=AA*UU*UU'*AA';
Pmusic（i）=abs（1/WW）;
end
Pmusic=10*log10（Pmusic）;
DOA（:
ii）=Pmusic';
end
aa=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）
升级会员 