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); %信號個數