本發(fā)明屬于信號處理領域，尤其涉及種基于jacobi迭代算法的高精度矩陣特征值分解實現(xiàn)方法。

背景技術：

在信號處理中，矩陣的特征值分解evd是一個應用廣泛的矩陣運算。如數(shù)據(jù)壓縮、噪聲去除、數(shù)值分析，包括近幾年興起的機器學習、深度學習其基本核心操作也包括矩陣特征值分解。實現(xiàn)矩陣特征值分解的常用方法有gauss變換、householder變換、jacobi迭代等，其中，jacobi迭代是精度較高的方法，并且很適合在fpga中實現(xiàn)。因此一種基于jacobi迭代算法的高精度矩陣特征值分解實現(xiàn)技術在實際工程中具有很高的應用價值。

經(jīng)典的jacobi迭代算法計算共軛矩陣a∈c^n×n的特征值分解如圖1所示，這種經(jīng)典的迭代算法雖然有較快收斂速度，但是該算法需要在矩陣a的眾多元素中選取aij，使得aij為非對角元素中絕對值最大的一個，再進行后面的計算操作。這樣每一步都要尋找絕對值最大的非對角元，比較費時也不適合在fpga實現(xiàn)，因此經(jīng)典的jacobi迭代算法在實際工程中并不實用。

目前實際工程中多數(shù)采用如圖2所示的循環(huán)jacobi迭代算法，通過逐行掃描遍歷法選取aij，這樣避免了尋找最大絕對值的非對角元的復雜繁瑣步驟。這樣選取aij的方式，在aij數(shù)值比較大時，fpga中使用cordic算法計算φ、θ誤差比較小，可以取得比較好的效果。但當aij比較小甚至接近0時，此時fpga中使用cordic算法計算φ、θ誤差比較大，就會導致后面計算a＝q^haq產(chǎn)生誤差，其中q∈c^n×n為復數(shù)域內(nèi)的平面旋轉(zhuǎn)矩陣。而計算過程需要多次的迭代運算，如果在迭代過程中出現(xiàn)多次aij比較小甚至接近0的情況，就會產(chǎn)生較大的累計誤差，從而使最終計算結果的精度相對較差。

技術實現(xiàn)要素：

發(fā)明的目的在于解決在循環(huán)jacobi迭代算法進行矩陣特征值分解過程中，因為逐行掃描遍歷法選取aij在aij比較小甚至接近于0時，導致fpga中使用cordic算法計算φ、θ誤差比較大，進而使迭代過程產(chǎn)生相對較大的累計誤差，導致計算結果誤差增大。即提供一種基于jacobi迭代算法的高精度矩陣特征值分解實現(xiàn)方法，在不明顯增加算法復雜度、fpga實現(xiàn)難度與增加資源消耗的情況下，提高實際工程中基于循環(huán)jacobi迭代算法的fpga實現(xiàn)矩陣特征值分解的計算精度。

一種基于jacobi迭代算法的高精度矩陣特征值分解實現(xiàn)方法，包括如下步驟：

s1、設數(shù)據(jù)矩陣a∈c^n×n為共軛矩陣，并且設定最大遍歷次數(shù)為t、最小清掃門限a、擴位門限b和算術左移位數(shù)m，其中，最小清掃門限a應小于要求計算結果的精度一個數(shù)量級，擴位門限b與算術左移位數(shù)m與fpga實現(xiàn)里數(shù)據(jù)位寬size有關，滿足b×2^m＜2^size-4保證計算結果不溢出，n為不為零的自然數(shù)，共軛矩陣a中的元素為aij，i＝1,2,3,...,n，j＝1,2,3,...,n，1≤t≤t且t為自然數(shù)；

s2、初始化遍歷次數(shù)計數(shù)器，令t＝0，

初始化特征向量初始矩陣，令v＝e，其中，e為單位陣；

s3、在s1所述共軛矩陣a中選取aij，初始化清掃元aij行列下標，令i＝1，j＝2；

s4、判斷aij是否滿足跳過清掃條件|real(aij)|＜a&|imag(aij)|＜a，若滿足則轉(zhuǎn)入s10，如不滿足則轉(zhuǎn)入s4；

s5、判斷aij是否滿足擴展條件|real(aij)|＜b&|imag(aij)|＜b，若滿足則轉(zhuǎn)入s6，若不滿足則轉(zhuǎn)入s7；

s6、進行位擴展，即計算a′ij＝aij×2^m，轉(zhuǎn)入s7；

s7、令a′ij＝aij，進入s8；

s8、計算根據(jù)所得計算

s9、計算a＝q^haq和v＝q^hv，其中，q∈c^n×n為復數(shù)域內(nèi)的平面旋轉(zhuǎn)矩陣，

即q的對角元素中除了qii＝e^jφcosθ、qjj＝e^-jφcosθ之外其他均為1，非對角元素中除了qij＝-e^jφsinθ、qji＝e^-jφsinθ其他元素均為0，θ為旋轉(zhuǎn)角度；

s10、判斷j＝n是否成立，是則進入s11，否則j＝j+1后跳轉(zhuǎn)到s4；

s11、判斷i＝n-1是否成立，是則進入s12，否則i＝i+1，j＝i+1后跳轉(zhuǎn)到s4；

s12、判斷t＝t是否成立，是則進入s13，否則t＝t+1后跳轉(zhuǎn)到s3；

s13、輸出迭代計算結果a和v，其中，a的對角元數(shù)位s1輸入數(shù)據(jù)矩陣a的特征值，v為對應的特征向量矩陣。

進一步地，s1所述遍歷次數(shù)t越大迭代次數(shù)越多，計算越精確，但計算時間越長，為了取得速度與精度的平衡，當n≤8時t＝3，當n＞8時t＝6。

進一步地，s1所述a＝e×10^-1。

本發(fā)明的有益效果是：

在不明顯增加算法復雜度、fpga實現(xiàn)難度與硬件資源消耗、計算消耗時間的情況下，提高基于循環(huán)jacobi迭代算法的fpga實現(xiàn)矩陣特征值分解的計算精度和計算速度，在實際工程中具有重要價值。

附圖說明

圖1為經(jīng)典jacobi迭代算法流程。

圖2為循環(huán)jacobi迭代算法流程。

圖3為本發(fā)明算法流程。

具體實施方式

下面將結合實施例和附圖，對本發(fā)明方法進行進一步說明。

本發(fā)明應用于估計信號與噪聲對應的蓋氏圓盤半徑，提高圓盤半徑計算精度，和計算速度。

實施例1、

接收陣列為8陣元組成的均勻線陣。

如圖3所示，考慮n＝1個載波頻率為的bpsk調(diào)制的遠場信號s(k)，以γ＝0°的方向入射到陣元數(shù)n＝8的均勻線陣上，且陣元間距為d＝0.5λ，其中，λ為信號波長，陣列接收噪聲是功率為σ²＝1的零均值高斯白噪聲，接收信號信噪比snr＝20db，快拍數(shù)為l＝1024。通過進行特征值分解估計信號與噪聲對應的蓋氏圓盤半徑i＝1,2,…,n-1。

在估計性能包括計算精度、計算速度和資源消耗，具體用下面指標評價：

1.計算精度：

(1).圓盤半徑計算精度：i＝1,2,…,n-1其中κi為圓盤半徑的理論值。εi越小表示計算精度越高。

(2).圓盤平均計算精度：越小表示平均計算精度越高。

2.計算速度：

(1).計算消耗的時鐘數(shù)nclk，越小表示計算消耗時間越少，計算速度越快。

3.資源消耗：

(1).寄存器消耗數(shù)量nreg，越小對應寄存器資源消耗越少。

(2).邏輯門消耗數(shù)量nlut，越小對應邏輯門資源消耗越少。

應用特征值分解估計信號與噪聲對應的蓋氏圓盤半徑包括，a.仿真接收信號數(shù)據(jù)建模、b.應用本發(fā)明進行特征值分解、c.計算圓盤半徑，具體為以下步驟：

a.仿真接收信號數(shù)據(jù)建模。

a1.由下式產(chǎn)生陣元數(shù)n＝8的陣列接收信號向量x(k)＝[x1(k)x2(k)…x8(k)]^h進入步驟a2。

x(k)＝a(γ)s(k)+n(k)，k＝1,2,…,l

式中，n(k)為8×1為均值為零、方差σ²＝1的高斯白噪聲矢量；遠場接收信號s(k)＝ass(k)，其中其幅度as＝10^snr/20；a(γ)＝[1e^-jφ…e^-j(n-1)φ]^t，為空間陣列的n×1維流型矩陣。

a2.由計算數(shù)據(jù)協(xié)方差矩陣r∈c^n×n，進入步驟a3。

a3.根據(jù)對數(shù)據(jù)協(xié)方差矩陣進行分塊，得到塊矩陣r′∈c^(n-1)×(n-1)，進入步驟b，其中，r′∈c^(n-1)×(n-1)，r∈c^(n-1)×1，rnn為數(shù)據(jù)矩陣r第n行第n列的元素。

b.應用本發(fā)明對塊矩陣r′進行特征值分解，計算r′的特征值矩陣d與特征向量矩陣v。

b1.進行初始化，具體方法為：

b11.設數(shù)據(jù)矩陣a＝r′為共軛矩陣，并且設定遍歷次數(shù)t＝5、最小清掃門限a＝10^-8、擴位門限b＝10^-5和算術左移位數(shù)m＝8，進入步驟b12。

其中t越大迭代次數(shù)越多計算越精確但計算時間越長，根據(jù)矩陣維數(shù)n進行選取，當n≤8時t＝4，n＞8時t＝6可以取得速度與精度的平衡；最小清掃門限a與計算精度要求e有關，a＝e×10^-1，比如要求計算精度為10^-5則a≈10^-6。擴位門限b與算術左移位數(shù)m與fpga實現(xiàn)里數(shù)據(jù)位寬size有關，滿足b×2^m＜2^size-4保證計算結果不溢出。

b12.初始化遍歷次數(shù)計數(shù)器和特征向量初始矩陣，t＝0、v＝e，其中，e為單位陣，進入步驟b13。

b13.初始化清掃元aij的行列下標，i＝1、j＝2，進入步驟b2。

b2.進行jaocbi旋轉(zhuǎn)，具體方法如下：

b21.在矩陣a中選取aij，進入步驟b22。

b22.判斷是否滿足跳過清掃條件|real(aij)|＜a&|imag(aij)|＜a，是則跳轉(zhuǎn)到步驟b3，否則進入步驟b23。

b23.判斷是否滿足擴展條件|real(aij)|＜b&|imag(aij)|＜b，是則進入步驟b24，否則進入步驟b25。

b24.進行位擴展，即a′ij＝aij×2^m，進入步驟b26。

b25.不進行位擴展，即a′ij＝aij，進入步驟b26。

b26.計算相角和模值，即和進入步驟b27。

b27.計算平面旋轉(zhuǎn)角，即進入步驟b28。

b28.進行jacobi旋轉(zhuǎn)，即計算a＝q^haq和v＝q^hv，其中q∈c^n×n為復數(shù)域內(nèi)的平面旋轉(zhuǎn)矩陣。

即q的對角元素中除了qii＝e^jφcosθ、qjj＝e^-jφcosθ之外其他均為1，非對角元素中除了qij＝-e^jφsinθ、qji＝e^-jφsinθ其他元素均為0，進入步驟b3。

b3.對迭代過程進行判斷。

b31.判斷j＝n是否成立，是則進入步驟b32，否則j＝j+1后跳轉(zhuǎn)到步驟b21。

b32.判斷i＝n-1是否成立，是則進入步驟b33，否則i＝i+1，j＝i+1后跳轉(zhuǎn)到步驟b21。

b33.判斷t＝t是否成立，是則進入步驟b4，否則t＝t+1后跳轉(zhuǎn)到步驟b13。

b4.輸出迭代計算結果a和v，其中a的對角元就是數(shù)據(jù)矩陣a的特征值，v為對應的特征向量矩陣，進入步驟c。

c.對數(shù)據(jù)協(xié)方差矩陣r進行酉變換，計算信號與噪聲對應的圓盤半徑。

c1.由下式構造酉變換矩陣t∈c^n×n，進入步驟c2。

其中，v∈c^(n-1)×(n-1)為前面計算塊矩陣r′的特征向量，滿足vv^h＝e，e為單位陣。

c2.進行酉變換得到圓盤半徑，即計算下式，進入步驟c3。

式中，λi，i＝1,2,…,n-1為塊矩陣r′的特征值。

c3.由ri＝|ρi|，i＝1,2,…,n-1計算圓盤半徑ri，進入步驟c4。

c4.計算圓盤半徑計算精度：i＝1,2,…,n-1其中κi為圓盤半徑的理論值，和圓盤平均計算精度：進入步驟c5。

c5.統(tǒng)計計算消耗的時鐘數(shù)nclk、寄存器消耗數(shù)量nreg和邏輯門消耗數(shù)量nlut，算法結束。

仿真結果為：

計算精度：

計算速度:nclk＝11710

資源消耗:nreg＝29104、nlut＝30254

此時，估計信號所對應的圓盤半徑計算精度為ε1≈10^-9；估計噪聲所對應的圓盤半徑計算精度為εi≈10^-4，i＝2,3,…,7；圓盤平均計算精度

實施例2、

經(jīng)典方法循環(huán)jacobi算法應用于估計信號與噪聲對應的蓋氏圓盤半徑的估計性能，作為實施例1的對比例。

實施例2的方法如附圖2所示，其余仿真條件與實施例1的相同，進行信號與噪聲對應的蓋氏圓盤半徑的估計。

實施例2的評價標準與實施例1的一致。

仿真結果為：

計算精度：

計算速度:n′clk＝17960

資源消耗:n′reg＝29101、n′lut＝29998

本發(fā)明此時，估計信號所對應的圓盤半徑計算精度為ε1≈10^-8；估計噪聲所對應的圓盤半徑計算精度為εi≈10^-1，i＝2,3,…,7；圓盤平均計算精度

綜上所述，對比實施例1與實施例2，本發(fā)明相對于經(jīng)典方法在增加(nreg-n′reg)/n′reg×％≈0.01％寄存器資源消耗，(nlut-n′lut)/n′lut×％≈0.85％查找表資源消耗的情況下，平均計算精度從10^-1提高到了10^-4數(shù)量級，提高了3個數(shù)量級，同時計算速度提高了|nclk-n′clk|/n′clk×％≈34.8％。

所以，本發(fā)明在基本不增加資源消耗的情況下，不僅可以提高計算精度，還可以提高計算速度，在實際工程中具有重要價值。