1.均匀线阵方向图
%8阵元均匀线阵方向图，来波方向为0度
clc;
clear all;
close all;
imag=sqrt(-1);
element_num=8;%阵元数为8
d_lamda=1/2;%阵元间距d与波长lamda的关系
theta=linspace(-pi/2,pi/2,200);
theta0=0;%来波方向
w=exp(imag*2*pi*d_lamda*sin(theta0)*[0:element_num-1]');
for j=1:length(theta)
a=exp(imag*2*pi*d_lamda*sin(theta(j))*[0:element_num-1]');
p(j)=w'*a;
end
figure;
plot(theta,abs(p)),grid on
xlabel('theta/radian')
ylabel('amplitude')
title('8阵元均匀线阵方向图')
当来波方向为45度时，仿真图如下：
8阵元均匀线阵方向图如下，来波方向为0度，20log（dB）
随着阵元数的增加，波束宽度变窄，分辨力提高：仿真图如下：
2.波束宽度与波达方向及阵元数的关系
clc
clear all
close all
ima=sqrt(-1);
element_num1=16; %阵元数
element_num2=128;
element_num3=1024;
lamda=0.03; %波长为0.03米
d=1/2*lamda; %阵元间距与波长的关系
theta=0:0.5:90;
for j=1:length(theta);
fai(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_num1*d));
psi(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_num2*d));
beta(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_num3*d)); end
figure;
plot(theta,fai,'r',theta,psi,'b',theta,beta,'g'),grid on
xlabel('theta');
ylabel('Width in radians')
title('波束宽度与波达方向及阵元数的关系')
仿真图如下：
3. 当阵元间距 时，会出现栅瓣，导致空间模糊。仿真图如下：
4. 类似于时域滤波，天线方向图是最优权的傅立叶变换 仿真程序和仿真图如下： clc
clear all close all
ima=sqrt(-1);
element_num=32; %阵元数 source_num=1; %信源数
d_lamda=1/2; %阵元间距与波长的关系 theta=linspace(-pi/2,pi/2,200); theta0=0; %来波方向
w=exp(ima*2*pi*d_lamda*sin(theta0)*[0:element_num-1]'); for j=1:length(theta);
a=exp(ima*2*pi*d_lamda*sin(theta(j))*[0:element_num-1]'); p(j)=w'*a; end figure;
subplot(1,2,1)
plot(theta,abs(p)),grid on xlabel('theta/radian') ylabel('amplitude')
title('按定义的方向图')
/2d λ
>
pfft=fftshift(fft(w,128));
subplot(1,2,2)
plot(linspace(-pi/2,pi/2,128),abs(pfft)),grid on
xlabel('theta/radian')
ylabel('FFT_amplitude')
title('最优权的傅里叶变换')
5.%最大信噪比准则方向图和功率谱
clc;
clear all;
close all;
ima=sqrt(-1);
element_num=8; %阵元数为8
d_lamda=1/2; %间距为半波长
theta=-90:0.5:90; %范围
theta0=0; %来波方向
theta1=20; %干扰方向
L=512; %采样单元数
for i=1:L
amp0=10*randn(1);
amp1=200*randn(1);
ampn=1;
s(:,i)=amp0*exp(ima*2*pi*1/2*sin(theta0*pi/180)*[0:element_num-1]');
j(:,i)=amp1*exp(ima*2*pi*1/2*sin(theta1*pi/180)*[0:element_num-1]');
n(:,i)=ampn*(randn(element_num,1)+ima*randn(element_num,1));
end
Rs=1/L*s*s'; %信号自相关矩阵
Rnj=1/L*(j*j'+n*n'); %干扰+噪声的自相关矩阵
[V,D]=eig(Rs,Rnj); %（Rs,Rnj）的广义特征值和特征向量[D,I]=sort(diag(D)); %排序
Wopt=V(:,I(8)); %最优权矢量
for j=1:length(theta)
a=exp(ima*2*pi*d_lamda*sin(theta(j)*pi/180)*[0:element_num-1]');
f(j)=Wopt'*a;
p(j)=a'*Rs*a+a'*Rnj*a;
end
F=20*log10(abs(f)/max(max(abs(f))));
P=20*log10(abs(p)/max(max(abs(p))));
subplot(121)
plot(theta,F);grid on;hold on
plot(theta0,-50:0,'.');plot(theta1,-50:0,'.')
xlabel('theta/0');ylabel('F in dB');
title('max-SNR 方向图');
axis([-90 90 -50 0]);
hold on
subplot(122)
plot(theta,P,'r');grid on
xlabel('theta/0');ylabel('功率in dB');
title('max-SNR功率谱')
仿真图如下：
6.
%ASC旁瓣相消----MSE准则
clc;close all;clear all
ima=sqrt(-1);
M=32; %辅助天线的数目
d_lamda=.5;
theta0=-30; %来波方向
theta1=60; %干扰方向
L=512; %采样单元数
s=zeros(1,512); %预划分一个区域
for ii=1:L;
amp0=1*randn(1); %信号的幅度随机产生，保证信号之间是不相关的
amp1=200*randn(1);
ampn=1;
jam(:,ii)=amp1*exp(ima*2*pi*0.5*sin(theta1*pi/180)*[0:M-1]')+ampn*(randn(M,1) +ima*randn(M,1)); %干扰+噪声
s(ii)=amp0*exp(ima*2*pi*0.5*sin(theta0*pi/180))+amp1*exp(ima*2*pi*0.5*sin(thet a1*pi/180))+ampn*(randn(1,1)+ima*randn(1,1));%接收信号（信号+干扰+噪声）s0(ii)=amp0*exp(ima*2*pi*0.5*sin(theta0*pi/180));
end
Rx=1/L*jam*jam'; %噪声自相关矩阵，相当于X（t）
r_xd=1/L*jam*s';
Wopt=pinv(Rx)*r_xd;
delta=s0-(s-Wopt'*jam);
delta1=abs(mean(delta.^2)-(mean(delta)).^2) %方差
theta=linspace(-pi/2,pi/2,200);
for jj=1:length(theta)
a=exp(ima*2*pi*.5*sin(theta(jj))*[0:M-1]');
f(jj)=Wopt'*a;
end
F=20*log10(abs(f)/(max(max(abs(f)))));
figure(1)
plot(theta*180/pi,F),grid on,hold on
plot(theta0,-50:0,'.')
plot(theta1,-50:0,'.')
xlabel('theta/o');
ylabel('F/dB');
title('MSE准则下的方向图')
axis([-90 90 -50 0]);%可为x轴和y轴设置一个极限范围
仿真图如下：
7. %线性约束最小方差(LCMV)准则
clc;
clear all ;
close all;
ima=sqrt(-1);
element_num=8; %阵元数
d_lamda=1/2; %阵元间距与波长的关系
theta=-90:0.5:90-0.3; %搜索范围
theta0=0; %三个信号源的来波方向
theta1=30;
theta2=60;
L=512; %采样单元数
for i=1:L;
amp0=10*randn(1);
amp1=100*randn(1);
amp2=10*randn(1);
ampn=10;
x(:,i)=amp0*exp(ima*2*pi*1/2*sin(theta0*pi/180)*[0:element_num-1]')+...
amp1*exp(ima*2*pi*1/2*sin(theta1*pi/180)*[0:element_num-1]')+...
amp2*exp(ima*2*pi*1/2*sin(theta2*pi/180)*[0:element_num-1]')+...
ampn*(randn(element_num,1)+ima*randn(element_num,1));
end
Rx=1/L*x* x';
steer1=exp(ima*2*pi*1/2*sin(theta0*pi/180)*[0:element_num-1]')
steer2=exp(ima*2*pi*1/2*sin(theta1*pi/180)*[0:element_num-1]')
steer3=exp(ima*2*pi*1/2*sin(theta2*pi/180)*[0:element_num-1]')
C=[steer1 steer2 steer3];
F=[1 0 1]'; %把三个方向都作为来波方向
w=inv(Rx)*C*(inv(C'*inv(Rx)*C))*F;
for j=1:length(theta);
a=exp(ima*2*pi*d_lamda*sin(theta(j)*pi/180)*[0:element_num-1]');
f(j)=w'*a;
p(j)=1/(a'*inv(Rx)*a);
end
f=10*log10(abs(f)/(max(max(abs(f)))));
figure(1)
subplot(121)
plot(theta,f),grid on,hold on
plot(theta0,-20:0,'.')
plot(theta1,-20:0,'.')
plot(theta2,-20:0,'.')
xlabel('theta/o');
ylabel('F/dB');
title('Capon beamforming方向图')
axis([-90 90 -20 0]);%可为x轴和y轴设置一个极限范围
P=10*log10(abs(p)/(max(max(abs(p)))));
subplot(122)
plot(theta,P),grid on,hold on
plot(theta0,-20:0,'.')
plot(theta1,-20:0,'.')
plot(theta2,-20:0,'.')
xlabel('theta/o');
ylabel('功率/dB');
title('Capon beamforming功率谱')
仿真图如下：
8. %Capon beamforming
Clc;
clear all ;
close all;
ima=sqrt(-1);
element_num=8; %阵元数
d_lamda=1/2; %阵元间距与波长的关系
theta=-90:0.5:90; %范围
theta0=0; %来波方向
theta1=20; %干扰方向
theta2=60; %干扰方向
L=1000; %采样单元数
for i=1:L;
amp0=10*randn(1);%信号的幅度随机产生，保证信号之间是不相关的amp1=200*randn(1);
amp2=200*randn(1);
ampn=3;
x(:,i)=amp0*exp(ima*2*pi*1/2*sin(theta0*pi/180)*[0:element_num-1]')+...
amp1*exp(ima*2*pi*1/2*sin(theta1*pi/180)*[0:element_num-1]')+...
amp2*exp(ima*2*pi*1/2*sin(theta2*pi/180)*[0:element_num-1]')+...
ampn*(randn(element_num,1)+ima*randn(element_num,1));
end
Rx=1/L*x* x';
R=inv(Rx);
steer=exp(ima*2*pi*1/2*sin(theta0*pi/180)*[0:element_num-1]');
w=R*steer/(steer'*R*steer);%Capon最优权矢量
for j=1:length(theta);
a=exp(ima*2*pi*d_lamda*sin(theta(j)*pi/180)*[0:element_num-1]');
f(j)=w'*a;
p(j)=1/(a'*R*a);
end
F=20*log10(abs(f)/(max(max(abs(f)))));
P=20*log10(abs(p)/(max(max(abs(p)))));%此处是功率的对数形式Figure;
subplot(121)
plot(theta,F),grid on,hold on
plot(theta0,-50:0,'.')
plot(theta1,-50:0,'.')
plot(theta2,-50:0,'.')
xlabel('theta/o');
ylabel('F/dB');
title('Capon beamforming方向图')
axis([-90 90 -50 0]);
subplot(122)
plot(theta,P),grid on
xlabel('theta/o');
ylabel('功率/dB');
title('Capon beamforming功率谱')
仿真图如下：
9.%不同方法估计协方差矩阵的Capon波束形成
clc;
clear all;
close all;
ima=sqrt(-1);
element_num=8; %阵元数为8
d_lamda=1/2; %间距为半波长
theta=-90:0.5:90; %范围
theta0=0; %来波方向
theta1=50; %干扰方向
L=512;%采样单元数
for i=1:L
amp0=10*randn(1);
amp1=50*randn(1);
ampn=0.5;
s(:,i)=amp0*exp(ima*2*pi*1/2*sin(theta0*pi/180)*[0:element_num-1]');
j(:,i)=amp1*exp(ima*2*pi*1/2*sin(theta1*pi/180)*[0:element_num-1]');
n(:,i)=ampn*exp(ima*2*pi*randn(1)*[0:element_num-1]');
end
Rx=1/L*(s+j+n)*(s+j+n)';%接收信号自相关矩阵
Rnj=1/L*(j+n)*(j+n)';%干拢+噪声的自相关矩阵
e=exp(ima*2*pi*d_lamda*sin(theta0*pi/180)*[0:element_num-1]');%来波方向信号Wopt_Rx=inv(Rx)*e/(e'*inv(Rx)*e);%采用接收信号的权矢量
Wopt_Rnj=inv(Rnj)*e/(e'*inv(Rnj)*e);%采用干拢+噪声信号的权矢量
for j=1:length(theta)
a=exp(ima*2*pi*d_lamda*sin(theta(j)*pi/180)*[0:element_num-1]');
f1(j)=Wopt_Rx'*a;
f2(j)=Wopt_Rnj'*a;
end
F1=20*log10(abs(f1)/max(max(abs(f1))));
F2=20*log10(abs(f2)/max(max(abs(f2))));
figure;
plot(theta,F1,theta,F2,'r'),grid on;hold on
plot(theta0,-50:0,'.');plot(theta1,-50:0,'.')
xlabel('theta/0');ylabel('F(1,2)/dB');
title('不同方法估计协方差矩阵的Capon波束形成');
axis([-90 90 -50 0]);
仿真图如下：
10．%多点约束的Capon波束形成和方向图
clc;
clear all ;
close all;
ima=sqrt(-1);
element_num=8; %阵元数
d_lamda=1/2; %阵元间距与波长的关系
theta=-90:0.3:90; %搜索范围
theta0=0; %来波方向
theta1=20; %干扰方向
theta2=50; %干扰方向
L=512; %采样单元数
Rx=zeros(element_num,element_num);%产生协方差矩阵
for i=1:L;
amp0=10*randn(1);
amp1=10*randn(1);
amp2=50*randn(1);
namp=0.5*randn(1);%噪声的幅度随机产生，保证噪声与信号之间是不相关的J(:,i)=amp1*exp(ima*2*pi*1/2*sin(theta1*pi/180)*[0:element_num-1]')+...
amp2*exp(ima*2*pi*1/2*sin(theta2*pi/180)*[0:element_num-1]')+...
namp*exp(ima*2*pi*randn(1)*[0:element_num-1]');
x(:,i)=amp0*exp(ima*2*pi*1/2*sin(theta0*pi/180)*[0:element_num-1]')+...
J(:,i);%表示接收信号
end
Rx=Rx+1/L*x*x';
R=inv(Rx);
w=amp0*exp(ima*2*pi*1/2*sin(theta0*pi/180)*[0:element_num-1]')+...
amp2*exp(ima*2*pi*1/2*sin(theta2*pi/180)*[0:element_num-1]')+...
amp1*exp(ima*2*pi*1/2*sin(theta1*pi/180)*[0:element_num-1]');
for j=1:length(theta);
a=exp(ima*2*pi*d_lamda*sin(theta(j)*pi/180)*[0:element_num-1]');
f(j)=w'*a;
p(j)=1/(a'*R*a);
end
F=10*log10(abs(f)/max(max(f)));%取对数的方向图
P=10*log10(abs(p)/max(max(p)));%取对数的功率谱估计
figure;
subplot(121)
plot(theta,F),grid on,hold on
plot(theta0,-50:0,'.');plot(theta1,-50:0,'.');plot(theta2,-50:0,'.');
xlabel('theta/o');
ylabel('F/dB');
title('capon beamforming方向图')
subplot(122)
plot(theta,P),grid on,hold on
plot(theta0,-50:0,'.');plot(theta1,-50:0,'.');plot(theta2,-50:0,'.');
xlabel('theta/o');
ylabel('P/dB');
title('capon beamforming功率谱')
仿真图如下：
大作业
%自适应波束形成方向图
clc;
clear all;
close all;
ima=sqrt(-1);
element_num=8; %阵元数
c=3e8;
f=500e6;
lamda=c/f;
d=1/2*lamda; %阵元间距与波长的关系
theta=-90:0.5:90; %范围
theta0=0; %来波方向
theta1=45; %干扰方向
theta2=60; %干扰方向
L=512; %采样单元数
for i=1:L;
amp0=10*randn(1);%信号的幅度随机产生，保证信号之间是不相关的amp1=100*randn(1);
amp2=100*randn(1);
ampn=10;
x(:,i)=amp0*exp(ima*2*pi*1/2*sin(theta0*pi/180)*[0:element_num-1]')+...
amp1*exp(ima*2*pi*1/2*sin(theta1*pi/180)*[0:element_num-1]')+...
amp2*exp(ima*2*pi*1/2*sin(theta2*pi/180)*[0:element_num-1]')+...
ampn*(randn(element_num,1)+ima*randn(element_num,1));%产生阵列信号end
Rx=1/L*x* x';
steer1=exp(ima*2*pi*1/2*sin(theta0*pi/180)*[0:element_num-1]')
steer2=exp(ima*2*pi*1/2*sin(theta1*pi/180)*[0:element_num-1]')
steer3=exp(ima*2*pi*1/2*sin(theta2*pi/180)*[0:element_num-1]')
C=[steer1 steer2 steer3];
F=[1 0 0]';%把0度方向作为来波方向
w=inv(Rx)*C*(inv(C'*inv(Rx)*C))*F;
for j=1:length(theta);
a=exp(ima*2*pi*1/2*sin(theta(j)*pi/180)*[0:element_num-1]');
f(j)=w'*a;
p(j)=1/(a'*inv(Rx)*a);
end
F=20*log10(abs(f)/(max(max(abs(f)))));
P=20*log10(abs(p)/(max(max(abs(p)))));%此处是功率的对数形式
figure
subplot(121)
plot(theta,P),grid on,hold on
xlabel('theta/o');
ylabel('F/dB');
title('阵列信号方向图')
axis([-90 90 -50 0]);
subplot(122)
plot(theta,F),grid on,hold on;
plot(theta0,-50:0,'.')
plot(theta1,-50:0,'.')
plot(theta2,-50:0,'.')
xlabel('theta/o');
ylabel('功率/dB');axis([-90 90 -50 0]);
title('自适应波束形成方向图')
仿真图如下：
当采样数L=2048时仿真图如下：
当采样数L=2048时仿真图如下：
当
L=5 12，幅度为
amp0=10*randn(1);
amp1=200*randn(1);
amp2=200*randn(1);
ampn=10;时仿真图如下：
当L=512，幅度为
amp0=100*randn(1); amp1=200*randn(1); amp2=200*randn(1); ampn=10;时仿真图如下：
多通道实时阵列信号处理系统的设计 由judyzhong于星期三, 11/11/2015 - 14:06 发表 作者：杨欣然,吴琼之,范秋香来源：电子科技 摘要：以全数字化信号产生和数字波束形成处理为基础的数字化阵列雷达已成为当代相控阵雷达技术发展的一个重要趋势，本文针对现代数字化阵列雷达对多通道数据采集和实时处理的需求，设计了一种基于FPGA的多通道实时阵列信号处理系统。可完成对20通道的中频数据采集，实时波束合成和数据传输功能，实验结果表明系统工作稳定、性能良好，具有良好的信噪比和通道一致性。 随着数字信号处理技术的不断进步和相应处理能力的不断提高，数字化阵列雷达以其动态范围大、扫描波束多、设计灵活性高的特点，已经逐步取代传统模拟阵列雷达成为高性能阵列雷达的主要研究方向。数字化阵列雷达需要将各个阵列天线接收的信号经过模拟下变频后经过AD采样并在数字域内进行信号处理，其典型的信号处理的方法包括数字波束形成(DBF)技术和波达方向估计(DOA)技术等。对于数字阵列雷达而言，对模拟下变频后的信号完成多
通道数据采集、数据处理和传输是系统的关键部分，对处理系统的同步性能、通道间幅相一致性均提出了很高要求，如进行DBF处理中通道间的不一致性将会影响波束合成后天线的方向图的特性，使增益下降、旁瓣电平升高。同时数字阵列雷达需要对各个通道采集下的数据实时地完成信号处理和数据传输功能，对处理系统的实时处理能力和信号吞吐能力提出了一定挑战。 本文介绍了一种针对DBF处理的多通道阵列信号处理系统的设计方案，以Xilinx的 XC7K325T FPGA为核心，完成了20通道的中频数据采集，并在FPGA内完成数字波束合成功能，可以同时完成8个波束指向的合成，并且将合成后的数据通过RapidIO结果传输至实时处理机进行进一步的处理。 1 系统设计原理和组成 本系统由数据采集模块和波束合成与传输模块两部分组成，系统框图如图1所示。 20路模拟中频输入通过SSMC连接器输入，由5片四通道A/D芯片AD9653采样后，通过LVDS接口串行输出到到FPGA的ISERDES输入模块中完成串并转换，并在FPGA内完
西工大《阵列信号处理》复习考点整理 考试形式： 一、8道问答题，每道题5分； 二、六道大题，包括PPT 上老师给出的那一道。 一 1. 均匀线列阵在波束扫描时，波束图怎么变化？ 当波束指向法线方向时，波束图具有最窄的主瓣宽度；随着阵元指向逐渐远离法线方向，主瓣一直指向所调方向并且展宽；除了指向法线方向外，主瓣都关于波束倾角轴不对称；当达到某一临界角时不能形成波束，但是在端射方向又可以形成波束。且在端射方向形成一个较宽的主瓣。 2．DI 是什么？ DI 表示指向性指数，其表达式为 D 为方向性，是阵列和孔径的一个常用性能度量。 ???=ππ φθθφθπφθ200 ),(sin 41) ,(P d d P D T T 3. DC 加权的特点 （1）旁瓣级给定时，主瓣宽度最小； （2）主瓣宽度给定时，旁瓣级最低； （3）等旁瓣级。 4. 频域快拍模型是什么，步骤是什么，常用的频域快拍取的时间有什么关系？ （1）记住《最优阵列处理技术》245页图 5.1 （2）步骤： ①把总的观测时间T 分为K 个不重叠的时间区域，区域长度为△T ； ②对时域快拍进行FT ； ③对频域向量（频域快拍）进行窄带波束形成； ④对上述频域信号进行IFT 。 （3）△T 的选择准则 ①△T 必须远大于平面波通过阵列的传播时间； ②△T 依赖于输入信号的带宽和信号的时域谱，16≥??T B （B*△T 足够大，选用频域快拍模型）。 5. 什么是均匀阵的瑞利限？ 常规波束形成分辨率的极限。表达式为 6. 空间白噪声的阵增益的相关计算。 阵列增益ωA 的定义为阵列的输出SNR 和一个阵元上的输入SNR 的比值。下标“ω”表示空域不相关的噪声输入。表达式如下：
信号子空间： 设N 元阵接收p 个信源，则其信号模型为：()()()()1 p i i i i x t s t a N t θ== +∑ 在无噪声条件下，()()()()()12,, ,P x t span a a a θθθ∈ 称()()()()12 ,, ,P span a a a θθθ为信号子空间，是N 维线性空间中的P 维子空间，记为P N S 。P N S 的正交补空间称为噪声子空间，记为N P N N -。 正交投影 设子空间m S R ∈，如果线性变换P 满足， 则称线性变换 P 为正交投影。 导向矢量、阵列流形 设N 元阵接收p 个信源，则其信号模型为：()()()()1 p i i i i x t s t a N t θ== +∑，其中矢量()i i a θ称为 导向矢量，当改变空间角θ，使其在空间扫描，所形成的矩阵称为阵列流形，用符号 A 表示，即 (){|(0,2)}a A θθπ=∈ 波束形成 波束形成（空域滤波）技术与时间滤波相类似，是对采样数据作加权求和,以增强特定方向信号的功率，即 ()()()()H H y t W X t s t W a θ==，通过加权系数W 实现对θ的选择。 最大似然 已知一组服从某概率模型 ()f X θ的样本集12,, ,N X X X ，其中θ为参数集合，使条件概率 ()12,,,N f X X X θ最大的参数θ估计称为最大似然估计。 不同几何形态的阵列的阵列流形矢量计算问题 假设有P 个信源，N 元阵列，则先建立阵列的几何模型求第i 个信源的导向矢量()i i a θ 选择阵元中的一个作为第一阵元，其导向矢量()1 [1]i a θ= 然后根据阵列的几何模型求得其他各阵元与第一阵元之间的波程差 n ?，则确定其导向矢量 ()2j n i a e πλ θ? =
1.均匀线阵方向图 %8阵元均匀线阵方向图，来波方向为0度 clc; clear all; close all; imag=sqrt(-1); element_num=8;%阵元数为8 d_lamda=1/2;%阵元间距d与波长lamda的关系 theta=linspace(-pi/2,pi/2,200); theta0=0;%来波方向 w=exp(imag*2*pi*d_lamda*sin(theta0)*[0:element_num-1]'); for j=1:length(theta) a=exp(imag*2*pi*d_lamda*sin(theta(j))*[0:element_num-1]'); p(j)=w'*a; end figure; plot(theta,abs(p)),grid on xlabel('theta/radian') ylabel('amplitude') title('8阵元均匀线阵方向图') 当来波方向为45度时，仿真图如下：
8阵元均匀线阵方向图如下，来波方向为0度，20log（dB）
随着阵元数的增加，波束宽度变窄，分辨力提高：仿真图如下：
2.波束宽度与波达方向及阵元数的关系 clc clear all close all ima=sqrt(-1); element_num1=16; %阵元数 element_num2=128; element_num3=1024; lamda=0.03; %波长为0.03米 d=1/2*lamda; %阵元间距与波长的关系 theta=0:0.5:90; for j=1:length(theta); fai(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_num1*d)); psi(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_num2*d)); beta(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_num3*d)); end figure; plot(theta,fai,'r',theta,psi,'b',theta,beta,'g'),grid on xlabel('theta'); ylabel('Width in radians') title('波束宽度与波达方向及阵元数的关系') 仿真图如下：
宽带信号中的三种二维平面阵DOA估计
宽带信号中的三种二维平面阵DOA 估计 一． 背景 目前关于阵列窄带信号的高分辨算法已比较成熟，但是随着信号处理技术的发展，信号环境日趋复杂，信号形式多样，信号密度日渐增大，窄带阵列探测系统的确定逐渐显示出来。 由于宽带信号具有目标回波携带的信息量大，有利于目标探测、参量估计和目标特征提取等特点，在有源探测系统中越来越多地得到应用。而在无源探测系统中，利用目标辐射的宽带连续谱进行目标检测是有效发现目标的一种重要手段。 ISM 方法把宽带信号在频域分解为J 个窄带分量，然后在每一个子带上直接进行窄带处理。因为信号为调频信号，所以信号在时域的分段实际上就是频域的分段。将信号分解为窄带信号后，我们就可以利用窄带算法进行处理，最后将各个结果进行加权综合，即可得到最终的结果。 二维DOA 估计是阵列信号处理中的重要内容，通过二维DOA 估计可以得到信号源在平面中的角度信息。一般采用L 型、面阵和平行阵或矢量传感器实现二维参数的估计，多数有效的二维DOA 估计算法是在一维DOA 估计的基础上，直接针对空间二维谱提出的，如二维MUSIC 算法以及二维CAPON 算法等。这两种算法可以产生渐进无偏估计，但要在二维参数空间搜索谱峰，计算量相当大。而采用二维ROOT MUSIC 算法可以减小计算量，但是需要付出精度下降的代价。 本次报告将结合宽带信号和二维DOA 估计算法，进行相关的算法介绍和仿真。 二． 算法介绍 1. 接收信号模型： 图 1 平面阵列示意图 如图1所示，设平面阵元数为M ×N ，信源数为K 。信源的波达方向为11(,),,(,)k k θφθφ ， 第i 个阵元与参考阵元之间的波程差为： 2(cos sin sin sin cos )/i i i x y z βπφθφθθλ=++ 设子阵1沿x 轴的方向矩阵为x A ，而子阵2的每个阵元相对于参考阵元的波程差就等于子阵1的阵元的波程差加上2sin sin /d πφθλ，所以接收信号为
%----------对角加载（LSMI 和SMI）方向图----------------------- %总结：这种算法主要给出了一种对角加载值的计算方法，对误差具有一定的稳健性，研究发现 %当数据协方差矩阵中含有信号分量会影响算法的性能。 clearall;clearall;clc; ratio_d_and_w=0.5; N_array=20;%阵列数 N_signal=60;% 样本数 ang1=0*pi/180;%所需信号的方向 SNR=5;%信噪比 ASd=sqrt(10.^(SNR/10)); ang2=40*pi/180;%干扰信号的方向 INR=45;%干噪比 ASi=sqrt(10.^(INR/10)); Sd=ASd*(randn(1,N_signal)+i*randn(1,N_signal));%Sd为所需信号 Si=ASi*(randn(1,N_signal)+i*randn(1,N_signal));%Si为干扰信号 Ni=randn(N_array,N_signal)+i*randn(N_array,N_signal);%Ni内噪声 Desired_Array=zeros(N_array,N_signal); Interferential_Array=zeros(N_array,N_signal); for LL=1:N_signal Interferential_Array(:,LL)=Si(LL)*test(ang2,N_array,ratio_d_and_w).'; Desired_Array(:,LL)=Sd(LL)*test(ang1,N_array,ratio_d_and_w).'; end X=zeros(N_array,N_signal); X= Interferential_Array +Ni; Rx=X*X'/N_signal; mm=std(diag(Rx));%对角加载值的确定下限 %mm=trace(Rx)/N_array;%对角加载值的确定上限 R1=Rx+mm*eye(size(Rx)); R=inv(R1); A_est=test(ang1,N_array,ratio_d_and_w); C= A_est; w_SMI=R*C/(C'*R*C);%对角加载 w_LSMI=inv(Rx)*C/(C'*inv(Rx)*C);%普通的Capon算法
阵列信号处理课程2011年作业 第1题 假定半波长间隔均匀分布线列阵的阵元数N ＝16，若入射平面波为62.5Hz 的正弦信号，信号持续时间为0.4s ，系统采样频率为1kHz ，阵列加权方式为均匀加权。分别给出 1. 当平面波信号分别从0,10, 20, 30, 40, 50, 60, 70, 80, 90, 100度方向入射时，指向90度的波束形成器的输出序列。 2. 当平面波信号分别从0:1:180度方向入射时，指向90度的波束形成器的输出序列经过平方求和后的分贝数输出。（把所有181个输出绘制在同一幅图中） 1）仿真图 图一：所求角度入射信号输出序列三维表示 注： 1. θ为信号入射角度，取值从0度到100度，每10°为一个间隔； 2. T 为整个阵元采样时间，对于不同的入射角度，t 的取值范围不同； 3. 输出信号幅度表示所有阵元的求和输出幅度，为有噪声情况。 结论： 0.8 t 输出信号幅度
从图一可以看出：①从90°入射的信号输出序列没有得到衰减，而其它角度入射的都得到了衰减；②从100°入射的信号和从80°入射的信号输出序列关于90°方向是对称的；③整个阵列对噪声有很好的抑制作用。 图二：入射信号0°到50°的输出序列 图三：入射信号60°到100°的输出序列 结论： 从图二和图三可以看出：①图一的所有结论；②90°方向入射信号没 0.10.20.30.40.5 -0.1 -0.0500.05 0.1 t A m p l i t u d e 0。 0.10.20.30.40.5 -0.1 -0.0500.05 0.1 t A m p l i t u d e 10。 0.10.20.30.40.5 -0.2 -0.100.1 0.2 t A m p l i t u d e 20。 0.10.20.30.40.5 -0.1 -0.0500.05 0.1 t A m p l i t u d e 30。 0.10.20.30.40.5 -0.1 -0.0500.05 0.1 t A m p l i t u d e 40。 0.10.20.30.40.5 -0.1 -0.0500.05 0.1 t A m p l i t u d e 50。 0.10.2 0.30.40.5 t A m p l i t u d e 60。 0.1 0.20.30.4 t A m p l i t u d e 70。 0.1 0.20.30.4 t A m p l i t u d e 80。 0.1 0.20.30.4 t A m p l i t u d e 90。 0.1 0.20.30.4 t A m p l i t u d e 100。
阵列信号处理仿真作业 需要解决的问题： 使用优化算法（可以使用遗传算法）挑选旁瓣相消的阵元 要求：（只需要选择一个突破点即可） ①可以针对不同类型的干扰、连片杂波、地杂波或密集型干扰等进行优化 ②也可以考虑存在阵列误差 下面我们针对第一个突破点进行仿真： 一、基本原理 图1给出了一般阵元级部分自适应处理的框图，通常称为多旁瓣相消器。 1 N H x 图1. 多旁瓣相消器结构部分自适应处理框图 如图1所示，整个天线阵的阵元加导向矢量权及用于压低旁瓣的锥削，可得到主通道输出0()m t ，0()m t 的方向图就指向目标方向，而从天线阵中选出M 个阵元作自适应单元，自适应单元加权为H x W ，于是得到主通道输出 00()()H m t t =W X ，辅助通道输出为()H x t W Y 。所以整个自适应信号处理器的输出为 0()()()H x e t m t t =-W Y (1) 其中0()m t 表示为主通道的输出；12[,,,]T M y y y =Y L 为选取的辅助单元接收的信 号；12[,,,]H T x M w w w ***=W L 为自适应权值；()H x t W Y 为形成辅助通道的输出。 在最小均方误差的准则下，求出的自适应权值就演变成为一个优化问题 220min ()min ()()H x E e t E m t t ?????-?????? W Y (2) 得
1 0()()()()H H x E t t E t m t -* ????=????W Y Y Y (3) 为了保证目标信号不损失，应对自适应权作约束，约束条件是在目标信号的 方向上，阵列自适应处理的增益为一常数。即在权值调整过程中，无论权x W 怎 样变化，对有用信号的增益不变。这样在使自适应阵输出()e t 的均方值2 ()E e t ?? ?? 最小时，能最大限度地抑制干扰且不损失有用信号能量。图1中应用式( 2) 的无 约束优化方程显然不合理，它不能保证有用信号增益不变。即由于辅助通道中包含有用信号的能量，就会导致辅助阵元中目标信号分量与主通道中目标信号相减，引起目标信号对消，导致目标增益下降。对此，应设法阻塞目标信号进入辅助支路，避免信号对消现象。一种方法是在优化方程中加入单位增益约束，强制目标方向增益不变，这样的优化方程求解比较麻烦，而且不适合自适应单元自动选取的算法。 这里，采用信号阻塞矩阵来抑制目标信号，使目标信号能量不能进入辅助通道。 0=X JX (4) 其中12[,,,]T N x x x =X L 为阵列信号；01020(1)0[,,,]T N x x x -=X L ；J 为信号阻塞矩阵(N -1)×N 。 在一般情况下，信号到达角为0θ，则信号阻塞矩阵J 可表示为 00001exp(())001exp(())01exp(())0 001exp(())j j j j ?θ?θ?θ?θ??--??--=??--??--??J L L L L 其中002()sin()d π?θθλ =，当00θ=?时，阻塞矩阵变为 11000110 0000 11-??-????=? ?????-?? J L L L L L L L L 这样就得到整个自适应系统框图（图2）。 图2中上支路为主通道, 保证信号完全通过，其加权为导向矢量权和压低旁瓣的锥削。下支路为辅助支路，信号阻塞矩阵阻止信号能量通过，将N 个阵元信号变成N -1个信号（降维），然后由辅助单元选择通路并选择参与自适应的单元，加自适应权后与上支路信号相减，得到自适应输出。
宽带信号中的三种二维平面阵DOA估计宽带信号中的三种二维平面阵DOA估计
一． 背景 目前关于阵列窄带信号的高分辨算法已比较成熟，但是随着信号处理技术的发展，信号环境日趋复杂，信号形式多样，信号密度日渐增大，窄带阵列探测系统的确定逐渐显示出来。 由于宽带信号具有目标回波携带的信息量大，有利于目标探测、参量估计和目标特征提取等特点，在有源探测系统中越来越多地得到应用。而在无源探测系统中，利用目标辐射的宽带连续谱进行目标检测是有效发现目标的一种重要手段。 ISM 方法把宽带信号在频域分解为J 个窄带分量，然后在每一个子带上直接进行窄带处理。因为信号为调频信号，所以信号在时域的分段实际上就是频域的分段。将信号分解为窄带信号后，我们就可以利用窄带算法进行处理，最后将各个结果进行加权综合，即可得到最终的结果。 二维DOA 估计是阵列信号处理中的重要内容，通过二维DOA 估计可以得到信号源在平面中的角度信息。一般采用L 型、面阵和平行阵或矢量传感器实现二维参数的估计，多数有效的二维DOA 估计算法是在一维DOA 估计的基础上，直接针对空间二维谱提出的，如二维MUSIC 算法以及二维CAPON 算法等。这两种算法可以产生渐进无偏估计，但要在二维参数空间搜索谱峰，计算量相当大。而采用二维ROOT MUSIC 算法可以减小计算量，但是需要付出精度下降的代价。 本次报告将结合宽带信号和二维DOA 估计算法，进行相关的算法介绍和仿真。 二． 算法介绍 1. 接收信号模型： 图 1 平面阵列示意图 如图1所示，设平面阵元数为M ×N ，信源数为K 。信源的波达方向为11(,),,(,)k k θφθφ ， 第i 个阵元与参考阵元之间的波程差为： 2(cos sin sin sin cos )/i i i x y z βπφθφθθλ=++ 设子阵1沿x 轴的方向矩阵为x A ，而子阵2的每个阵元相对于参考阵元的波程差就等于子阵1的阵元的波程差加上2sin sin /d πφθλ，所以接收信号为
万方数据
软件时空量，各阵元噪声满足空时白噪声的假设条件，即： Ｅ［ｎ（ｔ）ｎＨ（ｆ—ｆ）】＿盯２８（ｒ）ｘ Ｅ［ｎ（ｔ）ｎ７（卜ｆ）】－０（６） 阵列输出向量的二阶统计量用其外积的统计平 均表示，称之为阵列相关矩阵（将观测向量零均值化 则得到协方差矩阵）。定义为： Ｒ＝Ｅ［ｘ（ｔ）ｘＨ（ｆ）】－ＡＲｓＡＨ＋仃２，（７） 式中Ｒ＝Ｅ［ｓ（ｔ）ｓ爿（ｆ）】为信号的相关矩阵。 相关矩阵是阵列处理的基础，对Ｒ进行特征分图２单目标ＭＵＳＩＣ法的空间谱 解，根据信号子空间和噪声子空间的正交性可以实现仿真参数：（１）单目标情况：目标为２００Ｈ：的单频高分辨的目标方位估计。易证，Ｒ＝ＲＨ，这说明阵列协方正弦信号，目标方位角为６０。，噪声为零均值的高斯白差矩阵属于Ｈｅｒｍｉｔｉａｎ矩阵，其特征值为正值。令特征噪声，仿真分析的快拍数为１２８。 值为ｈｉ（ｉ＝ｌ…２一Ｍ），对应的特征向量为斗ｉ（ｉ＝１…２一Ｍ），协（２）两目标情况：目标１和目标２均为２００Ｈ：的方差矩阵的特征分解可写成：单频正弦信号，目标方位角分别为３０。和４５。，噪声为Ｒ＝ＵＡＵＨ＝ｙ．缸∥，（８）零均值的高斯白噪声，仿真分析的快拍数为１２８。 式中ｕ：【ｕ。，ｕｉ＝：１，－－．，ＨＭ］为由特征向量组成的酉矩。。仿妻竺果：单目标情况如图２所示，两目标情况阵；Ａ＝ｄｉａｇ［＆，五，．．．，知］为特征值构成的对角矩阵。如图啬霎磊染说明：空间谱中的峰值的高度并不表明将Ｒ的特征值按降序排列，根据特征值的大小可相应方位上的信号强度。增加阵元个数可以提高目标以将特征向量分成两部分，Ｕｓ＝［Ｕ。，ｕ：，．．．，ｕｄ为前Ｐ个最分辨力。 大特征值对应的特征向量构成的酉矩阵，其张成的空 间称为信号子空间，Ｕ。＝［ｕ吣ｕ嵋…，ｕ嗣为后Ｍ—Ｐ个最小 特征值对应的特征向量构成的酉矩阵，其张成的空间 称为噪声子空间。假设信号相关矩阵Ｒ。＝Ｅ【Ｓ（ｔ）ＳＨ（ｔ）】 非奇异，即各信号非相干，可以证明阵列方向矩阵Ａ 和信号子空间张成的子空间相同。又因为ｕ＝［ｕ。，Ｕｄ为 酉矩阵，所以有ｕｓＨＵ＃Ｏ。 由此可以定义ＭＵＳＩＣ算法的空间谱为： 删２蔬丽１（９）对以上空间谱进行峰值搜索可以得到波达方向的估计６；，ｉ＝ｌ…２．．，Ｐ。 实际中，Ｒ是未知的，可以由观测的数据向量估计，估计式为 食＝专善ｚ（力ｘＨ（力‘１。’对食进行特征分解得到噪声子空间的估计，进而得到ＭＵＳＩＣ空间谱和波达方向的估计。 ２Ｍａｔｌａｂ计算机仿真 下面对上面讨论的ＭＵＳＩＣ算法用Ｍａｔｌａｂ做计算机仿真。假设阵列为９阵元的等距均匀线列阵，阵元间距为信号中心频率对应的半波长，用该线阵来分别处理单个目标和两个目标信号源同时出现的情况。 图３两目标ＭＵＳＩＣ法的空间谱 ３结论 通过对ＭＵＳＩＣ算法的分析，从理论和系统仿真两方面证明将此法用于确定目标方位角的实用价值。是一种有效的测量目标方位角的方法。ＭＵＳＩＣ法对所有的特征向量重新加了权．噪声特征向量的权值为１．而信号特征向量的权值为０。对到达阵列的当前中的许多重要参数，如入射信号的个数，信号的入射方位、强度、入射波前的相关性以及噪声或干扰的强度等等，ＭＵＳＩＣ法都可以给出渐近无偏的估计。对于本文所讨论的空间谱估计的问题。ＭＵＳＩＣ法给出的谱要平滑得多，而且在信号方向上峰值又非常尖锐。除去不能分辨强相关或相干信号以外，ＭＵＳＩＣ法的主要缺点在于在搜索过程中使用了所有的噪声特征向量．从而导致较大的计算量。 参考文献： ［１】Ｒ．０．Ｓｃｈｍｉｄｔ：Ｍｕｌｔｉｐｌｅｅｍｉｔｔｅｒｌｏｃａｔｉｏｎａｎｄｓｉｇｎａｌｐａｒａｍｅｔｅｒ（转２９２页１ ＠㈨同邮局订阮８２?９４６ ３６０，，Ｌ／＿３０３— 万方数据
如对您有帮助，请购买打赏，谢谢您！ 信号子空间： 设N 元阵接收p 个信源，则其信号模型为：()()()()1 p i i i i x t s t a N t θ==+∑ 在无噪声条件下，()()()()()12,,,P x t span a a a θθθ∈ 称()()()()12,, ,P span a a a θθθ为信号子空间，是N 维线性空间中的P 维子空间，记为 P N S 。 P N S 的正交补空间称为噪声子空间，记为N P N N -。 正交投影 设子空间m S R ∈，如果线性变换P 满足， 则称线性变换P 为正交投影。 导向矢量、阵列流形 设N 元阵接收p 个信源，则其信号模型为：()()()()1 p i i i i x t s t a N t θ== +∑， 其中矢量()i i a θ称为导向矢量，当改变空间角θ，使其在空间扫描，所形成的矩阵称为阵列流形，用符号A 表示，即(){|(0,2)}a A θθπ=∈ 波束形成 波束形成（空域滤波）技术与时间滤波相类似，是对采样数据作加权求和,以增强特定方向信号的功率，即()()()()H H y t W X t s t W a θ==，通过加权系数W 实现对θ的 选择。 最大似然 已知一组服从某概率模型() f X θ的样本集12,,,N X X X ，其中θ为参数集合，使条件概 率()12,, ,N f X X X θ最大的参数θ估计称为最大似然估计。 不同几何形态的阵列的阵列流形矢量计算问题 假设有P 个信源，N 元阵列，则先建立阵列的几何模型求第i 个信源的导向矢量()i i a θ 选择阵元中的一个作为第一阵元，其导向矢量()1[1]i a θ= 然后根据阵列的几何模型求得其他各阵元与第一阵元之间的波程差n ?，则确定其导向矢量
阵列信号处理中的DOA（窄带） 空域滤波 波束形成：主要研究信号发射/接收过程中的信号增强。 空间谱估计 空域参数估计：从而对目标进行定位/给空域滤波提供空域参数。 测向波达方向估计(DOA) 空间谱：输出功率P关于波达角θ的函数，P(θ). 延迟——相加法/经典波束形成器注，延迟相加法和CBF法本质相同，仅仅是CBF法的最优权向量是归一化了的。
1、传统法常规波束形成CBF/Bartlett波束形成器 常规波束形成（CBF：Conventional Beam Former） Capon最小方差法/Capon 波束形成器/MVDR波束形成器 最小方差无畸变响应（MVDR：minimum variance distortionless response）Root-MUSIC算法
多重信号分类法解相干的MUSIC算法（MUSIC） 基于波束空间的MUSIC算法 2、[object Object]
TAM 旋转不变子空间法 LS-ESPRIT （ESPRIT） TLS-ESPRIT 确定性最大似然法（DML：deterministic ML） 3、最大似然法 随机性最大似然法（SML：stochastic ML）
4、综合法：特性恢复与子空间法相结合的综合法，首先利用特征恢复方案区分多个信号，估计空间特征，进而采用子空间法确定波达方向 最大似然估计法是最优的方法，即便是在信噪比很低的环境下仍然具有良好的性能，但是通常计算量很大。同子空间方法不同的是，最大似然法在原信号为相关信号的情况下也能保持良好的性能。 阵列流形矩阵（导向矢量矩阵）只要确定了阵列各阵元之间的延迟τ，就可以很容易地得出一个特定阵列天线的阵列流形矩阵A。 传统的波达方向估计方法是基于波束形成和零波导引概念的，并没有利用接收信号向量的模型（或信号和噪声的统计特性）。知道阵列流形 A 以后，可以对阵列进行电子导引，利用电子导引可以把波束调整到任意方向上，从而寻找输出功率的峰值。 ①常规波束形成(CBF)法 CBF法，也称延迟—相加法/经典波束形成器法/傅里叶法/Bartlett波束形成法,是最简单的DOA 估计方法之一。这种算法是使波束形成器的输出功率相对于某个信号为最大。 （参考自：阵列信号处理中DOA估计及DBF技术研究_赵娜）注意：理解信号模型
1．(1)关于接收天线阵列的假设。接收阵列由位于空间已知坐标处的无源阵元按一定的形式排列而成。假设阵元的接收特性仅与其位置有关而与其尺寸无关（认为其是一个点），并且阵元都是全向阵元，增益均相等，相互之间的互耦忽略不计。阵元接收信号时将产生噪声，假设其为加性高斯白噪声，各阵元上的噪声相互统计独立，且噪声与信号是统计独立的。 (2)关于空间源信号的假设。假设空间信号的传播介质是均匀且各向同性的，这时空间信号在介质中按直线传播，同时又假设阵列处在空间信号辐射的远场中，所以空间源信号到达阵列时可以看做是一束平行的平面波，空间源信号到达阵列各阵元在时间上的不同延时，可由阵列的几何结构和空间波的来向所决定。空间波的来向在三维空间中常用仰角和方位角来表征。其次，在建立阵列信号模型时，还常常要区分空间源信号是窄带信号还是宽带信号。所谓窄带信号是指相对于信号（复信号）的载频而言，信号包络的带宽很窄（包络是慢变的），因此在同一时刻，该类信号对阵列各阵元的不同影响仅在于因其到达各阵元的波程不同而导致的相位差异。 2．自适应波束形成亦称空域滤波，是阵列处理的一个主要方面，逐步成为阵列信号处理的标志之一，其实质是通过对各阵元加权进行空域滤波，来达到增强期望信号、抑制干扰的目的；而且可以根据信号环境的变化自适应嘚改变各阵元的加权因子。虽然阵列天线的方向图是全方向的，但阵列的输出经过加权求和后，可以被调整到阵列接收的方向增益聚集在一个方向，相当于形成了一个波束，这就是波束形成的物理意义所在。波束形成技术的基本思想是：通过将各阵元输出进行加权求和，将天线阵列波束导向到一个方向上，对期望信号得到最大输出功率的导向位置即给出波达方向估计。 3. ULA ：()1exp(2sin ) exp(2(1)sin )T k k k d d j j M θπθπθλλ?? =---???? α L 阵：(,)[(,),(,)]T x y a a a θφθφθφ=，其中 2sin cos 2(1)sin cos (,)[1,...],,T j d j M x a e e πθφπθφθφ---= 2s i n s i n 22s i n s i n 2(1 ...(,)[,,,] j d j d j M T y a e e e πθφπθφπθφθφ----= 面阵: 12()()()M D D D ?? ?? ??=??????? ?x y x y x y A A A A A A A ，其中1 1 2 2 1 1 2 2 2cos sin /2cos sin /2cos sin /2(1)cos sin /2(1)cos sin /2(1)cos sin /111 K K K K j d j d j d x j d M j d M j d M e e e e e e πθφλπθφλπθφλ πθφλ πθφλπθφλ---------?? ????=? ? ???? A
第25卷第4期 杭州电子科技大学学报V ol.25,N o.4 2005年8月Journal of Hangzh ou Dianzi University Aug.2005 阵列信号处理中几种关键技术的研究 王文勇1,陆安南2 (1.杭州电子科技大学通信工程学院,浙江杭州310018; 2.中国电子科技集团第36研究所,浙江嘉兴314001) 收稿日期:2004-09-17 作者简介:王文勇(1976-),男,安徽凤阳人,在读研究生,信号处理. 摘要:阵列信号处理技术的工程化研究日益成为其走向实际应用的重要步骤。该文首先介绍了阵 列信号处理实验系统的硬件组成,并在此基础上,分析了两种阵列信号处理的关键技术———谱估 计法中M USIC 算法和空间零点预处理波束形成技术,最后文章通过该实验系统对此两种处理技术 的有效性进行了验证。实验结果表明:M USIC 算法具有良好的测向精度,而空间零点预处理波束 形成技术具有较强信号分离能力。 关键词:阵列信号处理;空间谱估计;数字波束形成 中图分类号:T N914.53 文献标识码:A 文章编号:1001-9146(2005)04-0016-03 0 引 言 阵列信号处理是现代信号处理的一个重要分支,其本质是利用空间分散排列的传感器阵列和多通道接收机来获取信号的时域和空域等多维信息,以达到检测信号和提取其参数的目的。迄今为止,阵列信号处理的应用范围已经涉及诸如雷达、声纳、导航等领域。阵列信号处理的主要内容可分为波束形成技术、零点技术及空间谱估计技术等方面,它们都是基于对信号进行空间采样的数据进行处理,因此这些技术是相互渗透和相互关联的。波束形成技术的主要目的是要使阵列天线方向图的主瓣指向所需的方向,零点技术的主要目的是使天线的零点对准干扰方方向,前者是提高阵列输出所需要信号的强度,后者是减小干扰信号的强度,实质上都是提高阵列输出的信噪比的问题。而空间谱估计则主要研究信号到达方向(DOA )的问题。若将这几种技术结合,就会对空域信号处理的性能有很大的提高。 1 系统组成 一般阵列信号处理实验系统主要由3个部分组成:阵列天线、多通道接收机、以及阵列信号处理器。如图1所示 。 图1 阵列信号处理实验系统硬件基本结构 该文采用的阵列天线为9单元均匀圆阵,阵列半径为0.56m 。多通道数字接收机完成下变频、D/A 以及数字信号的Hilbert 变换。阵列信号处理器是实验系统的核心部分,进行的计算包括空间谱估计、波束形成及信号分离。由于考虑阵列信号处理的计算量巨大,其中包括有复矩阵的计算,为了满足对计算速度和实时性的要求,在本系统中阵列信号处理器采用两片高速数字信号处理专用芯片(DSP ),一片用来实现空间谱估计,另一片用来实现波束形成和信号分离。处理器的结构原理如图2所示。 由图2可见,信号既可以通过PCI 接口进入处理器,也可以通过DSP 的链路口进入处理器。利用一片接口芯片PCI9054实现和PCI 总线的通信。数字信号处理专用芯片通过两片FIFO 挂到本地总线的数据线上,这样就把双向接口设计为两个单项接口,即在当DSP 需要传输数据时,先把数据放入FIFO 中,然后通知主机已有数据输出,主机在任意时刻读出数据,完成DSP 到主机的数据传输;反之亦然。这样
阵列信号处理中的DOA（窄带） /接收过程中的信号增强。 参数估计：从而对目标进行定位/给空域滤波提供空域参数。 θ的函数，P(θ). /经典波束形成器注，延迟相加法和CBF法本质相同，仅仅 是CBF法的最优权向量是归一化了的。 CBF/Bartlett波束形成器 CBF：Conventional Beam Former） 最小方差法/Capon波束形成器/MVDR波束形成器 MVDR：minimum variance distortionless response） Root-MUSIC算法 多重信号分类法解相干的MUSIC算法 （MUSIC） 基于波束空间的MUSIC算法 TAM 旋转不变子空间法LS-ESPRIT （ESPRIT） TLS-ESPRIT 确定性最大似然法（DML：deterministic ML） 随机性最大似然法（SML：stochastic ML） 最大似然估计法是最优的方法，即便是在信噪比很低的环境下仍然具有良好的性能，但是通常计算量很大。同子空间方法不同的是，最大似然法在原信号为相关信号的情况下也能保持良好的性能。 计算量小，不需 进行谱峰搜索
阵列流形矩阵（导向矢量矩阵）只要确定了阵列各阵元之间的延迟τ，就可以很容易地得出一个 传统的波达方向估计方法是基于波束形成和零波导引概念的，并没有利用接收信号向量的模型（或信号和噪声的统计特性）。知道阵列流形A以后，可以对阵列进行电子导引，利用电子导引可以把波束调整到任意方向上，从而寻找输出功率的峰值。 ①常规波束形成(CBF)法 CBF法，也称延迟—相加法/经典波束形成器法/傅里叶法/Bartlett波束形成法,是最简单的DOA 估计方法之一。这种算法是使波束形成器的输出功率相对于某个信号为最大。 （参考自：阵列信号处理中DOA估计及DBF技术研究_赵娜） 注意：上式中，导向矩阵A表示第K个天线阵元对N个不同的信号 示第i个信号s(i)在M个不同的天线上的附加权值。 将式（2.6）的阵元接收信号，写成矢量形式为： X(t)=AS(t)+N(t) 其中，X(t)为阵列的M×1维快拍数据矢量，N(t)为阵列的M×1维噪声数据矢量，S(t)为信号空
阵列信号处理第二次理论作业 1．设一个辐射源距接收阵列的距离为r 0（该距离远大于天线的孔径），天线阵由M 个感应器构成，辐射源辐射的功率为P s ，噪声的平均功率为P n 。设阵列的时延可以使阵列的主瓣与信号的传播方向匹配，且阵列的加权系数为1。 1）在信号源处，信噪比是多少？ 2）在感应器处，信噪比是多少？ 3）当阵列的主瓣方向与传播方向匹配时，阵列输出的信噪比是多少？ 4）如信号源在接收阵列附近，发射的平均功率为P t ，信号以球面波的方式传播，到物体后备物体反射，且仍以球面波的方式传播，被物体反射的信号P s 为入射功率的ρ倍，物体到阵列中心位置的距离为r 0，计算阵列输出信号的信噪比（阵列的最大方向与目标反射信号的传播方向相同）？ 答： 1） 在信号源处，信噪比为0s n P SNR P = 2） 在感应器处，因信号功率被均匀分配到各个方向上，2 4s s P P R π'= 故：0 2 4s n P SNR SNR P R π''== 3） 当阵列主瓣方向匹配时，感应器处信噪比被加强，加强倍数为阵列增益。其值参考题4. 2 4s n P SNR SNR SNR G G M P R π''=?=?=? 4） 近场时，位于阵列中心的感应器处接收到的信号功率为 2244s t P P P R R ρ ππ?'= =，噪声功率仍为n P 。因此，该出的信噪比为： () 224t n n P P SNR P P R ρπ?= =? 计算得到阵列的输出信噪比为，中心点处信噪比乘以阵列增益G '。 () 2 24t n P G SNR SNR G P R ρπ'??'=?= ?阵列 2．阵列的增益在频域可以表示为：
基于MATLAB的阵列信号处理仿真方法 摘要:介绍如何使用MA TLAB构建阵列信号处理系统模型，包括相千信号模型，幅度和相位误差模型，针对不同模型，实现协方差矩阵产生方法，波达方向估计的子空间方法，自适应波束合成器的权值求解算法和方向图、阵列增益等系统参数的仿真这些仿真模型和方法，对于各种复杂的阵列信号处理研究，具有重要的基础作用 Array Signal Processing Simulation Based on MATLAB (Department of Communication Engineering, Chongqing University, Chongqing 400030, China) Abstract: MATLAB is used to develop the system model of array signal processing, such as the correlated signals, the error of amplitude and phase. To deal with these different models, MATLAB offers great instructions and factions, which make the simulating considerably convenient. The key programs were proposed, which could simulate the covariance matrix estimation, estimate the DOA (direction of arrival) of signal based on subspace, calculate the weight vector of adaptive beamformer, and get some basic system parameters. Some examples show the detail procedure and the programs are useful to the array signal processing research. Furthermore, these methods take an important function to many other complicated si mul anon s. Key words: array signal; covariance matrix; DOA; beamforming 引言 阵列信号处理源于60年代相控阵天线技术[}t}，目前， 广泛用于雷达、通信、超声波、语音、水下探测等等小同领 域。阵列信号处理的对象，是山阵列天线接收的采样快拍数 据，而利用各种信号处理的算法和工具，主要实现两个目的 f Z-41:一是进行空间谱估计，通过对信号分析，确定信号波 达方向;二是进行波束介成，利用自和非自算法，得到阵列 加权矢量。就实现手段而言，早期相控阵中，采用模拟方式， 而目前则是在数字域实现，利用FPGA或者DSP这些硬件 平台，完成信号处理算法。阵列信号处理的算法研究，}I-.}I 通过仿真进行。在MATLAB中，如何建立仿真模型，如何 实现各种基本研究参数评价指标仿真，对于复杂的阵列信号 处理具有重要意义。 1阵列信号处理的信号模型 山于阵列信号处理对象都是阵列天线接收信号，所以， 小同应用领域的信号模型是相同或相似的。对于一个远场窄 带零均值的入射信号，其N阵元构成阵列如图1所示。如 果用单位方向矢量“表示信号来向，矢量P表示阵元坐标， 则第Z个阵元接收信号相对于原点信号的时间延迟为 （1）