信號(hào)檢測(cè)與估計(jì)仿真作業(yè)

1、電 子 科 技 大 學(xué)信號(hào)檢測(cè)與估計(jì)計(jì)算機(jī)仿真作業(yè)2013年12月6日一 實(shí)驗(yàn)?zāi)康?.學(xué)習(xí)Matlab軟件在信號(hào)檢測(cè)與估計(jì)中的應(yīng)用2.學(xué)習(xí)MUSIC、ESPRIT、GEESE等的空間譜估計(jì)算法的原理，并通過仿真分析比較這三種算法的不同及性能特點(diǎn)3.通過仿真分析了解非平穩(wěn)噪聲和色噪聲對(duì)MUSIC、ESPRIT、GEESE方法性能的影響二 實(shí)驗(yàn)原理2.1最小錯(cuò)誤概率準(zhǔn)則出發(fā)點(diǎn)是如何使譯碼后的錯(cuò)誤概率PE為最小。其基本思路為：收到y(tǒng)j后，對(duì)于所有的后驗(yàn)概率P(x1|yj),P(x2|yj), ,P(xi|yj),，若其中P(x*|yj)具有最大值，則將x*判決為yj的估值。由于這種方法是通過尋找最大

2、后驗(yàn)概率來(lái)進(jìn)行譯碼的，故又常稱之為最大后驗(yàn)概率準(zhǔn)則。最大后驗(yàn)概率譯碼方法是理論上最優(yōu)的譯碼方法，但在實(shí)際譯碼時(shí)，既要知道先驗(yàn)概率又要知道后驗(yàn)概率，而后驗(yàn)概率的定量計(jì)算有時(shí)比較困難，需要尋找更為實(shí)際可行的譯碼準(zhǔn)則。2.2 MUSIC原理MUSIC算法是一種基于矩陣特征空間分解的方法。從幾何角度講，信號(hào)處理的觀測(cè)空間可以分解為信號(hào)子空間和噪聲子空間，顯然這兩個(gè)空間是正交的。信號(hào)子空間由陣列接收到的數(shù)據(jù)協(xié)方差矩陣中與信號(hào)對(duì)應(yīng)的特征向量組成，噪聲子空間則由協(xié)方差矩陣中所有最小特征值（噪聲方差）對(duì)應(yīng)的特征向量組成。MUSIC算法就是利用這兩個(gè)互補(bǔ)空間之間的正交特性來(lái)估計(jì)空間信號(hào)的方位。噪聲子空間的所有向

3、量被用來(lái)構(gòu)造譜，所有空間方位譜中的峰值位置對(duì)應(yīng)信號(hào)的來(lái)波方位。MUSIC算法大大提高了測(cè)向分辨率，同時(shí)適應(yīng)于任意形狀的天線陣列，但是原型MUSIC算法要求來(lái)波信號(hào)是不相干的。2.3 ESPRIT算法原理ESPRIT算法估計(jì)信號(hào)參數(shù)時(shí)要求陣列的幾何結(jié)構(gòu)存在所謂的不變性，這個(gè)不變性可以通過兩種手段得到：一是陣通過某些變換獲得兩個(gè)或兩個(gè)以上的相同子陣。由于這種算法在有效性和方面都有非常突出的表現(xiàn)，已經(jīng)被公認(rèn)為空間譜估計(jì)的一種經(jīng)典算法，隨著ESPRIT算法的深入研究，ESPRIT算法進(jìn)一步被廣大學(xué)者接受并推廣。2.4 GEESE算法基本原理信號(hào)子空間特征向量的廣義特征值法(GEESE)，可以在簡(jiǎn)化計(jì)算

4、的情況下解決ESPRIT算法中實(shí)際噪聲測(cè)量有誤差的問題。它利用信號(hào)子空間的一個(gè)顯著特征，那就是真實(shí)方向向量所張成的子空間與除了陣列輸出互相關(guān)矩陣的最小多重特征值之外的所有相應(yīng)特征向量所張成的子空間是一樣的。2.5非相關(guān)源的數(shù)學(xué)模型圖2.1 .陣元接收信號(hào)與位置的關(guān)系陣元接收信號(hào)與位置的關(guān)系如圖 2-1 所示。假設(shè)空間有M個(gè)陣元組成陣列，將陣元從1到M編號(hào)，并以陣元 1 作為參考點(diǎn)。由于各陣元無(wú)方向性，相對(duì)于基準(zhǔn)點(diǎn)的位置向量分別為令信源信號(hào)為 s(t),信號(hào)的載波為，則基準(zhǔn)點(diǎn)處的接收信號(hào)為,各陣元上的接收信號(hào)的表達(dá)式為式中k為波數(shù)向量，為入射信號(hào)傳播的傳播方向，單位向量，為信號(hào)相對(duì)于基準(zhǔn)點(diǎn)的延遲

5、時(shí)間，為信號(hào)傳播到基準(zhǔn)點(diǎn)r處的陣元相對(duì)于信號(hào)傳播到基準(zhǔn)點(diǎn)的滯后相位（弧度）。圖中為入射信號(hào)傳播方向角，k = k,。天線陣列中，信號(hào)的帶寬B一般比載波頻率小得多，所以就有，即信號(hào)在各陣元上的差異可以忽略不計(jì)，稱為窄帶信號(hào)。因此，陣列信號(hào)用向量形式可以表示為： 選定第一個(gè)陣元為基準(zhǔn)點(diǎn)，則方向向量為 式中，當(dāng)有P個(gè)信源時(shí)，波束的方向向量可分別表達(dá)為，。這P個(gè)方向向量組成的矩陣 稱為陣列的方向矩陣或響應(yīng)矩陣，它表示所有信源的方向信息。當(dāng)有N個(gè)窄帶信號(hào)入射到空間 M 個(gè)陣元上時(shí)候，接收的信號(hào)可以寫成如下的矢量形式： 式中，X(t)為陣列的M×1維快拍數(shù)據(jù)矢量，N(t)為陣列的 M ×

6、;1維噪聲數(shù)據(jù)矢量，S(t)為空間信號(hào)的N×1維矢，A空間陣列的M×N維響應(yīng)矩陣（導(dǎo)向矢量陣）。三實(shí)驗(yàn)過程3.1 實(shí)驗(yàn)1首先，當(dāng)M=1，實(shí)驗(yàn)1就變成了一般二元隨機(jī)振幅與隨機(jī)相位信號(hào)的波形檢測(cè)問題。這樣，就有如下模型：H0: xt=B1cos0t+0+nt, 0tTH1: xt=A1cos1t+0+nt, 0tT其中，噪聲n(t)是均值為零，功率譜密度為Pn=N02的高斯白噪聲；信號(hào)的振幅服從瑞麗分布，隨機(jī)相位服從均勻分布。由于功率信噪比早060dB范圍內(nèi)，取A0=100，N0=2，T=1us。根據(jù)教材的推導(dǎo)，得到錯(cuò)誤判決概率的表達(dá)式如下：PH0H1=N02N0+a2TPH1

7、H0=N02N0+a2T從上面可以看到，當(dāng)M=1時(shí)，模型很簡(jiǎn)單。當(dāng)M>1時(shí)，我們選取一個(gè)M值，用MATLAB軟件進(jìn)行仿真，得到仿真結(jié)果如下圖所示（MATLAB仿真程序見附錄1）：3.2 實(shí)驗(yàn)23.2.1. MUSIC算法過程及仿真1)收集陣列接收的數(shù)據(jù)樣本得到數(shù)據(jù)協(xié)方差矩陣 2)對(duì)協(xié)方差矩陣進(jìn)行特征分解，分解為特征值和特征向量的形式 式中，是對(duì)應(yīng)的特征向量所組成的矩陣，為的特征值。3)利用最小特征值的重?cái)?shù)M來(lái)估計(jì)信號(hào)源數(shù) 4)計(jì)算MUSIC 算法的功率譜 5)中極大值對(duì)應(yīng)的角度就是信號(hào)入射方向。在Matlab的命令窗口輸入仿真程序，見附錄2，得到如下結(jié)果：A = -10 23 45DOA

8、 = 44.9986 23.0019 -10.0519圖3.1 MUSIC算法非相關(guān)源的模擬測(cè)向仿真3.2.2 ESPRIT算法過程及仿真1)由快拍數(shù)據(jù)X可以得到數(shù)據(jù)相關(guān)矩陣R的估計(jì)2)對(duì)相關(guān)矩陣的估計(jì)做特征分解 （33）特征值矩陣，。3)利用小特征值的重?cái)?shù)估計(jì)得到信號(hào)源的估計(jì)數(shù)。4)將特征向量矩陣分解為子陣列矩陣得到： 5). 得到 6). 對(duì)進(jìn)行特征分解，得到K個(gè)特征值，就可以得到對(duì)應(yīng)K個(gè)信號(hào)的到達(dá)角。在Matlab的命令窗口輸入仿真程序，見附錄3，得到如下結(jié)果：A = -10 23 45DOA = 44.9986 23.0019 -10.0519圖3.2 ESPRIT算法非相關(guān)源的模擬測(cè)

9、向仿真3.2.3 GEESE算法：1)收集陣列接收的數(shù)據(jù)樣本得到數(shù)據(jù)協(xié)方差矩陣 2). 對(duì)進(jìn)行特征分解，得到和： 3). 選定J的值，將代入式中得到,: 4). 將，代入式得到；5). 將代入式中即得到信號(hào)源的波達(dá)方向: 在Matlab的命令窗口輸入仿真程序，見附錄4，得到如下結(jié)果：DOA = -30.0244 0.0000 30.0244 圖3.3 GEESE算法非相關(guān)源的模擬測(cè)向仿真3.2.4 三種算法性能比較MUSIC算法就是多重信號(hào)分類算法，它是一種信號(hào)參數(shù)估計(jì)算法，利用輸入信號(hào)協(xié)方差矩陣的特征結(jié)構(gòu)，給出的信息包括入射信號(hào)的數(shù)目、各個(gè)信號(hào)的波達(dá)方向、強(qiáng)度以及入射信號(hào)和噪聲間的互相關(guān)。E

10、SPRIT算法就是旋轉(zhuǎn)不變子空間算法，也是一種基于子空間的波達(dá)方向估計(jì)技術(shù)，與MUSIC算法不同的是，ESPRIT算法不需要精確知道陣列的方向向量，僅需各子需各子陣列之間保持一致，因此降低了對(duì)陣列校準(zhǔn)的嚴(yán)格性。GEESE算法是指信號(hào)子空間特征向量的廣義特征值法，可以在簡(jiǎn)化計(jì)算的情況下解決ESPRIT算法中實(shí)際噪聲測(cè)量有誤差的問題。它利用信號(hào)子空間的一個(gè)顯著特征，那就是真實(shí)方向向量所張成的子空間與除了陣列輸出互相關(guān)矩陣的最小多重特征值之外的所有相應(yīng)特征向量所張成的子空間是一樣的。這三種算法是空間譜估計(jì)中最經(jīng)典的算法。MUSIC算法估計(jì)值接近克拉美羅界算法(CRB)，對(duì)參數(shù)的少量偏差不敏感，更接近

11、實(shí)際應(yīng)用，具有較好的應(yīng)用前景,但需要對(duì)參數(shù)空間進(jìn)行搜索，計(jì)算量大。隨著信噪比的增加，MUSIC功率譜的峰值越高，估計(jì)精度越精確。在陣元數(shù)目不同，其他條件相同的情況下，陣元數(shù)目越大，旁瓣干擾越小，DOA估計(jì)越精確。在條件相同的情況下，相鄰信號(hào)（以50為例）的MUSIC功率譜隨著角度的增加而降低，信號(hào)源相關(guān)，MUSIC算法失效。色噪聲下，MUSIC算法方位估計(jì)不準(zhǔn)確。與MUSIC算法相比，ESPRIT算法還降低了計(jì)算量和存儲(chǔ)量，且避免了參數(shù)空間的搜索，計(jì)算量小于MUSIC算法，但是算法數(shù)據(jù)協(xié)方差矩陣中提取噪聲方差的估計(jì)，有時(shí)會(huì)使估計(jì)結(jié)果變壞，當(dāng)信號(hào)高度相關(guān)時(shí)估計(jì)性能同樣會(huì)變壞，且對(duì)所設(shè)的參數(shù)有較高

12、的要求，少量的誤差也會(huì)導(dǎo)致算法的失敗。在ESPRIT算法中隨著信噪比的增加，均方誤差越小，DOA估計(jì)效果越好。在陣元數(shù)目不同，其他條件相同的情況下，陣元數(shù)目越大，均方誤差越小，ESPRIT算法的估計(jì)精度越高。在條件相同的情況下，相鄰信號(hào)（以100為例）的均方差與信噪比關(guān)系隨著角度的增加而性能降低。ESPRIT 算法對(duì)相干信號(hào)的 DOA 估計(jì)失效。而GEESE算法，不僅計(jì)算量比較小，而且保證了算法的精度，但當(dāng)信號(hào)高度相關(guān)時(shí)性能仍然變壞。在GEESE算法中隨著信噪比的增加，均方誤差越小，DOA估計(jì)效果越好。在陣元數(shù)目不同，其他條件相同的情況下，陣元數(shù)目越大，均方誤差越小，GEESE 算法的估計(jì)精度

13、越高。在條件相同的情況下，相鄰信號(hào)（以50為例）的均方差與信噪比關(guān)系隨著角度的增加而性能降低，GEESE 算法對(duì)相干信號(hào)的DOA估計(jì)失效。3.3 實(shí)驗(yàn)33.3.1非平穩(wěn)噪聲下，三種算法對(duì)DOA估計(jì)影響的matlab仿真結(jié)果圖3.4 非平穩(wěn)噪聲下MUSIC算法對(duì)DOA估計(jì)的影響圖3.5非平穩(wěn)噪聲下ESPRIT算法對(duì)DOA估計(jì)的影響圖3.6非平穩(wěn)噪聲下GEESE算法對(duì)DOA估計(jì)的影響3.3.2色噪聲噪聲下，三種算法對(duì)DOA估計(jì)影響的matlab仿真結(jié)果圖3.7色噪聲下MUSIC算法對(duì)DOA估計(jì)的影響 圖3.8色噪聲下ESPRIT算法對(duì)DOA估計(jì)的影響圖3.9色噪聲下GEESE算法對(duì)DOA估計(jì)的影響

14、附錄1第1題仿真程序clc;clear;close allM=1 2 4 6 8 10; %M為接收機(jī)個(gè)數(shù)N=1 5 10; %N為射頻脈沖數(shù)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ì)算j個(gè)接收機(jī)和i個(gè)射頻脈沖數(shù)的最小錯(cuò)誤概率 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(' 估計(jì)誤差 '); title('16個(gè)接收單元最小誤判概率仿真'); legend('1個(gè)射頻脈沖數(shù)', '5個(gè)射頻脈沖數(shù)', '10個(gè)射頻脈沖數(shù)'); 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;%快拍數(shù)doa=-50 0 40 70/180*pi; %信號(hào)到達(dá)角w=1000 3000 2000 500'%信

22、號(hào)頻率M=8;%陣元數(shù)P=length(w); %信號(hào)個(gè)數(shù)，也可以用特征分解的大特征值數(shù)來(lái)決定l=150;%波長(zhǎng)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);%仿真信號(hào)(隨機(jī)信號(hào))s=10(snr/20)*sin(w*0:N-1);%仿真信號(hào)(正弦信號(hào))%s=10(snr/20)*exp(j*w*0:N-1);%仿真信號(hào)

23、（指數(shù)信號(hào)）%s=10(snr/20)*sawtooth(w*0:N-1);%仿真信號(hào)（鋸齒波信號(hào)）%s=2*exp(j*(w*1:N);%生成信號(hào)x=B*s;x=B*s+(1/sqrt(2)*(randn(M,N)+j*randn(M,N);%加了高斯白噪聲后的陣列接收信號(hào)x=B*s+awgn(x,snr);%加了高斯白噪聲后的陣列接收信號(hào)%v=1 1 1 1 1 1 1 1*randn;%Q=diag(v);%噪聲協(xié)方差矩陣對(duì)角線值相同%v=randn(1,M);%Q=diag(v);%噪聲協(xié)方差矩陣對(duì)角線值不同Q=randn(M);%噪聲協(xié)方差矩陣不為對(duì)角陣R=x*x'+Q; %

24、數(shù)據(jù)協(xié)方差矩陣%以下是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 非相關(guān)源測(cè)向仿真','fontsize',13,'fontweight','bold','fontname','隸書','color','red');grid on%MUSIC程序結(jié)束附錄3 ESPRIT算法程序N=1024;%快拍數(shù)doa=-10 23 45/180*pi; %信號(hào)到達(dá)角w=1000 3000 2000'%信號(hào)頻率M=8;

26、%陣元數(shù)P=length(w); %信號(hào)個(gè)數(shù)，也可以用特征分解的大特征值數(shù)來(lái)決定l=150;%波長(zhǎng)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);%仿真信號(hào)(隨機(jī)信號(hào))s=10(snr/20)*sin(w*0:N-1);%仿真信號(hào)(正弦信號(hào))%s=10(snr/20)*exp(j*w*0:N-1);%仿真信號(hào)（指數(shù)信號(hào)）%s=10(snr

27、/20)*sawtooth(w*0:N-1);%仿真信號(hào)（鋸齒波信號(hào)）%s=2*exp(j*(w*1:N);%生成信號(hào)%x=B*s;x=B*s+(1/sqrt(2)*(randn(M,N)+j*randn(M,N);x=B*s+awgn(x,snr);%加了高斯白噪聲后的陣列接收信號(hào)%v=1 1 1 1 1 1 1 1*randn;%Q=diag(v) %噪聲協(xié)方差矩陣對(duì)角線值相同v=randn(1,M);Q=diag(v);%噪聲協(xié)方差矩陣對(duì)角線值不同%Q=randn(M);%噪聲協(xié)方差矩陣不為對(duì)角陣R=x*x'+Q; %數(shù)據(jù)協(xié)方差矩陣%以下是ESPRIT程序Rxx=R(1:M-1,

28、1:M-1);%M-1維的自相關(guān)函數(shù)Rxy=R(1:M-1,2:M);%M-1維的互相關(guān)函數(shù)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值其對(duì)應(yīng)的角度即為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實(shí)際角度')ylabel('DOA估計(jì)角度')title('ESPRIT非相關(guān)源測(cè)向仿真','fontsize',13,'fontweight','bold','fontname' ,'color','red');grid on;%ESPRIT程序結(jié)束附錄4 GEESE算法程序N=1024;%快拍數(shù)doa=-30 0 30/180*pi; %信號(hào)到達(dá)角w=1000 3000 2000'%信號(hào)頻率M=8;%陣元數(shù)P=length(w); %信號(hào)個(gè)數(shù)