本發(fā)明涉及數(shù)字處理技術領域，尤其涉及基于cordic算法的三角函數(shù)計算與矢量變換的并行計算。

背景技術：

坐標旋轉(zhuǎn)數(shù)字計算機cordic(coordinaterotationdigitalcomputer，簡稱cordic)由j.volder于1959年提出的，它是一種用于計算常用超越函數(shù)的循環(huán)迭代算法，其基本思想是是通過一系列的只與運算基數(shù)有關的固定小角度的不斷偏擺來逼近所要旋轉(zhuǎn)的角度。為了擴展可計算函數(shù)的個數(shù)，1971年j.wahher提出了統(tǒng)一cordic算法，即把圓周坐標、線性坐標和雙曲坐標統(tǒng)一到同一個cordic迭代方程中。

cordic算法因為能夠?qū)⒍喾N難以用硬件電路直接實現(xiàn)的復雜運算分解為簡單的加減法和移位操作，所以很適合用數(shù)字電路來實現(xiàn)，因而其應用范圍也就極其廣泛，比如數(shù)控振蕩器、正余弦函數(shù)發(fā)生器、數(shù)字頻率合成器等。

下面介紹旋轉(zhuǎn)模式圓周坐標下的cordic算法：

在平面坐標系下，將向量(x1，y1)旋轉(zhuǎn)到如附圖1的向量(x2，y2)，兩個向量之間滿足如式1所示的運算關系。

對所旋轉(zhuǎn)θ角度進行分解，將其分成n個遞減的小旋轉(zhuǎn)角θi，即其中θi，δi為判決算子，用于確定旋轉(zhuǎn)的方向。當順時針旋轉(zhuǎn)θi角度時，δi為-1；當逆時針旋轉(zhuǎn)θi角度時，δi為1.每次小的角度旋轉(zhuǎn)都有

再令θi＝tan^-12^-i，即tanθi＝2^-i。

這樣，平面坐標系中的向量(x1，y1)經(jīng)過n次圓周旋轉(zhuǎn)之后達到同一坐標平面中的向量(x2，y2)，該旋轉(zhuǎn)過程可表示為

其中k為伸縮因子

其中，k值可以提前計算出，k＝0.607529350088。

此時若令x1＝k，y1＝0，可得出x2＝cosθ，y2＝sinθ。

通過上述的理論推導，平面向量(x1，y1)的旋轉(zhuǎn)計算問題就可以由如下基本計算式迭代n次實現(xiàn)。此時第i次的迭代公式就轉(zhuǎn)變?yōu)橄率剑?/p>

由上式可知道，傳統(tǒng)的cordic算法中，每一步的旋轉(zhuǎn)方向都要根據(jù)上一步的結(jié)果進行判斷，即δi＝sign(zi)。要想得到n位的精度，需要至少迭代n次，需要n個時鐘周期才能完成一次計算。為了提高運算速度，最有效的辦法就是減少z通路的計算，提高運算的并行度。

tso-bingjuang，shen-fuhsiao等人于2004年提出并行cordic算法。其思想是根據(jù)輸入角度的n位二進制表示，將其分為高m位和低n-m位，先對高位實施binary-to-bipolar(bbr)轉(zhuǎn)換，一次產(chǎn)生高位的旋轉(zhuǎn)次序，然后將高位產(chǎn)生的角度誤差疊加在低位，然后在對低位采用bbr轉(zhuǎn)換，一次產(chǎn)生低位的旋轉(zhuǎn)次序。為了維持統(tǒng)一的矯正因子k，對于一次高位旋轉(zhuǎn)2^-i的弧度，采用若干次統(tǒng)一方向的微旋轉(zhuǎn)進行補償(此方法稱為microrotationanglerecoding)，使其在第i次旋轉(zhuǎn)產(chǎn)生的誤差可以在低位進行補償而不會產(chǎn)生高位的精度損失。該方法的雖然實現(xiàn)了cordic的并行處理，但增加了很多補償?shù)奈⑿D(zhuǎn)(在輸入角度位數(shù)為32位時，需要額外8次微旋轉(zhuǎn))，造成x/y路徑的延遲加大。公布號為cn102073471a，已授權(quán)生效的發(fā)明專利“一種處理器cordic迭代運算方法及電路”將三次迭代運算并行展開，目的是想一次計算三次迭代，但是這三次迭代的方向無法在開始迭代前確定，還是需要在計算通路中確定，實際上并沒有使運算加速多少。目前大部分cordic算法都是采用流水線計算或者迭代結(jié)構(gòu)。得到一個計算的值往往需要很多時鐘周期，非常不利于cordic算法在fpga上的快速實現(xiàn)。

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

針對現(xiàn)有技術的缺點，提出一種新型的低z通路判斷，低延遲，高精度的并行cordic計算方法和裝置。

旋轉(zhuǎn)弧度可用二進制表示為其中i∈{0,1}，θ≤π/4。

將輸入的弧度拆分成三個部分：高位，中位，低位：

m是滿足2^-m-tan^-12^-m<2^-n的最小值，可得出m＝[(n-log23)/3]。當i≥m時，tan2^-i＝2^-i。l為n/2，此時模矯正因子

根據(jù)現(xiàn)有技術對θ^h進行bbr轉(zhuǎn)換，可得到：

根據(jù)(9)式，因為bi-1∈{0,1}，所以ri∈{-1,1}，由此可根據(jù)輸入弧度的高m-1位直接得到1到m次旋轉(zhuǎn)方向及r1,r2,…,rm的值。

此時，理想的情況下1到m次的旋轉(zhuǎn)的角度都為2^-i，i＝1,2,…,m。但此時i＝1,2,…,m-1.硬件上難以實現(xiàn)tan2^-i。

對弧度2^-i進行量化處理可由下式表示：

為第i次旋轉(zhuǎn)的角度值。如圖2所示，實線為第i次旋轉(zhuǎn)目標角度，虛線為實際旋轉(zhuǎn)角度，若旋轉(zhuǎn)超過目標角度，角度誤差-εi，若旋轉(zhuǎn)未達到目標角度，角度誤差為εi。

此時第i次變換矩陣變?yōu)椋?/p>

σj的取值滿足式12的前提下，取

1到m-1次旋轉(zhuǎn)由式(11)替代，由式(10)可得θ^h產(chǎn)生的剩余角度為：

結(jié)合式(12)和式(16)可知所以σj滿足式(12)就可以將高位產(chǎn)生的誤差傳遞到低位而不會丟失精度。變換矩陣(11)硬件上可由簡單的移位和相加實現(xiàn)。

進一步的經(jīng)過前m次轉(zhuǎn)換，剩余弧度為:

對θ^m′繼續(xù)采用與θ^h同樣的轉(zhuǎn)換，得到:

其中rm＝1-2b′m-1,ri＝(2bi-1-1)，i＝m+1，m+2,…，l.

根據(jù)式(18)可得rm,rm+1,…,rl次旋轉(zhuǎn)的值，因為i在這個區(qū)間內(nèi)，滿足tan^-12^-i＝2^-i，轉(zhuǎn)換矩陣變?yōu)閭鹘y(tǒng)矩陣如下：

經(jīng)過第二階段l-m次旋轉(zhuǎn)，剩余弧度為:

在i∈[l,n]，由于模矯正因子cosθi＝1，可以得出若后續(xù)的某些旋轉(zhuǎn)沒有進行，也能保持統(tǒng)一的模校正因子，所以無需進行變換，此時可根據(jù)θ^l″進行簡單的單向旋轉(zhuǎn)，對于bi″＝0的不進行旋轉(zhuǎn)，bi″＝1的進行單向旋轉(zhuǎn)。直接根據(jù)二進制所在位的值，直接確定剩余旋轉(zhuǎn)的次序。更進一步的，當i>n/2，任意的兩次相鄰的旋轉(zhuǎn)由下式表示：

式(21)中2^-(m+n)<2^-n，已超出了最小精度，所以式(21)可簡化為：

由式(22)可得n/2到n次旋轉(zhuǎn)可壓縮成下面矩陣表示：

式(23)中，若剩余弧度為負，則s為-1，否則為1。式(23)可由兩次無符號乘法，兩次加法實現(xiàn)。

上述中旋轉(zhuǎn)方向的預測可分為兩步，第一步根據(jù)θ^h得出r1,r2,…,rm的值，第二步根據(jù)θ^m′+θ^l′，對θ^m′和θ^l′采用不同的預測方法，可同時得出rm,rm+1,…,rl,rl,…,rn的值。所以整體的并行cordic運算只要兩拍。

根據(jù)上述推導，偽旋轉(zhuǎn)產(chǎn)生的總校正因子可以提前計算出來，不會因為輸入角度不同而不同，保留了傳統(tǒng)cordic算法的優(yōu)點，且不需要額外的微旋轉(zhuǎn)進行補償。

附圖說明：

圖1是cordic算法旋轉(zhuǎn)示意圖。

圖2是第i次旋轉(zhuǎn)誤差產(chǎn)生示意圖。

圖3是式(11)的旋轉(zhuǎn)矩陣結(jié)構(gòu)示意圖。

圖4是本實施例32位系統(tǒng)第一旋轉(zhuǎn)矩陣硬件結(jié)構(gòu)示意圖。

圖5是本實施例總體結(jié)構(gòu)示意圖。

具體實施方式：

本發(fā)明的具體實施步驟為：

(1)對輸入的xin，yin，zin進行變換，使輸入角度值調(diào)整到(0，π/4)區(qū)間內(nèi)。對調(diào)整后的弧度拆分成三個部分：高位，中位，低位：

(2)對θ的前m-1次進行預測得出r1,r2,…,rm的值，對前m-1旋轉(zhuǎn)的目標角度進行量化，產(chǎn)生新的旋轉(zhuǎn)矩陣。將前m次旋轉(zhuǎn)誤差累計到低位，得到新的θ^m′+θ^l′，根據(jù)θ^m′+θ^l′，對θ^m′和θ^l′采用不同的預測方法，可同時得出rm,rm+1,…,rl,rl,…,rn的值。

(3)對m到l次旋轉(zhuǎn)采用傳統(tǒng)矩陣，對l到n次旋轉(zhuǎn)利用乘法直接求得。最后進行象限恢復輸出計算得到的三角函數(shù)值。

本實施例以32位的系統(tǒng)為例進行說明，對于一個32位系統(tǒng)，最大的弧度為2π，所以采用3位二進制表示整數(shù)和29位二進制表示小數(shù)。對于輸入弧度不在(0，π/4)區(qū)間中的，根據(jù)現(xiàn)有的區(qū)間折疊技術，很容易將輸入的弧度調(diào)整到(0，π/4)中，這里不在贅述。

調(diào)整后的弧度可表示為其中bi∈{0,1}，θ≤π/4。

將輸入的弧度拆分成三個部分：高位，中位，低位：

對于θ^h進行bbr轉(zhuǎn)換，得到：

根據(jù)式(27)可一次得到r1,r2,…,r10的值。

前10次的微旋轉(zhuǎn)角度和矩陣為：

θ1＝tan^-1(2^-1+2^-5+2^-6)(28)

θ2＝tan^-1(2^-2+2^-8+2^-10)(30)

θ3＝tan^-1(2^-3+2^-11+2^-13)(32)

θ4＝tan^-12^-4(34)

θ5＝tan^-12^-5(35)

θ6＝tan^-12^-6(37)

θ7＝tan^-12^-7(38)

θ8＝tan^-12^-8(39)

θ9＝tan^-12^-9(40)

θ10＝tan^-12^-10＝2^-10(41)

4到10次采用傳統(tǒng)的旋轉(zhuǎn)矩陣

由上述各個微旋轉(zhuǎn)產(chǎn)生的誤差為：

ε1＝|2^-1-θ1|＝0.00044＝o(2^-12)(42)

ε2＝|2^-2-θ2|＝0.00043＝o(2^-12)(43)

ε3＝|2^-3-θ3|＝0.00008＝o(2^-14)(44)

ε4＝|2^-4-θ4|＝0.000081＝o(2^-14)(45)

ε5＝|2^-5-θ5|＝0.00001＝o(2^-17)(46)

ε6＝|2^-6-θ6|＝0.000001＝o(2^-20)(47)

ε7＝|2^-7-θ7|＝o(2^-23)(48)

ε8＝|2^-8-θ8|<o(2^-23)(49)

ε9＝|2^-9-θ9|<o(2^-23)(50)

由θ^h產(chǎn)生的誤差最大弧度為

預先將每一個微旋轉(zhuǎn)的誤差值制表保存。每一步旋轉(zhuǎn)的誤差值用保留進位加法器疊加產(chǎn)生新的

對中位θ^m′進行數(shù)學轉(zhuǎn)換可得到

r10＝1-2b′9,ri＝(2bi-1-1)，i＝11,12,…,15(51)

可直接得到r10,r11,r12,r13,r14,r15的值。對i∈[10,15]采用傳統(tǒng)的旋轉(zhuǎn)矩陣。

后15到29次旋轉(zhuǎn)矩陣可由下式表示：

根據(jù)式的原碼表示，若為負，s＝-1，反之，s＝1。ri與位權(quán)2^-i所在位的值相等，ri∈{0,1}。15到29次旋轉(zhuǎn)直接由兩次乘法，兩次加法實現(xiàn)。上述中旋轉(zhuǎn)方向的預測可分為兩步，第一步根據(jù)θ^h得出r1,r2,…,r9的值，第二步根據(jù)θ^m′+θ^l′，對θ^m′和θ^l′采用不同的預測方法，可同時得出r10,…,r15,r15,…,r29的值。

下面結(jié)合附圖5進行說明：

附圖5中1為前處理模塊，采用區(qū)間折疊技術，對輸入的xin，yin，zin進行變換，使輸入角度值調(diào)整到(0，π/4)區(qū)間內(nèi)，并輸出象限信息。

附圖5中3為第一角度預測模塊，對調(diào)整后的弧度拆分成三個部分：高位，中位，低位：

對θ^h進行數(shù)學轉(zhuǎn)換，重編碼后一次產(chǎn)生r1,r2,…,r10。各個值的確定方法如下：由于調(diào)整后的弧度一直為正，所以r1＝1，ri＝1-2bi-1。確定r1,r2,…,r10值后，結(jié)合α1,α2,…,α9的值，每次旋轉(zhuǎn)后產(chǎn)生的角度誤差與θ^m+θ^l相加形成θ^m′+θ^l′。

上面提到的α1,α2,…,α9由每次旋轉(zhuǎn)后的誤差值確定，若實際轉(zhuǎn)過的角度大于目標角度，則當前次αi＝-1，若實際轉(zhuǎn)過的角度小于目標角度，則當前次αi＝1。

附圖5中2為第一旋轉(zhuǎn)模塊包含10次旋轉(zhuǎn)。第一級旋轉(zhuǎn)的角度為θ1＝tan^-1(2^-1+2^-5+2^-6)，轉(zhuǎn)換矩陣為：第二級旋轉(zhuǎn)的角度為θ2＝tan^-1(2^-2+2^-8+2^-10)，轉(zhuǎn)換矩陣為第三級旋轉(zhuǎn)的角度為θ3＝tan^-1(2^-3+2^-11+2^-13)，轉(zhuǎn)換矩陣為后面四級到十級旋轉(zhuǎn)角度依次為θ4＝tan^-12^-4，θ5＝tan^-12^-5，θ6＝tan^-12^-6，θ7＝tan^-12^-7，θ8＝tan^-12^-8，θ9＝tan^-12^-9，θ10＝tan^-12^-10＝2^-10，轉(zhuǎn)換矩陣為傳統(tǒng)的矩陣。r1,r2,…,r10由第一角度預測模塊得到后，開始由第一級到第十級依次計算，x/y通路的延遲嚴格保持一樣。在x/y加法的計算中，根據(jù)現(xiàn)有技術，可選擇csa加法器合成的4-2壓縮加法器加速計算，計算完的x/y值輸入第二選擇模塊。

附圖5中5為第二角度預測模塊，第一角度預測模塊累加的剩余角度傳入第二角度預測模塊。對剩余角度的第10到第14位進行與第一預測模塊高位值相同的等式變換，產(chǎn)生r10,r11,r12,r13,r14,r15的值，同時產(chǎn)生低n/2位的剩余角度θ^l″，根據(jù)θ^l″的值產(chǎn)生r15,…,r29的值。r10＝1-2b′9,r11＝2b′10-1,r12＝2b′11-1,r13＝2b′12-1,r14＝2b′13-1,r15＝2b′14-1。r15,…,r29與θ^l″的原碼對應，r15＝b″15,…,，r29＝b″29。

附圖5中4為第二旋轉(zhuǎn)計算模塊包括6次傳統(tǒng)旋轉(zhuǎn)計算和一次乘法計算。第二旋轉(zhuǎn)模塊接收第一旋轉(zhuǎn)模塊的x/y的輸出，接收第二角度預測模塊輸出的r9,r10,…,r15,r15,…,r29的值。r9,r10,…,r15用于確定第10級到第15級傳統(tǒng)旋轉(zhuǎn)的方向。r15,…,r29的值確定15位的乘數(shù)。

附圖5中6為后處理模塊，后處理模塊對第二旋轉(zhuǎn)計算模塊的x/y根據(jù)1中輸出的象限信息進行象限的恢復。最后輸出x/y的值。