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

背景技術(shù)：

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

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

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

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

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

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

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

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

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

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

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

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

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

s7、令a′ij＝aij，進(jìn)入s8；

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

s9、計(jì)算a＝q^haq和v＝q^hv，其中，q∈c^n×n為復(fù)數(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是否成立，是則進(jìn)入s11，否則j＝j(luò)+1后跳轉(zhuǎn)到s4；

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

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

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

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

進(jìn)一步地，s1所述a＝e×10^-1。

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

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

附圖說明

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

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

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

具體實(shí)施方式

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

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

實(shí)施例1、

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

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

在估計(jì)性能包括計(jì)算精度、計(jì)算速度和資源消耗，具體用下面指標(biāo)評價(jià)：

1.計(jì)算精度：

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

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

2.計(jì)算速度：

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

3.資源消耗：

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

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

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

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

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

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

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

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

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

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

b1.進(jìn)行初始化，具體方法為：

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

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

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

b13.初始化清掃元aij的行列下標(biāo)，i＝1、j＝2，進(jìn)入步驟b2。

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

b21.在矩陣a中選取aij，進(jìn)入步驟b22。

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

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

b24.進(jìn)行位擴(kuò)展，即a′ij＝aij×2^m，進(jìn)入步驟b26。

b25.不進(jìn)行位擴(kuò)展，即a′ij＝aij，進(jìn)入步驟b26。

b26.計(jì)算相角和模值，即和進(jìn)入步驟b27。

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

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

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

b3.對迭代過程進(jìn)行判斷。

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

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

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

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

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

c1.由下式構(gòu)造酉變換矩陣t∈c^n×n，進(jìn)入步驟c2。

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

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

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

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

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

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

仿真結(jié)果為：

計(jì)算精度：

計(jì)算速度:nclk＝11710

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

此時(shí)，估計(jì)信號(hào)所對應(yīng)的圓盤半徑計(jì)算精度為ε1≈10^-9；估計(jì)噪聲所對應(yīng)的圓盤半徑計(jì)算精度為εi≈10^-4，i＝2,3,…,7；圓盤平均計(jì)算精度

實(shí)施例2、

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

實(shí)施例2的方法如附圖2所示，其余仿真條件與實(shí)施例1的相同，進(jìn)行信號(hào)與噪聲對應(yīng)的蓋氏圓盤半徑的估計(jì)。

實(shí)施例2的評價(jià)標(biāo)準(zhǔn)與實(shí)施例1的一致。

仿真結(jié)果為：

計(jì)算精度：

計(jì)算速度:n′clk＝17960

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

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

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

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