1、課程設計報告實驗名稱：ESPRIT算法研究實驗日期：姓 名：學 號： 哈爾濱工業大學（威海）一、 設計任務實現空間譜估計算法，并考察算法性能。二、 方案設計1） 由均勻線陣形式，確定陣列的導向矢量；2） 由陣列導向矢量，對接收信號進行建模仿真；3） 由ESPRIT算法實現信號DOA估計；4） 考察算法性能與信噪比，采樣率，觀測時間等參數的關系。三、 設計原理3.1空間譜估計數學模型空間譜估計就是利用空間陣列實現空間信號的參數估計的一項專門技術。整個空間譜估計系統應該由三部分組成：空間信號入射、空間陣列接收及參數估計。相應地可分為三個空間，即目標空間、觀察空間及估計空間，也就是說空間譜估計系統由

2、這三個空間組成，其框圖見圖1。圖1 空間譜估計的系統結構對于上述的系統結構，作以下幾點說明。(1)目標空間是一個由信號源的參數與復雜環境參數張成的空間。對于空間譜估計系統，就是利用特定的一些方法從這個復雜的目標空間中估計出信號的未知參數。(2)觀察空間是利用空間按一定方式排列的陣元，來接收目標空間的輻射信號。由于環境的復雜性，所以接收數據中包括信號特征(方位、距離、極化等)和空間環境特征(噪聲、雜波、干擾等)。另外由于空間陣元的影響，接收數據中同樣也含有空間陣列的某些特征(互耦、通道不一致、頻帶不一致等)。這里的觀察空間是一個多維空間，即系統的接收數據是由多個通道組成，而傳統的時域處理方法通常

3、只有一個通道。特別需要指出的是：通道與陣元并不是一一對應，通道是由空間的一個、幾個或所有陣元合成的(可用加權或不加權)，當然空間某個特定的陣元可包含在不同的通道內。(3)估計空間是利用空間譜估計技術(包括陣列信號處理中的一些技術，如陣列校正、空域濾波等技術)從復雜的觀察數據中提取信號的特征參數。從系統框圖中可以清晰的看出，估計空間相當于是對目標空間的一個重構過程，這個重構的精度由眾多因素決定，如環境的復雜性、空間陣元間的互耦、通道不一致、頻帶不一致等。3.2 陣列信號處理首先，考慮N個遠場的窄帶信號入射到空間某陣列上，陣列天線由M個陣元組成，這里假設陣元數等于通道數，即各陣元接收到信號后經過各

4、自的傳輸信道送到處理器，也就是說處理器接收來自M個通道的數據。 (3.2-1)式中，是接受信號的幅度，是接收信號的相位，是接收信號的頻率。在窄帶遠場信號源的假設下，有 (3.2-2)根據式(3.2-1)和式(3.2-2)，顯然有下式成立： (3.2-3)則可以得到第L個陣元接收信號為 (3.2-4)式中，為第L個陣元對第i個信號的增益，表示第L個陣元在t時刻的噪聲，表示第i個信號到達第L個陣元時相對參考陣元的時延。將M個陣元在特定時刻接收的信號排列成一個列矢量，可得(3.2-5)在理想情況下，假設陣列中各陣元是各向同性的且不存在通道不一致、互耦等因素的影響，則式(3.2-4)中的增益可以省略(

5、即歸一化1)，在此假設下式(3.2-5)可以簡化為 (3.2-6)將式(3.2-6)寫成矢量形式如下： (3.2-7)式中，為陣列的維快拍數據矢量，為陣列的維噪聲數據矢量，為空間信號的維矢量，A為空間陣列的維流型矩陣(導向矢量陣)，且 (3.2-8)其中導向矢量 (3.2-9)式中，c為光速，為波長。由上述的知識可知，一旦知道陣元間的延遲表達式，就很容易得出待定空間陣列的導向矢量或陣列流型。下面推導一下空間陣元間的延遲表達式。假設空間任意兩個陣元，其中一個為參考陣元(位于原點)，另一個陣元的坐標為(x ,y, z)，兩陣元的幾何關系見圖，圖中“×”表示陣元。圖2 空間任意兩陣元的幾何

6、關系由幾何關系可以推導出兩陣元的波程差為 (3.2-10)這里的波程差其實就是位于x軸上兩陣元間的延遲、位于y軸上兩陣元間的延遲和位于z軸上兩陣元間的延遲之和。根據式(3.2-10)的結論，下面給出實際環境中常用的幾種陣列及陣元間的相互延遲表達式。(1)平面陣設陣元的位置為，以原點為參考點，另假設信號入射參數為，分別表示方位角與俯仰角，其中方位角表示與x軸的夾角。(2)線陣設 陣元的位置為，以原點為參考點，另假設信號入射參數為，表示方位角，其中方位角表示與y軸的夾角（即與線陣法線的夾角），則有 (3.2-11)(3)均勻圓陣 設以均勻圓陣的圓心為參考點，則有 (3.2-12)其中方位角表示與x

7、軸的夾角，r為圓半徑。3.3旋轉不變子空間算法原理3.3.1信號模型 算法介紹前，首先對信號進行建模。為了推導分析的方便，將波達方向的數學模型做如下理想狀態的假設：1) 陣列形式為線性均勻陣，陣元間距不大于信號波長的二分之一。2) 存生兩個完全相同的子陣，且兩個子陣的間距是己知的。 3) 噪聲序列為一零均值高斯過程，各陣元間噪聲相互獨立，噪聲與信號也相互獨立。4) 空間信號為零均值平穩隨機過程，通常為窄帶遠場信號。5) 信號源數小于子陣陣列元數，信號取樣數大于子陣陣列元數，以確保子陣陣列流型的各列線性獨立。6) 組成陣列的各傳感器為各向同性陣元，且無互耦以及通道不一致的干擾。圖3-1均勻線陣的

8、數學模型示意圖下圖給出了均勻線陣的數學模型示意圖:3.3.2 算法原理對于均勻線陣，相鄰子陣間存在一個固定間距，這個固定間距反映出各相鄰子陣間的一個固定關系，即子陣間的旋轉不變性，而ESPRIT算法正是利用了這個子陣間的旋轉不變性實現陣列的DOA估計。ESPRIT算法最基本的假設是存在兩個完全相同的子陣，且兩個子陣的間距是已知的。由于兩個子陣的結構完全相同，且子陣的陣元數為m，對于同一個信號而言，兩個子陣的輸出只有一個相位差，=1，2， N。下面假設第一個子陣的接收數據為，第二個子陣的接收數據為，根據前面所述的陣列模型可知 (3.1) (3.2)式中，子陣1的陣列流型=A，子陣2的陣列流型=

9、，且式中 （3.3）從上面的數學模型可知，需要求解的是信號的方向，而信號的方向信息包含在A和中，由于是一個對角陣，所以下面只考慮這個矩陣，即 （3.4）由上可知。只要得到兩個子陣間的旋轉不變關系，就可以方便地得到關于信號到達角的信息。下面的任務就是從式（3.1）和式（3.2）中得到兩個子陣間的關系。先將兩個子陣的模型進行合并，即 （3.5）在理想條件下，可得上式的協方差矩陣 (3.6) 對上式進行特征分解可得 （3.7）顯然上式中得到的特征值有如下關=，US為大特征值對應的特征矢量張成的信號子空間，為小特征值對應矢里張成的噪聲子空間。對于實際的快拍數據，式(3.7)應修正如下： (3.8) 由

10、前面的知識可知，上述的特征分解中大特征矢量張成的信號子空間與陣列流型張成的信號子空間是相等的。即 (3.9)此時，存在一個惟一的非奇異矩陣T，使得 （3.10）顯然，上述的結構對兩個子陣都成立，所以有 （3.11） 很顯然 ，由子陣1的大特征矢量張成的子空間、由子陣2的大特征矢量張成的子空間與陣列流型A張成的子空間三者相等，即 （3.12） 另外，由兩個子陣列在陣列流型上的關系可知 （3.13）再利用式（3.11）可知兩個子陣列的信號子空間的關系如下: （3.14）式（3.13）反映了兩個子陣列的陣列流型間的旋轉不變性，而式（3.14）反映了兩個子陣的陣列接收數據的信號子空間的旋轉不變性。如果

11、陣列流型A是滿秩矩陣，則由式（3.14）可以得到 （3.15）所以上式中的特征值組成的對角陣一定等于，而矩陣T的各列就是矩陣特征矢量。所以一旦得到上述的旋轉不變關系矩陣，就可以直接利用式（3.4）得到信號的入射角度。3.4 標準的旋轉不變子空間算法有上節的知識可知， ESPRIT算法的基本原理就是利用式（3.14）的旋轉不變性，常規的旋轉不變子空間算法就是利用上述的基本原理求解信號的入射角度信息。下面就分析解這個等式的兩種最經典、應用最廣泛方法:最小二乘（LS）法和總體最小二乘（TLS）法。3.4.1 最小二乘法由最小二乘的數學知識，我們知道式（3.14）的最小二乘解的方法等價于，約束條件 （

12、3.16）因此最小二乘法的基本思想就是使校正項盡可能小，而同時保證滿足約束條件。為了得到LS解，將式（3.14）代入式（3.16）即得 （3.17）對上式進行展開可得 = （3.18）上式對求導并令其等于0，可得 （3.19）上式的解顯然有兩種可能:(1) 當滿秩時，也就是子陣1的信號子空間的維數等于信號源數時，則上式的解是唯一的，可得上式的最小二乘解 （3.20） (2) 當不滿秩，即時，也就是信號源間存在相干或相差時，則存在很多解，但我們卻無法區別對應于方程的各個不同的解，可以稱這些解是不可辨識的，解的不可辨識性是我們需要解相干的原因所在。下面給出LS-ESPRIT算法的求解步驟：1由兩個

13、子陣的接收數據,，分別得到兩個子陣的數據協方差矩陣；2對矩陣對R, 進行特征分解，從而得到兩個數據矩陣的信號子空間和；3按式（3.20）得到矩陣，然后對其進行特征分解.得到N個特征值，就可得到對應的N個信號的到達角。當考慮嗓聲影響時，上述基于最小二乘算法的估計都是有偏的，這就是為什么需要考慮總體最小二乘ESPRIT算法的原因。3.4.2 總體最小二乘法 我們知道，普通最小二乘的基本思想是用一個范數平方為最小的擾動去于擾信號子空間，目的是校正中存在的嗓聲。顯然這就存在一個問題：如果同時擾動和，并使擾動范數的平方保持最小，是否可以同時校正和中存在的嗓聲?答案是肯定的，這就是總休最小二乘（TLS）的

14、思想。它考慮的是如下矩陣方程的解: (3.27)顯然上式可以改寫成 （3.28）所以TLS的解等價于 （3.29）定義如下一個矩陣，再結合上述分析過程。我們發現其實就是尋找一個的酉矩陣F，便得矩陣F與正交，也就說明了由F張成的空間與或列矢量張成的空間正交。所以矩陣F可從的特征分解中得到。因為 （3.30）式中的是由特征值構成的對角矩陣，E是與其相應的特征矢量構成的矩陣。即 （3.31）令是由對應特征值為0的特征矢量構成的矩陣.它屬子噪聲子空間，所以只要選擇矩陣F使之等于、，即可滿足上面提到的要求。即有 (3.32)可得 （3.33）如果令，則 （3.34）上式說明的特征值即是對角線元素。這說明

15、通過構造一個矩陣就可得到有關信號角度的信息.而這個矩陣的構造可通過式（3.30）得到，即 （3.35）下面直接給出TLS-ESPRIT算法的求解步驟：1由兩個子陣的接收數據，, 由式(3.8)得到數據協方差矩陣R；2通過矩陣對于的廣義特征分解，得到維數為的信號子空間；3由構造矩陣，并按式防（3.30）進行特征分解得到矩陣E，然后再按式(3.31)將矩陣分為四個小的矩陣；4按式（3.35）得到矩陣，然后對其進行特征分解,得到N個特征值，就可得到對應的N個信號的到達角。通過分析，我們可以得到標準ESPRIT算法的計算過程如下：(1)通過特征值或奇異值分解（EVD或SVD）分別估計兩個存在旋轉不變關

16、系的子陣的信號子空；(2)用上述的LS、TLS等方法求解式（3.14）所示的不變等式；(3)計算的特征值，其中如式（3.3）所示。然后利用式(3.4)求解人射信號的角度信息。就ESPRIT算法而言，TLS算法與LS算法性能基本一致，只是在低信噪比情況下TLS算法性能略好。四、 仿真結果 主要分析各個參數對估計誤差的影響，誤差函數定義如式(1)：4.1 信噪比 SNR對估計誤差的影響分析 首先對信噪比 SNR離散化取值，然后求得不同信噪比下的誤差，從而繪制出誤差隨信噪比改變的函數曲線如圖 2 所示，圖 2 中信噪比 SNR從- 15 取到 15，間隔為 1，運行次數為 100 次，其余條件如題中

17、所述。由圖 2 可知，隨著信噪比的增大，估計誤差會越來越小，即估計精度會越來越高。當待估計的信號方位角相差比較小時，估計的誤差也會相應的增大。另外，若兩信號為相干信號，則此方法將不能對其進行正確的估計。 4.2 陣元數 L對估計誤差的影響分析與 4.1節類似，首先對陣元數 L離散化取值，然后求得不同陣元數下的誤差， 從而繪制出誤差隨陣元數改變的函數曲線如圖 3 所示，圖 3 中陣元數從 K+1 取到K+25，間隔為 1，運行次數為 100 次，其余條件如題中所述。由于陣元數 L 需大于信號個數 K才能正確估計，故取值中含有信號個數 K。由圖 3 可知，隨著陣元數的增加，估計誤差會越來越小，即估

18、計精度會越來越高，但當陣元數大到一定程度后，對估計精度的影響則會慢慢的減小。 4.3 采樣點數 N對估計誤差的影響分析與 4.1節類似，首先對采樣點數 N離散化取值，然后求得不同采樣點數下的誤差，從而繪制出誤差隨采樣點數改變的函數曲線如圖 4 所示，圖 4 中采樣點數從 10 取到 200，間隔為 5，運行次數為 100 次，其余條件如題中所述。由圖 4可知，隨著采樣點數的增加，估計誤差會越來越小，即估計精度會越來越高。 估計誤差（角度）4.4 兩信號之間的角度差（GAP）對估計誤差的影響分析由于采用 ESPRI T算法對 DOA進行估計，若兩信號的方位距離較近時，雖然能得出估計結果，但估計的

19、精度會大受影響。因此，為了分析兩信號之間的不同間隔會對估計精度造成多大的影響，繪制不同 GAP下的估計誤差曲線如圖 5所示。處理方法與 4. 1 節類似，圖 5 中 GAP（單位為度）從 0. 1 取到 5，間隔為 0. 1，獨立運行次數為 100 次，其余條件如題中所述。由圖 5 可知，GAP越大估計越準確，但當 GAP大到一定程度后則估計精度趨于穩定。 4.5 單信號 DOA不同分布對估計誤差的影響分析信號波達方向（DOA）的取值區間為-90度到 90度，若只考慮只有一個信號的情況，則當信號的 DOA不同時，估計誤差也會不一樣。因此，為了分析不同的 DOA會對估計精度造成多大的影響，繪制不

20、同 DOA下的估計誤差曲線如圖 6所示。處理方法與 4. 1 節類似，圖 6 中 GAP從- 80 度取到 80 度，間隔為 5 度，獨立運行次數為 100 次，其余條件如題中所述。由圖 6 可知，DOA越靠近 0 度估計越準確，越靠近正負 90 度估計誤差越大。且仿真結果表明，當 DOA在正負 90 附近時，估計誤差太大，因此，為了不影響估計結果顯示效果，故在圖中未繪制正負 90 度附近的估計誤差。 2.6 減與不減噪聲方差(Rn)對估計誤差的影響分析由于有噪聲的影響，因此在估計信號自相關矩陣 R時，若將無信號時的自相關矩陣 Rn減去，即相當與減去估計出噪聲方差，則估計的精度會有所提高。結合

21、信噪比 SNR對估計誤差的影響，繪制減與不減噪聲方差兩種情況下估計誤差隨 SNR的變化曲線如圖 7所示，圖 7 中 SNR從-15dB到 5dB，間隔為 1dB，獨立運行次數為 100次。仿真結果表明，若減 Rn，主要是在低信噪比時對估計精度的改善較大，當信噪比較大時二者幾乎一樣。 五、 程序清單%本文件名為 drawTLSesprit.m %分析基于總體最小二乘的 ESPRIT算法（TLS-ESPRIT）的 DOA估計的性能%clear;clc;close all; %清除變量，清屏，關閉所有繪圖窗口% 調用格式:estimated,error=TLSesprit(p,L,K,SNR,DOA

22、);% 估計結果(弧度，矢量：p行 1列)：estimated% 估計誤差(弧度，標量：均方誤差)：error% 信號個數：p% 陣元數：L% 快拍數：K% 信噪比：SNR% 波達方向(弧度，矢量：p行 1列)：DOA% p=2; L=8; K=100; SNR=5; DOA=pi*(-10/180) pi*(20/180);%顯示估計結果%M=100; %設定獨立重復運行次數DOA=pi*(0/180) pi*(30/180); %波達方向(弧度，矢量：p行 1列)p=length(DOA); L=8; K=100; SNR=5; %參數設置，estimated,error=TLSespri

23、t(p,L,K,SNR,DOA); %函數調用polar(estimated,1 1,'r*'); %在極坐標中顯示估計結果(必須先轉化為弧度)h=title('');set(h,'string','TLS-ESPRIT: 估計值: ',num2str(estimated);h1=xlabel('');set(h1,'string','信號 DOA(度): ',num2str(DOA*180/pi);% %陣元數 L對估計誤差的影響分析% Ln=p+1:1:p+25; %陣元數 L需

24、大于信號個數 p才能正確估計% for n=1:length(Ln)% L=Ln(n);% for k=1:M% estimated,error=TLSesprit(p,L,K,SNR,DOA);% errorm(k)=error; %將每次的估計誤差存入變量 errorm中，便于求均值% end% errorn(n)=sum(errorm)/M; %求多次運行后的估計誤差的均值% end% figure(2);plot(Ln,errorn*180/pi,'r:*','LineWidth',2); %繪制曲線，并適當標注% xlabel('陣元數 L&#

25、39;);ylabel('估計誤差(° )');title('陣元數 L對估計誤差的影響');% % 結論：陣元數 L越大估計越準確，但當 L大到一定程度后則估計精度趨于穩定% %快拍數 K 對估計誤差的影響分析% Kn=10:10:200; %對快拍數離散化取值% for n=1:length(Kn)% K=Kn(n);% for k=1:M% estimated,error=TLSesprit(p,L,K,SNR,DOA);% errorm(k)=error; %將每次的估計誤差存入變量 errorm中，便于求均值% end% errorn1(n)

26、=sum(errorm)/M; %求多次運行后的估計誤差的均值% end% figure(3);plot(Kn,errorn1*180/pi,'r:*','LineWidth',2); %繪制曲線并適當標注% xlabel('快拍數 K');ylabel('估計誤差(° )');title('快拍數 K 對估計誤差的影響');% % 結論：快拍數 K越大估計越準確，但當 K 大到一定程度后則估計精度趨于穩定% %信噪比 SNR對估計誤差的影響分析% SNRn=-15:1:15; %對信噪比 SNR離散化取

27、值% for n=1:length(SNRn)% SNR=SNRn(n);% for k=1:M % estimated,error=TLSesprit(p,L,K,SNR,DOA);% errorm(k)=error; %將每次的估計誤差存入變量 errorm中，便于求均值% end% errorn(n)=sum(errorm)/M; %求多次運行后的估計誤差的均值% end% figure(4);plot(SNRn,errorn*180/pi,'r:*','LineWidth',2);%繪制曲線并適當標注(誤差：角度)% xlabel('SNR

28、9;);ylabel('估計誤差(° )');title('SNR對估計誤差的影響');% %結論：信噪比 SNR越大估計越準確，但當信噪比 SNR大到一定程度后則估計精度趨于 穩定%兩信號之間的角度差（GAP）的大小對估計誤差的影響分析%GAPn=0.1:0.1:5; %對兩信號之間的角度差（GAP）離散化取值for n=1:length(GAPn)GAP=GAPn(n); %每次循環只取其中一個值DOA=pi*(0/180) pi*(GAP/180);for k=1:M %M為獨立重復運行次數estimated,error=TLSesprit(p,

29、L,K,SNR,DOA);errorm(k)=error; %將每次的估計誤差存入變量 errorm中，便于求均值enderrorn(n)=sum(errorm)/M; %求多次運行后的估計誤差的均值endfigure(5);plot(GAPn,errorn*180/pi,'r:*','LineWidth',2);%繪制曲線并適當標注(誤差：角度) xlabel('GAP(° )');ylabel('估計誤差(° )');% title('兩信號之間的角度差（GAP）對估計誤差的影響');%結論

30、：GAP越大估計越準確，但當 GAP大到一定程度后則估計精度趨于穩定% %單個信號時，信號波達方向分布不同時對估計誤差的影響分析% DOAn=pi*(-80/180):(5/180):pi*(80/180); %對信號波達方向離散化取值（80度到 90度時誤差太大，因此未取）% for n=1:length(DOAn)% DOA=DOAn(n);p=1; %每次循環只取其中一個值,信號個數 p設為為 1% for k=1:M %M為獨立重復運行次數% estimated,error=TLSesprit1(p,L,K,SNR,DOA); %調用 TLSesprit1（一個信號的情況）% erro

31、rm(k)=error; %將每次的估計誤差存入變量 errorm中，便于求均值% end% errorn(n)=sum(errorm)/M; %求多次運行后的估計誤差的均值% end% figure(6);plot(DOAn*180/pi,errorn*180/pi,'b:*','LineWidth',2);%繪制曲線并適當標注(誤差：角度)% xlabel('DOA(° )');ylabel('估計誤差(° )'); % % title('DOA（單信號）不同分布對估計誤差的影響');% %

32、 %結論：越靠近 0度估計越準確，越靠近正負 90度估計誤差越大% %估計相關矩陣 R時，減與不減 Rn（無信號時的噪聲自相關矩陣）對估計誤差的 影響分析(結合信噪比 SNR對估計誤差的影響曲線)% SNRn=-15:1:5; %對信噪比 SNR離散化取值% for n=1:length(SNRn)% SNR=SNRn(n);% for k=1:M% estimated,error=TLSesprit(p,L,K,SNR,DOA);errorm(k)=error; %調用減Rn的函數% estimated,error=TLSespritRn(p,L,K,SNR,DOA);errormRn(k)

33、=error; %調用不減 Rn的函數% end% errorn(n)=sum(errorm)/M; errornRn(n)=sum(errormRn)/M;% end% figure(7);h=plot(SNRn,errorn*180/pi,SNRn,errornRn*180/pi,'r-.'); %繪制減與不減 Rn時的估計誤差曲線% legend('R=R-Rn','R'); set(h,'LineWidth',2); %用圖示在圖中標明哪條為減或不減 Rn的曲線% xlabel('SNR(dB)');ylabel('估計誤差(° )');% % title('減與不減 Rn對估計誤差的影響');% % %結論：若減 Rn，主要是在低信噪比時對估計精度的改善較大，當信噪比較大 時二者幾乎一樣%本文件名為 TLSesprit.m %基于總體最小二乘的 ESPRIT算法（TLS-ESPRIT）的 DOA估計函數% function estimated,error=TLSesprit(p,L,K,SNR,DOA)% 調用格式:estimated,error=TLSesprit(p,L,K,SNR,DOA);% 估計結果(