esprit與music演算法性能

③ WCDMA數據卡的關鍵技術

空時處理技術通過在空間和時間上聯合進行信號處理可以非常有效地改善系統特性。隨著第三代移動通信系統對空中介面標準的支持以及軟體無線電的發展，空時處理技術必將融入自適應數據機中，從而達到優化系統設計的目的。採用空時處理的方法，系統的發送端或接收端使用多個天線，同時在空間和時間上處理信號，它所達到的效果是僅靠單個天線的單時間處理方法所不能實現的：可以在一個給定BER質量門限下，增加用戶數；在小區給定的用戶數下，改善BER特性；可以更有效地利用信號的發射功率等等。
1、空時處理方法
在單用戶的情況下，空時處理方法的分類如圖1所示。
圖1
由於移動台一般不適於用多天線接收，在基站採用多個天線進行發射分集，可以使移動台的接收效果和移動台用多個接收天線時的效果相比擬，所以本文主要圍繞基站的空時處理技術展開討論。
2、波束成形技術
波束成形技術（Beamforming, BF）可分為自適應渣慧波束成形、固定波束和切換波束成形技術。固定波束即天線的方向圖是固定的，把IS－95中的3個120°扇區分割即為固定波束。切換波束是對固定波束的擴展，將每個120°的扇區再分為多個更小的分區，每個分區有一固定波束，當用戶在一扇區內移動時，切換波束機制可自動將波束切換到包含最強信號的分區，但切換波束機制的致命弱點是不能區分理想信號和干擾信號。自適應波束成形器可依據用戶信號在空間傳播的不同路徑，最佳地形成方向圖，在不同到達方向上給予不同的天線增益，實時地形成窄波束對准用戶信號，而在其他方向盡量壓低旁瓣，採用指向性接收，從而提高系統的容量。由於移動台的移動性以及散射環境，基站接收到的信號的到達方向是時變的，使用自適應波束成形器可以將頻率相近但空間可分離的信號分離開，並跟蹤這些信號，調整天線陣的加權值，使天線陣的波束指向理想信號的方向。
自適應波束成形的關鍵技術是如何較精確地獲得信道參數呢？對於上行鏈路，根據形成波束所用的信息可以將波束成形技術分成以下3類。
(1)基於空間結構的BF
基於空間結構的BF如基於輸入信號到達方向的BF（DOB），包括3類：基於最大信干噪比（SINR）的BF；基於最大似然（ML）准則的BF；基於最小均方誤差（MMSE）准則的BF。多址干擾的抑制依賴於信號的到達方向如陪答（DOA），所以DOB中的一個重要部分是信號的DOA估計。DOA估計方法有離散付里葉變換、MVDR（Minimum Variance Distortionless Response）估計器、線性預測、最大包絡法（MEM）、ML濾波器以及可變特徵結構的方法，其中包亂舉括MUSIC（Multiple Signal Classification）和ESPRIT法（Estimation of Signal Parameters via Rotational Invariance Technique）。
(2)基於訓練序列的BF
基於訓練序列的BF即時間參考BF（TRB），適用於多徑豐富且信道特性連續變化的環境，根據演算法可以分為塊自適應演算法（BAA）和采樣自適應演算法（SAA）兩類。BAA演算法包括特徵濾波器（EF）法、Stanford法、最大比合並（MRC）法和第一維納濾波器解（FWFS）、第二維納濾波器解（SWFS）。SAA演算法包括最小均方（LMS）演算法、歸一最小均方（NLMS）演算法、遞歸最小平方（RLS）演算法和共軛梯度法（CGM）。TRB技術要求同步精確，當時延擴展小時可以得到較好的性能。
(3)基於信號結構的BF（SSBF）
基於信號結構的BF（SSBF）即利用接收信號的時間或空間結構和特性來構造BF，可利用SSBF需要存儲例如恆包絡調制信號的恆模（CM）特性、信號的周期平穩性或數字調制信號的FA（Finite Alphabet）特性等知識，這種BF方法可以應用於不同的傳播條件，但需要考慮收斂性和捕獲問題。
對於下行鏈路而言，不同的復用方式可採用不同的解決方法：TDD方式，由於上下行鏈路採用相同的頻率，在保證信道參數在相鄰的上下行數據幀中幾乎沒有變化的情況下可以直接利用上行估計得到的信道參數，但這只適用於慢速移動的系統；FDD方式，由於上下行鏈路的頻率間隔一般都大於相關帶寬，因此上下行的瞬時信道幾乎是不相關的，此時採用反饋信道是最好的方法。
需要強調的一點是發送機的波束成形技術和接收機的波束成形技術是截然不同的，接收波束成形可在每個接收機獨立實現而不會影響其他鏈路，而發送波束成形會改變對其他所有接收機的干擾，所以要在整個網路內部聯合使用發送波束成形技術。
3、接收分集
由於CDMA系統通常有較多的多址干擾分量，而天線陣可以去除M－1個（M為天線數）干擾的特性並不能明顯地改善接收機 的SINR，所以在一般情況下，更好的方法是利用接收分集的方法，估計接收信號的形式，並確定匹配濾波器的加權系數。接收分集技術中的分集天線其實是空間域內的分集合並器，而不是BF。對於寬頻CDMA信號，信號帶寬一般大於信道相干帶寬，所以在時間域採用RAKE接收機，將信號在空間/時間上利用各種合並准則進行合並，這就是所謂的2D－RAKE接收機。一般的合並方式有：選擇合並（SC）即選擇具有最大信號功率的多徑；最大比合並（MRC）即每一路有一加權，根據各支路信噪比（SNR）來分配加權的權重，SNR大的支路權重大，SNR小的支路權重小。當每個分離多徑上的干擾不相關時，MRC方法可使合並信號的SINR最大；等增益合並（EGC）即選擇每一路的加權值都相等；Wiener濾波（OPT）即無論多徑之間的干擾是否相關，均可抑制干擾並使合並器輸出端的SINR最大，因此Wiener濾波的方法要好於最大比合並法，又稱為優化合並。
在空間和時間上利用不同的合並准則可以對系統起到不同的改善效果，理論證明，在理想功率控制和理想信道估計的條件下，空時聯合域優化合並方式對系統性能的改善最好。
4、發送分集技術
當發送方不能獲得信道參數時，空時發送分集可改善前向鏈路性能，這種機制是將發送天線的空間分集轉化為接收機可以利用的其他形式的分集，如延遲發送分集和空時編碼技術。空時編碼技術是同時從空間和時間域考慮設計碼字，它的基本原理是在多個天線上同時發送信息比特流所產生的向量，利用發送天線所發送序列的正交性，用兩個發送天線、一個接收天線所獲得的分集增益與一個發送天線、兩個接收天線的MRC接收機的一樣。
根據是否需要從接收機到發射機的反饋電路，發送分集技術可以分為開環和閉環兩種類型，前者發射機不需要任何信道方面的知識。開環發送分集方式有空時發送分集（STTD）、正交發送分集（OTD）、時間切換發送分集（TSTD）、延遲發送分集（DTD）以及分層的空時處理和空時柵格編碼；閉環發送分集方式有選擇發送分集（STD）。發送分集各方式具體如下。
(1)正交發送分集（OTD）
經過編碼和交織後的數據分成兩個不同的子流在兩個不同的天線上同時發送。為保證正交性，這兩個子流所用的Walsh碼是不同的。
(2)時間切換發送分集（TSTD）
在某一時刻每個用戶只使用一個天線，使用偽隨機碼機制在兩個天線之間切換。
(3)選擇發送分集（STD）
由於在TSTD方式中，瞬時使用的發送天線並不一定能在接收端得到最大的信噪比，所以使用一個反饋電路來選擇能提供使接收端得到最大信噪比的天線。
(4)空時發送分集（STTD）
空時發送分集是按如圖2所示的方法將數據編碼之後在兩個天線上發送出去。
(5)延遲發送分集（DTD）
用多個天線在不同時刻發送同一原始數據信號的多個復本，人為地產生多徑。
(6)分層空時結構(Bell LAyered Space－Time architecture，BLAST）
首先將原始信息比特分解成n個並行的數據流(稱為層)，送入不同的編碼器，再將編碼器的輸出調制以後使用相同的Walsh碼通過不同的天線發送出去。接收機側使用一個BF（迫零或MMSE准則）來分離不同的編碼數據流，然後將數據送入不同的解碼器，解碼器的輸出再重新組合建立原始的信息比特流。由於在波束成形處理中，MMSE和迫零方法都沒有充分利用接收機天線陣的分集潛力，所以提出了改進方案將接收處理也進行分級。即首先使用ViterbiMLSE演算法譯出最強的信號，然後將該強信號從接收的天線信號中去除後再檢測第二強的信號，如此反復直到檢測出最弱的信號。
該機制中，層到天線的映射並不是固定的，而是每np個碼符號之後周期性地改變，如圖3所示。這種映射關系保證了這些數據流最大可能地在不同的天線上被發送出去。
(7)空時柵格編碼
根據秩准則和行列式准則設計碼字，使設計出的碼字得到最大分集增益和編碼增益。以四進制相移鍵控（QPSK）四狀態空時柵格編碼為例，假定使用兩根天線發射，則星座圖和格形如圖4所示。
最右邊的元素編號S1S2的涵義是：從第一根天線發射出去的字元為S1，從第二根天線發射的字元為S2。
5、結束語
空時處理技術已顯示出非常誘人的發展前景，第三代移動通信標准中也支持空時處理技術，標準的出台為我們繼續研究物理可實現的空時處理技術提供了可能性，但將此技術實用化還存在許多亟待解決的方法和技術問題，有待於我們進一步研究。

④ 什麼叫窄帶，寬頻，單音信號

信號帶寬遠小於中心頻率的是窄帶信號，帶寬能和中心頻率相比擬或著是遠大於中心頻率的信號是寬頻信號。
窄帶信號的功率集中在中心頻率附近，兩者的功率譜密度和頻譜密度圖灶啟有很大的差距。處理方法也有很大差距。
寬頻信號與窄帶信號是相對的，不滿足窄帶條件的信號就稱為寬頻信號。目前，窄帶信號的定義也不盡相同。若信號帶寬為B
，時寬為T，中心頻率為f0，則窄帶信號的定義有：
定義1：
B
<<
f
0，即相對帶寬B/f
0<<
1
,一般B/f
0<
0.1。
定義2：
2v/c<<1/TB,其中v
是陣列與目標的相對徑向運動速度，c
是信號在介質中的傳播速度。
定義3：(
N－1)
d/c<<1/B，其中N
是陣元數目，d
是陣元間距。
定義4：2
Bτθ
+1≈1
,其中τθ是整個陣列以及延遲線的沒租延時之和。
定義5：該信號空間協方差矩陣在沒有雜訊時的第二個特徵值小於雜訊功率。
「定義1」是對窄帶信號的直觀理解,同時也是窄帶實信號可有效表示為其復解析形式的充分條件,在很多文獻中均以該定義來區分信號是寬頻信號還是窄帶信號。
「定義2」是指在存在相對運動的系統中,在信號的持續時間T
內,相對於信號的距離解析度,目標沒有明顯的位移,此時信號可視為窄帶信號,否則信號就是寬頻信號。
「定義3」是指在陣列信號處理中,如果信號帶寬的倒數遠遠大於信號掠過陣列孔徑的最大傳播時間,就稱為窄帶信號,否則為寬頻信號。
「定義4」和「定義5」是從陣列采樣數據自相關矩陣特徵值的角度來定義窄帶信號的,窄帶情況下,陣列采樣數據自相關矩陣的大特徵值個數等於信號個數。
可見，枯辯兆窄帶信號定義的非一致性決定了寬頻信號定義的非絕對性，不同的處理場合，應使用不同的定義確定信號是否是寬頻的。
至於某些窄帶信號的分析結果與演算法不適合寬頻信號，原因也是有差別的，不能給出統一的解釋。例如，在DOA估計中，子空間方法中的MUSIC、ESPRIT
方法以及它們的改進演算法，都是基於空間信號是窄帶這一假設推導出的。因為當信號是窄帶時，信號的延時近似等於信號的相移。當接收到的信號為寬頻信號時，近似條件不滿足，演算法性能會嚴重下降。
雖然寬頻信號沒有官方或某一組織的嚴格定義,但對於超寬頻信號,
FCC
(
Federal
Communica2tions
Commission)
給出了嚴格的定義.「超寬頻」的定義是:信號的相對帶寬(信號頻譜的帶寬與其中心頻率之比)
大於等於20
%
,或者絕對帶寬大於等於500MHz。該定義沒有界定信號的時域波形特徵。

⑤ 從經典譜估計到現代譜估計

諧波分析最早可追溯到古代對時間的研究。18世紀伯努利（Bernoulli）、歐拉（Euler）和拉格朗日（Lagrange）等人對波動方程及其正弦解進行了研究，19世紀初葉，傅立葉（Fourier）證明了在有限時間段上定義的任何函數都可以用正弦和餘弦分量的無限諧波的總和來表示。1898年舒斯特（Schuster）以傅立葉分析為基礎來擬合待分析信號，研究太陽黑子數的周期變化，並提出了周期圖的概念。1930年維納（Wiener）發表了經典性論文《廣義諧波分析》，對平穩隨機過程的自相關函數和功率譜密度作了精確的定義，證明了二者之間存在著傅氏變換（以下簡稱傅氏變換）的關系，從而為功率譜分析奠定了堅實的統計學基礎。由於1934年辛欽（Khintchine）也獨立地證明了自相關函數和功率譜之間的傅氏變換關系，即維納-辛欽（Wiener-Khintchine）定理。根據這個定理（詳見第一章），平穩離散隨機信號x（n）的自相關函數r_xx（m）

r_xx（m）=E［x（m+n）x^*（n）］ （4-1）

與功率譜P_xx（e^jω）之間構成一傅氏變換對，即

地球物理信息處理基礎

若x（n）還是各態遍歷性的，則其自相關函數可由它的一個采樣時間序滲銷列用時間平均的方法求出，即

地球物理信息處理基礎

在大多數應用中x（n）是實信號，於是上式可寫成

地球物理信息處理基礎

實際上，一般只能在時域觀測到隨機信號的有限個采樣值（例如N個值），可表示為

x_N（n）=｛x（0），x（1），…，x（N-1）｝=｛x（n），n=0，1，…，N-1｝

其自相關函數只能由這N個采樣數據進行估計，常用有偏估計

地球物理信息處理基礎

這是一種漸近一致估計，稱之為采樣自相關函數。

用采樣自相關函數的傅氏變換作為功率譜的估計，這種方法是布萊克曼（R.Blackman）和杜基（J.Tukey）在1958年提出來的，稱為功率譜估計的自相關法（簡稱BT法）。此方法需要先求出有限個觀測數據估計自相關函數，然後再根據式（4-2）計算出功率譜。在快速傅氏變換（FFT）演算法提出之前，這是一種最流行的功率譜估計方法。

1965年庫利（Cooley）和杜基（Tukey）完善了著名的FFT演算法，使計算傅氏變換的速度提高了兩個數量級，運算量顯著降低，這樣DFT變換很快在各領域，特別是在工程實踐中得到了孝運廣泛應用。由式（4-5）知，

為x（n）與x（-n）的卷積運算，因為

地球物理信息處理基礎

若x（n）的傅氏變換為X（e^jω），則x（-n）的傅氏變換等於X^*（e^jω）。對式（4-5）兩端取傅氏變換，得到

地球物理信息處理基礎

這表明：通過對隨機數據直接進行離散傅氏變換，然後取其幅值的平方，再對多樣本進行此種運算並取平均值作為功率譜的估計，即舒斯特的周期圖，這種譜估計受到了人們的普遍重視巧喊梁，因為它不需要計算自相關函數，而直接計算功率譜。

周期圖和自相關法以及它們的改進方法稱為功率譜估計的經典方法，周期圖和自相關法是經典功率譜估計的兩個基本方法。由於FFT的出現，周期圖和自相關法往往被結合起來使用，其步驟如下：

（1）對x_N（n）補N個零，求

；

（2）由

作傅氏變換，得

，這時｜m｜≤M=N-1；

（3）對

加窗函數v（m），這時｜m｜≤M＜＜N-1，得

；

（4）利用

，求

的傅氏變換，即

地球物理信息處理基礎

由周期圖法得到的功率譜

，其估計方差並不隨樣本長度的增加而趨於零，

，出人意料的是，不管數據記錄有多長，周期圖和自相關法得到的估計都不是功率譜的良好估計。事實上，隨著記錄長度增加，這兩種估計的隨機起伏反而會更加嚴重！此外，它們存在著以下兩個難以克服的固有缺點。

（1）頻率解析度（區分兩個鄰近頻率分量的能力）不高。因為它們的頻率解析度（赫茲）反比於數據記錄長度（秒）（即Δf=k/T_p=k/NT，k為常數，T_p=NT為數據的記錄長度，T為采樣周期），而實際應用中一般不可能獲得很長的數據記錄，即觀察到的數據只能是有限個，而觀察不到的數據被認為是0。這樣，如果只有N個觀測數據，而對於N以外的數據，信號仍有較強的相關性，那麼估計出的功率譜就會出現很大的偏差。

（2）對於有限的觀測數據，相當於將信號在時域內乘以矩形窗函數，因而在頻域內則相當於使真正的功率譜與sinc函數進行卷積，由於sinc函數不同於δ 函數，它有主瓣和旁瓣，這樣使卷積後的功率譜不同於真正的功率譜。sinc函數的主瓣不是無限窄的，引起功率譜向附近頻域擴展，造成譜的模糊，降低譜的解析度；同時，由於sinc函數的旁瓣存在，導致能量向旁瓣中「泄漏」（稱之為旁瓣泄漏），即引起頻譜間的干擾，信號強的功率譜旁瓣影響信號弱的功率譜檢測，嚴重時，會使主瓣產生很大失真，檢測不出弱信號，或者把旁瓣誤認為是信號，造成假信號。為了對經典功率譜估計進行改進，可以採用各種不同的窗函數，但其結果都是以增加主瓣寬度來換取旁瓣的壓低，因此功率譜解析度低是經典功率譜估計的致命缺點。

為了克服以上缺點，人們曾做過長期努力，提出了平均、加窗平滑等辦法，在一定程度上改善了經典功率譜估計的性能。實踐證明，對於長數據記錄來說，以傅氏變換為基礎的經典功率譜估計方法，的確是比較實用的。但是，經典方法始終無法根本解決頻率解析度和功率譜估計穩定性之間的矛盾，特別是在數據記錄很短的情況下，這一矛盾顯得尤為突出。這就促進了現代功率譜估計方法研究的展開。

現代功率譜估計方法主要是以隨機過程（Stochastic Process）的參數模型（Parameter Model）為基礎的，稱之為參數模型方法。雖然說現代功率譜估計技術的研究和應用主要起始於60年代，但實際上，時間序列模型在非工程領域早已被採用，如Yule在1927年、Walker在1931年都曾使用過自回歸模型預測描述經濟的時間序列的發展趨勢，而Prony則早在1795年就曾採用指數模型去擬合在氣體化學實驗中獲得的數據。在統計學和數值分析領域中，人們也曾採用過模型方法。

現代功率譜估計的提出主要是針對經典功率譜估計（周期圖和自相關法）的解析度和方差性能不好的問題而提出的。1967年Burg在地震學研究中受到線性預測濾波的啟發，提出了最大熵譜估計方法，在提高解析度方面作了最有意義的探索。1968年Parzen正式提出了自回歸譜估計方法。1971年Van der Bos證明了一維最大熵譜估計與自回歸譜估計等效。1972年出現的譜估計的Prony方法在數學上與自回歸方法有某些類似。目前以自回歸滑動平均模型為基礎的譜估計已經比自回歸模型譜估計具有更高的頻率解析度和更好的性能。1973年Pisarenko提出的諧波分解方法提供了可靠的頻率估計方法。1981年Schmidt提出了譜估計的多信號分類（MUSIC）演算法等。因此，現代功率譜分析主要有ARMA譜分析、最大似然、熵譜估計和特徵分解四種方法。ARMA譜分析是一種建模方法，即通過平穩線性信號過程建立模型來估計功率譜密度；熵譜估計包括最大熵譜和最小交叉法；特徵分解也叫特徵構造法和子空間法，包括Pisarenko諧波分解法、Prony法、MUSIC法和ESPRIT法（用旋轉不變技術估計參數方法）。

現代功率譜估計研究仍側重於一維功率譜分析，而且大部分是建立在二階矩（相關函數、方差、功率譜密度）基礎上的。但由於功率譜密度是頻率的實函數，缺少相位信息，因此，建立在高階譜基礎上的譜估計方法正引起人們的注意，特別是雙譜估計和三譜估計的研究受到了高度的重視。其它如多維譜估計、多通道譜估計等的研究也正在發展中。人們希望這些新方法能更多地在提取信息、估計相位和描述非線性等方面獲得應用。

⑥ 高分求一段漢譯英，用軟體直接翻譯的勿擾。

With the development of science and technology, now in many areas of need to know the exact direction of useful signal, the frequency of parameters, the precise measurement of the signal to put a very high demand. Spectrum space technology in many areas has been widely used, and array signal processing plays a very important role.
Often, the signal space is measured up to one direction Wei Boda (DOA) estimates. In practical applications, often at the same time to signal the number of parameters of the estimates, such as signal DOA, pitch angle and frequency. This paper first introced for the two-dimensional space spectrum of the estimated MUSIC algorithm, ESPRIT algorithm for the basic idea and principle of algorithm, and then focus on the two algorithms based on the position - two-dimensional frequency spectrum estimation. Then this paper were used MATLAB language MUSIC algorithm and ESRPIT algorithm estimated two-dimensional and computer simulation results are analyzed and compared.

This paper first introced the concept of radio frequency identification technology, the rise and development. Radio frequency identification system with the basic structure and principles, described the antenna transceiver and the working principle of the process. Subsequent analysis of the transmitting power, receive power chip, chip power consumption and receiving antennas to receive the echo of the relationship between power and to examine in detail the chip to receive power and receiving antennas to receive the echo power as a range of changes, as well as the decision The size of these power factors, the actual design of radio frequency identification hardware to provide reference data

參考

⑦ 在陣列信號處理中，為什麼雜訊子空間與導向矢量是相互正交的

陣列信號處理經過多年來的發展，時域、空域上所蘊含的信息已得到充分的挖掘。近年來，新的研究工作在極化上得到開展。作為一種既能感知極化信息又能感知時空域信息的接收裝置，電磁矢量感測器陣列接收到的數據具有高維結構。然而，傳統的方法處理高維數據是將其轉化為矩陣數據，這使得數據的高維信息未得到有效利用。

利用張量這一處理高維數據的專用工具，研究了非完全極化波信號的BTD（Block
term
decompositions）模型、基於Tucker分解的參數估計方法、基於CP（Canoical
Decomposition）分解的參數估計方法和基於BTD分解的DOA估計方法。本文圍繞如...
展開
陣列信號處理經過多年來的發展，時域、空域上所蘊含的信息已得到充分的挖掘。近年來，新的研究工作在極化上得到開展。作為一種既能感知極化信息又能感知時空域信息的接收裝置，電磁矢量感測器陣列接收到的數據具有高維結構。然而，傳統的方法處理高維數據是將其轉化為矩陣數據，這使得數據的高維信息未得到有效利用。
本文利用張量這一處理高維數據的專用工具，研究了非完全極化波信號的BTD（Block
term
decompositions）模型、基於Tucker分解的參數估計方法、基於CP（Canoical
Decomposition）分解的參數估計方法和基於BTD分解的DOA估計方法。本文圍繞如何充分利用電磁矢量感測器陣列中蘊含的高維信息，
展開了以下幾個部分工作：

1、改進了基於張量的Tucker分解的MUSIC演算法。通過張量分解改善了估計的雜訊子空間與真實信號導向矢量的之間的正交性，提升了對電磁矢量感測器陣列參數估計的分辨力。相對於矩陣方法對雜訊子空間的估計，張量法利用到了數據的各個維度的信息，能從被雜訊覆蓋的數據中還原出更接近真實值的數據。模擬表明，利用張量分解估計的雜訊子空間進行MUSIC參數估計能夠提升來波信號的分辨力。這部分說明了張量各維度之間的整體性，即張量數據的高維特徵，是一種值得利用的信息。

2、改進了基於張量的CP分解對陣列進行盲估計的方法。相比於Tucker分解的估計方法，CP分解不僅能反映數據各維度之間的整體性，更重要的是還可以盲估計出各維度上的組成成分。利用這種盲估計特性，可以進一步將對應組成成分具有的結構特點作為約束去影響CP分解的每一步迭代。即這種盲估計特性為充分利用陣列自身具有的結構特點創造了條件。另外，本文提出了基於參數化的加入結構約束的方法。即在每一步迭代中，利用盲估計出來的組成成分估算參數，並根據參數和信號模型生成滿足模型特徵的組成成分。模擬表明，加入了結構信息的約束之後，參數估計誤差有所減小。
3、研究了基於
BTD分解的非完全極化波參數估計方法。對於非完全極化波的電磁矢量感測器陣列接收模型，同一方向的來波總是可以等效為兩路不相乾的信號被一個6×2矩陣導向。這時候，CP分解唯一性所需的條件不滿足，本文首次採用了符合這種陣列模型的BTD分解。BTD分解可以盲估計出空域陣列導向矢量和極化部分的一個列空間。另外，根據電磁場中的坡印廷定理，本文研究了從這個列空間中求解方向向量的方法。最後，融合空域陣列導向矢量和方向向量的估計得到來波方向。最後與ESPRIT演算法進行了模擬對比。結果表明，對於非完全極化波的接收模型，基於BTD分解的方法能夠得到更精確的估計結果。

⑧ DOA估計演算法

學號：20000300055

姓名：王鐸澎

嵌牛導讀：文章對DOA演算法進行了簡單的介紹。

嵌牛正文：https://blog.csdn.net/zhangziju/article/details/100730081?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522160689878119725222413438%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=160689878119725222413438&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~_landing_v2~default-1-100730081.pc_first_rank_v2_rank_v28&utm_term=Musicsuanfa&spm=1018.2118.3001.4449

DOA估計演算法

DOA（Direction Of Arrival）波達方向定位技術主要有ARMA譜分析、最大似然法、熵譜分析法和特徵分解法，特徵分解法主要有MUSIC演算法、ESPRIT演算法WSF演算法等。

MUSIC (Multiple Signal Classification)演算法，即多信號分類演算法，由Schmidt等人於1979年提出。MUSIC演算法是一種基於子空間分解的演算法，它利用信號子空間和雜訊子空間的正交性，構建空間譜函數，通過譜峰搜索，估計信號的參數。對於聲源定位來說，需要估計信號的DOA。MUSIC演算法對DOA的估計有很高的解析度，且對麥克風陣列的形狀沒有特殊要求，因此應用十分廣泛。

運用矩陣的定義，可得到更為簡潔的表達式：

X = A S + N X=AS+NX=AS+N

式中

X = [ x 1 ( t ) , x 2 ( t ) , . . . x M ( t ) ] T X=[x_1(t),x_2(t),...x_M(t)]^TX=[x1​(t),x2​(t),...xM​(t)]T，

S = [ S 1 ( t ) , S 2 ( t ) , . . . S D ( t ) ] T S=[S_1(t),S_2(t),...S_D(t)]^TS=[S1​(t),S2​(t),...SD​(t)]T，

A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T A=[a(\theta_1),a(\theta_2),...a(\theta_D)]^TA=[a(θ1​),a(θ2​),...a(θD​)]T，

N = [ n 1 ( t ) , n 2 ( t ) , . . . n M ( t ) ] T N=[n_1(t),n_2(t),...n_M(t)]^TN=[n1​(t),n2​(t),...nM​(t)]T。

X XX為陣元的輸出，A AA為方向響應向量，S SS是入射信號，N NN表示陣列雜訊。

其中 φ k = 2 π d λ s i n θ k \varphi_k=\frac{2\pi d}{\lambda}sin\theta_kφk​=λ2πd​sinθk​有

A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T = [ 1 1 ⋯ 1 e − j φ 1 e − j φ 2 ⋯ e − j φ D ⋮ ⋮ ⋱ ⋮ e − j ( M − 1 ) φ 1 e − j ( M − 1 ) φ 2 ⋯ e − j ( M − 1 ) φ D ] A=[a(\theta_1),a(\theta_2),...a(\theta_D)]^T=\left[

1e−jφ1⋮e−j(M−1)φ11e−jφ2⋮e−j(M−1)φ2⋯⋯⋱⋯1e−jφD⋮e−j(M−1)φD11⋯1e−jφ1e−jφ2⋯e−jφD⋮⋮⋱⋮e−j(M−1)φ1e−j(M−1)φ2⋯e−j(M−1)φD

\right]A=[a(θ1​),a(θ2​),...a(θD​)]T=⎣⎢⎢⎢⎡​1e−jφ1​⋮e−j(M−1)φ1​​1e−jφ2​⋮e−j(M−1)φ2​​⋯⋯⋱⋯​1e−jφD​⋮e−j(M−1)φD​​⎦⎥⎥⎥⎤​

對x m ( t ) x_m(t)xm​(t)進行N點采樣，要處理的問題就變成了通過輸出信號x m ( t ) x_m(t)xm​(t)的采樣{ x m ( i ) = 1 , 2 , . . . , M } \{ x_m (i)=1,2,...,M\}{xm​(i)=1,2,...,M}估計信號源的波達方向角θ 1 , θ 2 . . . θ D \theta_1,\theta_2...\theta_Dθ1​,θ2​...θD​,由此可以很自然的將陣列信號看作是雜訊干擾的若干空間諧波的疊加，從而將波達方向估計問題與譜估計聯系起來。

對陣列輸出X做相關處理，得到其協方差矩陣

R x = E [ X X H ] R_x=E[XX^H]Rx​=E[XXH]

其中H HH表示矩陣的共軛轉置。

根據已假設信號與雜訊互不相關、雜訊為零均值白雜訊，因此可得到：

R x = E [ ( A S + N ) ( A S + N ) H ] = A E [ S S H ] A H + E [ N N H ] = A R S A H + R N R_x=E[(AS+N)(AS+N)^H] =AE[SS^H]A^H+E[NN^H]=AR_SA^H+R_NRx​=E[(AS+N)(AS+N)H]=AE[SSH]AH+E[NNH]=ARS​AH+RN​

其中R s = E [ S S H ] R_s=E[SS^H]Rs​=E[SSH]稱為信號相關矩陣

R N = σ 2 I R_N=\sigma^2IRN​=σ2I是雜訊相關陣

σ 2 \sigma^2σ2是雜訊功率

I II是M × M M\times MM×M階的單位矩陣

在實際應用中通常無法直接得到R x R_xRx​，能使用的只有樣本的協方差矩陣：

R x ^ = 1 N ∑ i = 1 N X ( i ) X H ( i ) \hat{R_x}=\frac{1}{N} \sum_{i=1}^{N}X(i)X^H (i)Rx​^​=N1​∑i=1N​X(i)XH(i)，R x ^ \hat{R_x}Rx​^​是R x R_xRx​的最大似然估計。

當采樣數N → ∞ N\to\inftyN→∞，他們是一致的，但實際情況將由於樣本數有限而造成誤差。根據矩陣特徵分解的理論，可對陣列協方差矩陣進行特徵分解，首先考慮理想情況，即無雜訊的情況：R x = A R s A H R_x=AR_sA^HRx​=ARs​AH，對均勻線陣，矩陣A由

A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T = [ 1 1 ⋯ 1 e − j φ 1 e − j φ 2 ⋯ e − j φ D ⋮ ⋮ ⋱ ⋮ e − j ( M − 1 ) φ 1 e − j ( M − 1 ) φ 2 ⋯ e − j ( M − 1 ) φ D ] A=[a(\theta_1),a(\theta_2),...a(\theta_D)]^T=\left[

1e−jφ1⋮e−j(M−1)φ11e−jφ2⋮e−j(M−1)φ2⋯⋯⋱⋯1e−jφD⋮e−j(M−1)φD11⋯1e−jφ1e−jφ2⋯e−jφD⋮⋮⋱⋮e−j(M−1)φ1e−j(M−1)φ2⋯e−j(M−1)φD

\right]A=[a(θ1​),a(θ2​),...a(θD​)]T=⎣⎢⎢⎢⎡​1e−jφ1​⋮e−j(M−1)φ1​​1e−jφ2​⋮e−j(M−1)φ2​​⋯⋯⋱⋯​1e−jφD​⋮e−j(M−1)φD​​⎦⎥⎥⎥⎤​

所定義的范德蒙德矩陣，只要滿足θ i ≠ θ j , i ≠ j \theta_i\neq \theta_j,i\neq jθi​​=θj​,i​=j，則他的各列相互獨立。

若R s R_sRs​為非奇異矩陣R a n k ( R s ) = D Rank(R_s)=DRank(Rs​)=D，各信號源兩兩不相干，且M > D M>DM>D，則r a n d ( A R s A H ) = D rand(AR_sA^H)=Drand(ARs​AH)=D，

由於R x = E [ X X H ] R_x=E[XX^H]Rx​=E[XXH]，有：

R s H = R x R_s^H=R_xRsH​=Rx​

即R s R_sRs​為Hermite矩陣，它的特性是都是實數，又由於R s R_sRs​為正定的，因此A R s A … … H AR_sA……HARs​A……H為半正定的，它有D個正特徵值和M − D M-DM−D個零特徵值。

再考慮有雜訊存在的情況

R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2IRx​=ARs​AH+σ2I

由於σ 2 > 0 \sigma^2>0σ2>0，R x R_xRx​為滿秩陣，所以R x R_xRx​有M個正實特徵值λ 1 , λ 2 . . . λ M \lambda_1,\lambda_2...\lambda_Mλ1​,λ2​...λM​

分別對應於M個特徵向量v 1 , v 2 . . . v M v_1,v_2...v_Mv1​,v2​...vM​。又由於R x R_xRx​為Hermite矩陣，所以各特徵向量是正交的，即：v i H v j = 0 , i ≠ j v_i^Hv_j=0,i\neq jviH​vj​=0,i​=j與信號有關的特徵值只有D個,分別等於矩陣A R s A H AR_sA^HARs​AH的各特徵值與σ 2 \sigma^2σ2之和，其餘M − D M-DM−D個特徵值為σ 2 \sigma^2σ2，即σ 2 \sigma^2σ2為R RR的最小特徵值，它是M − D M-DM−D維的，對應的特徵向量v i , i = 1 , 2 , . . . , M v_i,i=1,2,...,Mvi​,i=1,2,...,M中，也有D個是與信號有關的，另外M − D M-DM−D個是與雜訊有關的，可利用特徵分解的性質求出信號源的波達方向θ k \theta_kθk​。

MUSIC演算法的原理及實現

通過對協方差矩陣的特徵值分解，可得到如下結論：

將矩陣R x R_xRx​的特徵值進行從小到大的排序，即λ 1 ≥ λ 2 ≥ . . . ≥ λ M > 0 \lambda_1 \geq \lambda_2\geq...\geq\lambda_M>0λ1​≥λ2​≥...≥λM​>0，其中D個較大的特徵值對應於信號，M − D M-DM−D個較小的特徵值對應於雜訊。

矩陣R x R_xRx​的屬於這些特徵值的特徵向量也分別對應於各個信號和雜訊，因此可把R x R_xRx​的特徵值（特徵向量）劃分為信號特徵（特徵向量）與雜訊特徵（特徵向量）。

設λ i \lambda_iλi​為R x R_xRx​的第i ii個特徵值，v i v_ivi​是與λ i \lambda_iλi​個相對應的特徵向量，有：

R x v i = λ i v i R_xv_i=\lambda_iv_iRx​vi​=λi​vi​

再設λ i = σ 2 \lambda_i=\sigma^2λi​=σ2是R x R_xRx​的最小特徵值R x v i = σ 2 v i i = D + 1 , D + 2... M R_xv_i=\sigma^2v_i i=D+1,D+2...MRx​vi​=σ2vi​i=D+1,D+2...M，

將R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2IRx​=ARs​AH+σ2I代入可得σ 2 v i = ( A R s A H + σ 2 I ) v i \sigma^2v_i=(AR_sA^H+\sigma^2I)v_iσ2vi​=(ARs​AH+σ2I)vi​，

將其右邊展開與左邊比較得：

A R s A H v i = 0 AR_sA^Hv_i=0ARs​AHvi​=0

因A H A A^HAAHA是D ∗ D D*DD∗D維的滿秩矩陣，( A H A ) − 1 (A^HA)^{-1}(AHA)−1存在；

而R s − 1 R_s^{-1}Rs−1​同樣存在，則上式兩邊同乘以R s − 1 ( A H A ) − 1 A H R_s^{-1}(A^HA)^{-1}A^HRs−1​(AHA)−1AH，

有：

R s − 1 ( A H A ) − 1 A H A R s A H v i = 0 R_s^{-1}(A^HA)^{-1}A^HAR_sA^Hv_i=0Rs−1​(AHA)−1AHARs​AHvi​=0

於是有

A H v i = 0 ， i = D + 1 , D + 2 , . . . , M A^Hv_i=0，i=D+1,D+2,...,MAHvi​=0，i=D+1,D+2,...,M

上式表明：雜訊特徵值所對應的特徵向量（稱為雜訊特徵向量）v i v_ivi​，與矩陣A AA的列向量正交，而A AA的各列是與信號源的方向相對應的，這就是利用雜訊特徵向量求解信號源方向的出發點。

用各雜訊特徵向量為例，構造一個雜訊矩陣E n E_nEn​:

E n = [ v D + 1 , v D + 2 , . . . v M ] E_n=[v_{D+1},v_{D+2},...v_{M}]En​=[vD+1​,vD+2​,...vM​]

定義空間譜P m u ( θ ) P_{mu}(\theta)Pmu​(θ):

P m u ( θ ) = 1 a H ( θ ) E n E n H a ( θ ) = 1 ∥ E n H a ( θ ) ∥ 2 P_{mu}(\theta)=\frac{1}{{a^H}(\theta)}E_nE_n^Ha(\theta)=\frac{1}{\Vert E_n^Ha(\theta)\Vert^2}Pmu​(θ)=aH(θ)1​En​EnH​a(θ)=∥EnH​a(θ)∥21​

該式中分母是信號向量和雜訊矩陣的內積，當a ( θ ) a(\theta)a(θ)和E n E_nEn​的各列正交時，該分母為零，但由於雜訊的存在，它實際上為一最小值，因此P m u ( θ ) P_{mu}(\theta)Pmu​(θ)有一尖峰值，由該式，使θ \thetaθ變化，通過尋找波峰來估計到達角。

MUSIC演算法實現的步驟

1.根據N個接收信號矢量得到下面協方差矩陣的估計值：

R x = 1 N ∑ i = 1 N X ( i ) X H ( i ) R_x=\frac{1}{N} \sum_{i=1}^NX(i)X^H(i)Rx​=N1​∑i=1N​X(i)XH(i)

對上面得到的協方差矩陣進行特徵分解

R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2IRx​=ARs​AH+σ2I

2.按特徵值的大小排序 將與信號個數D DD相等的特徵值和對應的特徵向量看做信號部分空間，將剩下的M − D M-DM−D個特徵值和特徵向量看做雜訊部分空間，得到雜訊矩陣E n E_nEn​:

A H v i = 0 ， i = D + 1 , D + 2 , . . . , M A^Hv_i=0，i=D+1,D+2,...,MAHvi​=0，i=D+1,D+2,...,M

E n = [ v D + 1 , v D + 2 , . . . v M ] E_n=[v_{D+1},v_{D+2},...v_{M}]En​=[vD+1​,vD+2​,...vM​]

3.使θ \thetaθ變化 ，按照式

P m u ( θ ) = 1 a H ( θ ) E n E n H a ( θ ) P_{mu}(\theta)=\frac{1}{{a^H}(\theta)E_nE_n^Ha(\theta)}Pmu​(θ)=aH(θ)En​EnH​a(θ)1​

來計算譜函數，通過尋求峰值來得到波達方向的估計值。

clear; close all;

%%%%%%%% MUSIC for Uniform Linear Array%%%%%%%%

derad = pi/180;      %角度->弧度

N = 8;              % 陣元個數       

M = 3;              % 信源數目

theta = [-30 0 60];  % 待估計角度

snr = 10;            % 信噪比

K = 512;            % 快拍數

dd = 0.5;            % 陣元間距

d=0:dd:(N-1)*dd;

A=exp(-1i*2*pi*d.'*sin(theta*derad));  %方向矢量

%%%%構建信號模型%%%%%

S=randn(M,K);            %信源信號，入射信號

X=A*S;                    %構造接收信號

X1=awgn(X,snr,'measured'); %將白色高斯雜訊添加到信號中

% 計算協方差矩陣

Rxx=X1*X1'/K;

% 特徵值分解

[EV,D]=eig(Rxx);                  %特徵值分解

EVA=diag(D)';                      %將特徵值矩陣對角線提取並轉為一行

[EVA,I]=sort(EVA);                %將特徵值排序 從小到大

EV=fliplr(EV(:,I));                % 對應特徵矢量排序

% 遍歷每個角度，計算空間譜

for iang = 1:361

    angle(iang)=(iang-181)/2;

    phim=derad*angle(iang);

    a=exp(-1i*2*pi*d*sin(phim)).';

    En=EV(:,M+1:N);                  % 取矩陣的第M+1到N列組成雜訊子空間

    Pmusic(iang)=1/(a'*En*En'*a);

end

Pmusic=abs(Pmusic);

Pmmax=max(Pmusic)

Pmusic=10*log10(Pmusic/Pmmax);            % 歸一化處理

h=plot(angle,Pmusic);

set(h,'Linewidth',2);

xlabel('入射角/(degree)');

ylabel('空間譜/(dB)');

set(gca, 'XTick',[-90:30:90]);

grid on;

實現結果