信号检测与估计仿真作业

1、电 子 科 技 大 学信号检测与估计计算机仿真作业2013年12月6日一 实验目的1.学习Matlab软件在信号检测与估计中的应用2.学习MUSIC、ESPRIT、GEESE等的空间谱估计算法的原理，并通过仿真分析比较这三种算法的不同及性能特点3.通过仿真分析了解非平稳噪声和色噪声对MUSIC、ESPRIT、GEESE方法性能的影响二 实验原理2.1最小错误概率准则出发点是如何使译码后的错误概率PE为最小。其基本思路为：收到yj后，对于所有的后验概率P(x1|yj),P(x2|yj), ,P(xi|yj),，若其中P(x*|yj)具有最大值，则将x*判决为yj的估值。由于这种方法是通过寻找最大

2、后验概率来进行译码的，故又常称之为最大后验概率准则。最大后验概率译码方法是理论上最优的译码方法，但在实际译码时，既要知道先验概率又要知道后验概率，而后验概率的定量计算有时比较困难，需要寻找更为实际可行的译码准则。2.2 MUSIC原理MUSIC算法是一种基于矩阵特征空间分解的方法。从几何角度讲，信号处理的观测空间可以分解为信号子空间和噪声子空间，显然这两个空间是正交的。信号子空间由阵列接收到的数据协方差矩阵中与信号对应的特征向量组成，噪声子空间则由协方差矩阵中所有最小特征值（噪声方差）对应的特征向量组成。MUSIC算法就是利用这两个互补空间之间的正交特性来估计空间信号的方位。噪声子空间的所有向

3、量被用来构造谱，所有空间方位谱中的峰值位置对应信号的来波方位。MUSIC算法大大提高了测向分辨率，同时适应于任意形状的天线阵列，但是原型MUSIC算法要求来波信号是不相干的。2.3 ESPRIT算法原理ESPRIT算法估计信号参数时要求阵列的几何结构存在所谓的不变性，这个不变性可以通过两种手段得到：一是阵通过某些变换获得两个或两个以上的相同子阵。由于这种算法在有效性和方面都有非常突出的表现，已经被公认为空间谱估计的一种经典算法，随着ESPRIT算法的深入研究，ESPRIT算法进一步被广大学者接受并推广。2.4 GEESE算法基本原理信号子空间特征向量的广义特征值法(GEESE)，可以在简化计算

4、的情况下解决ESPRIT算法中实际噪声测量有误差的问题。它利用信号子空间的一个显著特征，那就是真实方向向量所张成的子空间与除了阵列输出互相关矩阵的最小多重特征值之外的所有相应特征向量所张成的子空间是一样的。2.5非相关源的数学模型图2.1 .阵元接收信号与位置的关系阵元接收信号与位置的关系如图 2-1 所示。假设空间有M个阵元组成阵列，将阵元从1到M编号，并以阵元 1 作为参考点。由于各阵元无方向性，相对于基准点的位置向量分别为令信源信号为 s(t),信号的载波为，则基准点处的接收信号为,各阵元上的接收信号的表达式为式中k为波数向量，为入射信号传播的传播方向，单位向量，为信号相对于基准点的延迟

5、时间，为信号传播到基准点r处的阵元相对于信号传播到基准点的滞后相位（弧度）。图中为入射信号传播方向角，k = k,。天线阵列中，信号的带宽B一般比载波频率小得多，所以就有，即信号在各阵元上的差异可以忽略不计，称为窄带信号。因此，阵列信号用向量形式可以表示为： 选定第一个阵元为基准点，则方向向量为 式中，当有P个信源时，波束的方向向量可分别表达为，。这P个方向向量组成的矩阵 称为阵列的方向矩阵或响应矩阵，它表示所有信源的方向信息。当有N个窄带信号入射到空间 M 个阵元上时候，接收的信号可以写成如下的矢量形式： 式中，X(t)为阵列的M×1维快拍数据矢量，N(t)为阵列的 M ×

6、;1维噪声数据矢量，S(t)为空间信号的N×1维矢，A空间阵列的M×N维响应矩阵（导向矢量阵）。三实验过程3.1 实验1首先，当M=1，实验1就变成了一般二元随机振幅与随机相位信号的波形检测问题。这样，就有如下模型：H0: xt=B1cos0t+0+nt, 0tTH1: xt=A1cos1t+0+nt, 0tT其中，噪声n(t)是均值为零，功率谱密度为Pn=N02的高斯白噪声；信号的振幅服从瑞丽分布，随机相位服从均匀分布。由于功率信噪比早060dB范围内，取A0=100，N0=2，T=1us。根据教材的推导，得到错误判决概率的表达式如下：PH0H1=N02N0+a2TPH1

7、H0=N02N0+a2T从上面可以看到，当M=1时，模型很简单。当M>1时，我们选取一个M值，用MATLAB软件进行仿真，得到仿真结果如下图所示（MATLAB仿真程序见附录1）：3.2 实验23.2.1. MUSIC算法过程及仿真1)收集阵列接收的数据样本得到数据协方差矩阵 2)对协方差矩阵进行特征分解，分解为特征值和特征向量的形式 式中，是对应的特征向量所组成的矩阵，为的特征值。3)利用最小特征值的重数M来估计信号源数 4)计算MUSIC 算法的功率谱 5)中极大值对应的角度就是信号入射方向。在Matlab的命令窗口输入仿真程序，见附录2，得到如下结果：A = -10 23 45DOA

8、 = 44.9986 23.0019 -10.0519图3.1 MUSIC算法非相关源的模拟测向仿真3.2.2 ESPRIT算法过程及仿真1)由快拍数据X可以得到数据相关矩阵R的估计2)对相关矩阵的估计做特征分解 （33）特征值矩阵，。3)利用小特征值的重数估计得到信号源的估计数。4)将特征向量矩阵分解为子阵列矩阵得到： 5). 得到 6). 对进行特征分解，得到K个特征值，就可以得到对应K个信号的到达角。在Matlab的命令窗口输入仿真程序，见附录3，得到如下结果：A = -10 23 45DOA = 44.9986 23.0019 -10.0519图3.2 ESPRIT算法非相关源的模拟测

9、向仿真3.2.3 GEESE算法：1)收集阵列接收的数据样本得到数据协方差矩阵 2). 对进行特征分解，得到和： 3). 选定J的值，将代入式中得到,: 4). 将，代入式得到；5). 将代入式中即得到信号源的波达方向: 在Matlab的命令窗口输入仿真程序，见附录4，得到如下结果：DOA = -30.0244 0.0000 30.0244 图3.3 GEESE算法非相关源的模拟测向仿真3.2.4 三种算法性能比较MUSIC算法就是多重信号分类算法，它是一种信号参数估计算法，利用输入信号协方差矩阵的特征结构，给出的信息包括入射信号的数目、各个信号的波达方向、强度以及入射信号和噪声间的互相关。E

10、SPRIT算法就是旋转不变子空间算法，也是一种基于子空间的波达方向估计技术，与MUSIC算法不同的是，ESPRIT算法不需要精确知道阵列的方向向量，仅需各子需各子阵列之间保持一致，因此降低了对阵列校准的严格性。GEESE算法是指信号子空间特征向量的广义特征值法，可以在简化计算的情况下解决ESPRIT算法中实际噪声测量有误差的问题。它利用信号子空间的一个显著特征，那就是真实方向向量所张成的子空间与除了阵列输出互相关矩阵的最小多重特征值之外的所有相应特征向量所张成的子空间是一样的。这三种算法是空间谱估计中最经典的算法。MUSIC算法估计值接近克拉美罗界算法(CRB)，对参数的少量偏差不敏感，更接近

11、实际应用，具有较好的应用前景,但需要对参数空间进行搜索，计算量大。随着信噪比的增加，MUSIC功率谱的峰值越高，估计精度越精确。在阵元数目不同，其他条件相同的情况下，阵元数目越大，旁瓣干扰越小，DOA估计越精确。在条件相同的情况下，相邻信号（以50为例）的MUSIC功率谱随着角度的增加而降低，信号源相关，MUSIC算法失效。色噪声下，MUSIC算法方位估计不准确。与MUSIC算法相比，ESPRIT算法还降低了计算量和存储量，且避免了参数空间的搜索，计算量小于MUSIC算法，但是算法数据协方差矩阵中提取噪声方差的估计，有时会使估计结果变坏，当信号高度相关时估计性能同样会变坏，且对所设的参数有较高

12、的要求，少量的误差也会导致算法的失败。在ESPRIT算法中随着信噪比的增加，均方误差越小，DOA估计效果越好。在阵元数目不同，其他条件相同的情况下，阵元数目越大，均方误差越小，ESPRIT算法的估计精度越高。在条件相同的情况下，相邻信号（以100为例）的均方差与信噪比关系随着角度的增加而性能降低。ESPRIT 算法对相干信号的 DOA 估计失效。而GEESE算法，不仅计算量比较小，而且保证了算法的精度，但当信号高度相关时性能仍然变坏。在GEESE算法中随着信噪比的增加，均方误差越小，DOA估计效果越好。在阵元数目不同，其他条件相同的情况下，阵元数目越大，均方误差越小，GEESE 算法的估计精度

13、越高。在条件相同的情况下，相邻信号（以50为例）的均方差与信噪比关系随着角度的增加而性能降低，GEESE 算法对相干信号的DOA估计失效。3.3 实验33.3.1非平稳噪声下，三种算法对DOA估计影响的matlab仿真结果图3.4 非平稳噪声下MUSIC算法对DOA估计的影响图3.5非平稳噪声下ESPRIT算法对DOA估计的影响图3.6非平稳噪声下GEESE算法对DOA估计的影响3.3.2色噪声噪声下，三种算法对DOA估计影响的matlab仿真结果图3.7色噪声下MUSIC算法对DOA估计的影响 图3.8色噪声下ESPRIT算法对DOA估计的影响图3.9色噪声下GEESE算法对DOA估计的影响

14、附录1第1题仿真程序clc;clear;close allM=1 2 4 6 8 10; %M为接收机个数N=1 5 10; %N为射频脉冲数for i=1:1 Estimate=zeros(3,11); figure for j=1:3 for snr=0:60 all_err=0; for con=1:1100 Error=mreceiver(M(5),N(j),snr); %计算j个接收机和i个射频脉冲数的最小错误概率 all_err=all_err+sum(Error); end Estimate(j,snr+1)=all_err/2200.0; end end hold on plo

15、t(Estimate(1,:),'r'); plot(Estimate(2,:),'b-'); plot(Estimate(3,:),'m'); xlabel(' SNR =060dB '); ylabel(' 估计误差 '); title('16个接收单元最小误判概率仿真'); legend('1个射频脉冲数', '5个射频脉冲数', '10个射频脉冲数'); grid on hold off; endfunction da=anc_mac(x,y)

16、if(size(x)=size(y)error('the size of x is not equal with y'); endsu=0;a,b=size(x); for i=1:a for j=1:b su=su+x(i,j)*y(i,j); end end da=su;% m is the number of receiverfunction err=mreceiver(m,n,snr)S1=cos(0:2*pi/40:6*pi-2*pi/40);S2=cos(0:2*pi/10:(24*pi)-2*(pi/10);Gains=ones(m,1);a=raylrnd(Ga

17、ins);b=raylrnd(Gains);x=unifrnd(zeros(m,n),2*pi*ones(m,n);y=unifrnd(zeros(m,n),2*pi*ones(m,n);noise1=10(-snr/20)*wgn(m,120*n,0);noise2=10(-snr/20)*wgn(m,120*n,0);S_in1=zeros(m,120*n);S_in2=zeros(m,120*n);% signals receivedfor i=1:m for j=1:n S_in1(i,(j-1)*120+1:j*120)=a(i,1)*cos(0:2*pi/40:(120-1)*2*

18、pi/40)+x(i,j)*ones(1,120)+noise1(i,(j-1)*120+1:j*120); S_in2(i,(j-1)*120+1:j*120)=b(i,1)*cos(0:2*pi/10:(120-1)*2*pi/10)+y(i,j)*ones(1,120)+noise2(i,(j-1)*120+1:j*120); endendPcout11=zeros(m,1);Pcout12=zeros(m,1);Pcout21=zeros(m,1);Pcout22=zeros(m,1);% the structure of the macfor i=1:m for j=1:n Pcou

19、t11(i,1)=Pcout11(i,1)+anc_mac(S_in1(i,(j-1)*120+1:j*120),S1); Pcout12(i,1)=Pcout12(i,1)+anc_mac(S_in1(i,(j-1)*120+1:j*120),S2); Pcout21(i,1)=Pcout21(i,1)+anc_mac(S_in2(i,(j-1)*120+1:j*120),S1); Pcout22(i,1)=Pcout22(i,1)+anc_mac(S_in2(i,(j-1)*120+1:j*120),S2); endend % get the absolute value of the m

20、acPcout11=abs(Pcout11);Pcout12=abs(Pcout12);Pcout21=abs(Pcout21);Pcout22=abs(Pcout22);for i=1:m Pcout11(i,1)=Pcout11(i,1)*Pcout11(i,1)/(2*a(i,1)*a(i,1)+2)-log(a(i,1)*a(i,1)+1); Pcout12(i,1)=Pcout12(i,1)*Pcout12(i,1)/(2*a(i,1)*a(i,1)+2)-log(a(i,1)*a(i,1)+1); Pcout21(i,1)=Pcout21(i,1)*Pcout21(i,1)/(2*

21、b(i,1)*b(i,1)+2)-log(b(i,1)*b(i,1)+1); Pcout22(i,1)=Pcout22(i,1)*Pcout22(i,1)/(2*b(i,1)*b(i,1)+2)-log(b(i,1)*b(i,1)+1);end% get the decesion of the detectionerr=sum(Pcout11)<sum(Pcout12),sum(Pcout21)>sum(Pcout22);附录2 MUSIC算法程序N=1024;%快拍数doa=-50 0 40 70/180*pi; %信号到达角w=1000 3000 2000 500'%信

22、号频率M=8;%阵元数P=length(w); %信号个数，也可以用特征分解的大特征值数来决定l=150;%波长d=l/2;%阵元间距snr=15;%信噪比%阵列流形矩阵B=zeros(P,M);for k=1:P B(k,:)=exp(-j*2*pi*d*sin(doa(k)/l*0:M-1);end B=B'%阵列流形矩阵%s=10(snr/20)*(1/sqrt(2)*(randn(4,N)+j*randn(4,N);%仿真信号(随机信号)s=10(snr/20)*sin(w*0:N-1);%仿真信号(正弦信号)%s=10(snr/20)*exp(j*w*0:N-1);%仿真信号

23、（指数信号）%s=10(snr/20)*sawtooth(w*0:N-1);%仿真信号（锯齿波信号）%s=2*exp(j*(w*1:N);%生成信号x=B*s;x=B*s+(1/sqrt(2)*(randn(M,N)+j*randn(M,N);%加了高斯白噪声后的阵列接收信号x=B*s+awgn(x,snr);%加了高斯白噪声后的阵列接收信号%v=1 1 1 1 1 1 1 1*randn;%Q=diag(v);%噪声协方差矩阵对角线值相同%v=randn(1,M);%Q=diag(v);%噪声协方差矩阵对角线值不同Q=randn(M);%噪声协方差矩阵不为对角阵R=x*x'+Q; %

24、数据协方差矩阵%以下是MUSIC的程序U,V=eig(R);UU=U(:,1:M-P);theta=-90:0.5:90;%谱峰搜索for ii=1:length(theta) AA=zeros(1,length(M); for jj=0:M-1 AA(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/l); end Pmusic(ii)=abs(1/(AA*UU*UU'*AA');endPmusic=10*log10(Pmusic/max(Pmusic)+eps);plot(theta,Pmusic)xlabel('Angle

25、theta/degree')ylabel('P(theta) /dB')title('MUSIC 非相关源测向仿真','fontsize',13,'fontweight','bold','fontname','隶书','color','red');grid on%MUSIC程序结束附录3 ESPRIT算法程序N=1024;%快拍数doa=-10 23 45/180*pi; %信号到达角w=1000 3000 2000'%信号频率M=8;

26、%阵元数P=length(w); %信号个数，也可以用特征分解的大特征值数来决定l=150;%波长d=l/2;%阵元间距snr=15;%信噪比%阵列流形矩阵B=zeros(P,M);for k=1:P B(k,:)=exp(-j*2*pi*d*sin(doa(k)/l*0:M-1);endB=B'%s=10(snr/20)*(1/sqrt(2)*(randn(3,N)+j*randn(3,N);%仿真信号(随机信号)s=10(snr/20)*sin(w*0:N-1);%仿真信号(正弦信号)%s=10(snr/20)*exp(j*w*0:N-1);%仿真信号（指数信号）%s=10(snr

27、/20)*sawtooth(w*0:N-1);%仿真信号（锯齿波信号）%s=2*exp(j*(w*1:N);%生成信号%x=B*s;x=B*s+(1/sqrt(2)*(randn(M,N)+j*randn(M,N);x=B*s+awgn(x,snr);%加了高斯白噪声后的阵列接收信号%v=1 1 1 1 1 1 1 1*randn;%Q=diag(v) %噪声协方差矩阵对角线值相同v=randn(1,M);Q=diag(v);%噪声协方差矩阵对角线值不同%Q=randn(M);%噪声协方差矩阵不为对角阵R=x*x'+Q; %数据协方差矩阵%以下是ESPRIT程序Rxx=R(1:M-1,

28、1:M-1);%M-1维的自相关函数Rxy=R(1:M-1,2:M);%M-1维的互相关函数b=zeros(1,M-2);eye(M-2);b=b zeros(M-1,1);Cxx=Rxx-min(eig(Rxx)*eye(M-1);Cxy=Rxy-min(eig(Rxx)*b;a=eig(Cxx,Cxy);%找出最接近1的a值其对应的角度即为a1=abs(abs(a)-1);for i=1:P c,d=min(a1); a1(d)=1000; bb(i)=a(d); a(d)=1000;end%if P>1 disp('The angles of signals are'

29、;)%else %disp('The angle of signal is')%endA=doa*180/pidoa=A;DOA=asin(angle(bb)/pi)/pi*180DOA=sort(DOA);stem(doa,DOA,'r-')xlabel('DOA实际角度')ylabel('DOA估计角度')title('ESPRIT非相关源测向仿真','fontsize',13,'fontweight','bold','fontname' ,'color','red');grid on;%ESPRIT程序结束附录4 GEESE算法程序N=1024;%快拍数doa=-30 0 30/180*pi; %信号到达角w=1000 3000 2000'%信号频率M=8;%阵元数P=length(w); %信号个数