鄧凌煊 , 安軍社
(1.中國科學院大學 北京 100190;2.中國科學院 空間科學與應用研究中心,北京 100190)
隨著航天技術的發展,航天任務對于導航計算機的性能要求越來越高。導航計算機除了要對傳感器數據進行采集,與控制系統進行實時通訊,還要能進行實時的計算[1]。盡管目前航天任務中使用的處理器芯片性能越來越強,但大多數CPU并沒有處理常用超越函數(sin,cos,arctan,exp,sqrt,ln 等)的專用指令。而通過純軟件循環迭代求解超越函數往往需要CPU數十甚至上百個周期[2],這極大降低了導航計算機的實時性。本文提出了一種高實時性、低復雜度的CORDIC協處理器核,提供了高吞吐率的超越函數運算能力,從而提高了導航計算機的并行運算能力。此IP核使用verilog編寫,由于其資源占用率低,可以非常容易地被集成入各種航天用FPGA中。
CORDIC算法最早由Volder提出[3],用于求解一般三角函數,之后由Walther改進[4],使得CORDIC可以用于計算雙曲函數和進行乘除運算。由于幾乎所有的通用CPU都具有硬件乘法除法功能,因此對數坐標模式所提供的乘除功能很少被實際使用,故CORDIC算法的主要應用是三角函數和雙曲函數的運算。CORDIC的基本思想是通過一系列預定大小角度的旋轉,使得輸入向量被旋轉到所期望的位置,從而求解出一系列函數。CORDIC的旋轉方程可以表示為

其中m為1時CORIDC工作在圓坐標系下,m為-1時工作在雙曲坐標系下。dn為1時向量按順時針方向旋轉,當dn為-1時按逆時針方向旋轉。
雖然每一步迭代只需要進行加減操作和固定位的移位操作,但由(1)可知,CORIDC每一次旋轉都會改變向量的模。所以為了得到正確的向量分量x和y,需要補償旋轉后的向量,而由于向量模的變化與旋轉的方向無關,這可以通過一個常數乘法器實現。補償常數為:

CORDIC每次旋轉的方向由工作模式決定。在旋轉模式下,CORDIC試圖讓向量與x軸之間的夾角z趨近于0。在向量模式下,CORDIC試圖讓向量的y分量趨近于0。故每次旋轉的方向由運行模式和向量的yz分量符號決定:

綜上,各個模式下CORDIC的迭代結果以及可實現的函數如表1所示。

表1 各模式CORDIC旋轉結果Tab.1 Rotation result of CORDIC under all modes
Cordic能輸出正確結果的前提是旋轉結束時向量能夠被旋轉到預期的位置(在誤差范圍內),即旋轉收斂。可以通過遞歸證明的是[5],對于圓坐標模式,對于任意n總有

令剩余角度為rn+1,

總rn<2·arctan2-n有。但是對于雙曲坐標模式,類似的結論并不成立。然而可以證明的是

對于m=3i+1,i=1,2,3…成立。故對于雙曲坐標模式,為了保證算法收斂,需要將第4,13,40...次循環重復。此外,由于在每一級旋轉的角度是預先固定的,故CORDIC所能將向量旋轉的角度有限,即系的收斂域分別為A∈[-1.743,1.743],對于雙曲坐標模式為A∈[-1.055,1.055]。
為了保證協處理器能夠提供足夠的計算吞吐量并使協處理器的計算延遲可預測,本文使用流水線實現整個CORDIC核[6]。協處理器的流水線結構如圖1所示。

圖1 CORDIC協處理器流水線框圖Fig.1 Pipeline structure of CORDICcoprocessor
協處理器核從輸入FIFO獲得初始輸入數據,包括3個坐標分量以及1個控制命令字。輸入參數的格式為1位符號位,2位整數位,小數位的位數作為IP核的參數可以在例化時調整。這樣的輸入格式使得此IP核可以容易地被應用于使用定點運算的許多DSP處理器。對于浮點數運算,實際上CPU可以通過簡單的移位縮放操作使得輸入范圍外的xyz分量落到協處理器可接受的范圍內,這是由于規格化浮點數的尾數本來就在區間[1,2)內。相對的,已有的很多CORDIC協處理器實現[7]使用了浮點數進行中間運算,然而這不僅顯著地增加了資源的使用,而且使得每一個CORDIC旋轉需要通過多級流水線完成,增大了每個運算的延遲。
此外,為了降低CPU和協處理器之間交互次數,本IP核允許CPU在計算某些函數時不對所有的輸入寄存器進行寫入。對于輸入參數少于3個的函數,協處理器自動生成其他分量的輸入。比如對于cos(x),CPU只需要對協處理器的a0寄存器和控制字寄存器寫入即可觸發cos(x)的運算,CORDIC協處理器會自動把x分量初始化成2.0,y分量初始化成0,z分量初始化成a0。對于有效位數較小的配置如18位,可以進一步將控制命令字和a0放到同一個32位寄存器中,則對于單輸入函數,CPU只需向一個地址寫入數據即可完成操作。18位數據精度時的輸入寄存器格式如圖2所示。

圖2 18位配置時的CORDIC核輸入寄存器Fig.2 Input register of CORDICcore under 18 bits configuration
如表1所示,CORDIC的運算結果并不直接對應所要求的函數,故需要對與輸入參數進行處理。例如對于ln(a)和sqrt(a)運算,需要令 x=a+1,y=a-1,對于 cos(x),sin(x)等運算,需要生成相應的其他分量輸入。此外,由于雙曲坐標的性質,arctanh1并不存在,故雙曲坐標模式只能從i=1開始迭代,而圓坐標系可以從i=0開始迭代,這導致了兩種模式的旋轉過程不同。為了能用同一個流水線實現2種模式的操作,本文令所有模式都從i=1開始迭代。但這樣會導致在圓坐標模式下的收斂域過小,只有。解決的辦法是在預處理單元加入象限折疊,即通過三角函數關系,將[-π,π]上的向量折疊到[0,]上,再在后處理單元對結果進行修正。
旋轉單元是CORDIC協處理器的核心,實現(1)所描述的向量旋轉操作。其結構如圖3所示。

圖3 CORDIC旋轉單元的原理圖Fig.3 Schematic sheet of CORDICrotation unit
此模塊首先通過控制命令字、yz分量計算出旋轉的方向,然后對向量作旋轉。由于圓坐標和雙曲坐標模式下旋轉的角度不同,所以需要根據控制命令字進行選擇。由于使用了流水線結構,移位操作實際上通過布線而靜態確定,不需要專門的移位器。模塊使用3個可進行加減運算的ALU單元分別對xyz分量進行修正。Rom存儲了圓坐標和雙曲坐標模式下的旋轉角度,其精度根據IP核的配置而定,可以簡單地通過在實例化verilog模塊時指定參數來進行配置。
由圖2可知,旋轉單元的輸出結果并不直接對應CORDIC協處理器所提供的函數。后處理單元對CORDIC旋轉單元的輸出結果進行處理,實現所需要的函數。具體的處理如表2所示。

表2 命令字對應的輸出處理Tab.2 Control words and corresponding post processing operations
此外,由于前處理單元對圓坐標輸入進行了象限折疊,后處理單元需要進行相應的修正。例如對于圓坐標旋轉模于是相應的后處理單元需要將x和y的值相交換。
由(4)可知,CORDIC每次旋轉都會改變向量的模,故需要對最終的xy分量進行補償。由于向量模變化只與坐標系模式有關,故補償單元可以用一個常數乘法器實現。常數乘法器可以通過用華萊士樹把移位后幾個向量相加來實現。本文在線下通過程序窮舉找出了一組加減操作數最少乘數編碼方式,對于18位的配置,可以使用一個9輸入的華萊士樹和一個加法器實現,這使得該核在缺乏硬件乘法器的基于flash的Actel FPGA上也可以輕松使用。華萊士樹中一位的結構如圖4所示。

圖4 9位華萊士樹結構Fig.4 The structure of a 9 input Wallace tree
為了證明所提出的IP核的實用性,本文選取了迭代次數16-20,擴展位數為5位的幾種配置進行了綜合。綜合平臺使用了航天電子系統常用的2種 FPGA:A3P3000E和xc4vf40。綜合結果如表3所示。
可見本IP核具有較高的性能和較低的資源占用率,可以較容易地被集成,且隨著迭代次數和精度的提高,資源的增長趨勢穩定,速度的下降并不明顯。

表3 本文提出的CORIDC協處理器實現與基于冗余數實現方案的綜合對比Tab.3 Comparison between proposed CORDIC coprocessor and redundant arithmetic based
一直以來有很多為加速CORDIC運算而提出的方法,其中很多使用了冗余數,試圖通過冗余數進位傳播有限的特點來降低加法器的延遲。但本文通過對比發現,基于冗余數的算法僅在A3P3000E的實現上有一定速度優勢,而在virtex4上沒有速度優勢,還在資源利用率上有很大劣勢。原因是當前的大多數FPGA針對常見運算邏輯,如加法器的實現作出了特殊優化,使用了快速進位鏈降低加法器的進位延遲。相對的,由于FPGA缺乏對冗余數邏輯的支持,冗余數加法器需要由基本單元綜合生成,具有較大的布線延遲,且需要消耗大量的LUT或邏輯門資源。此外,由于冗余數使用2倍于補碼長度的編碼方式來表示數值,冗余數需要使用更多的寄存器寄存每級流水線的中間結果。綜上,對于有效位數較少的配置,本文提出的CORDIC協處理器具有更好的性能與資源利用更具優勢。
根據參考文獻[8]可知,CORDIC計算的誤差主要來自2部分:由于每次迭代時用來修正z分量的角度有限而產生的近似誤差和由于計算位數有限而產生的xy分量截斷誤差。為了盡量降低誤差,首先應該保證旋轉過程中用來修正z的數值是經過舍入得到的,而不是簡單通過截斷得到的,這樣可以減小由于使用近似的旋轉角度而產生的舍入誤差。其次可以考慮在最終階段運算結果前進行舍入。運算結果的舍入可以在補償單元通過增加華萊士樹的輸入實現。即增加一個輸入變量,其大小為最低有效位的一半,符號與補償單元的輸入相同。除此之外,需要根據所需的精度合理選擇CORDIC協處理器的數據位數,增加擴展位,即對于n位的CORDIC運算,在迭代運算時使用n+m位,在迭代結束后舍棄m位,這樣做可以保證前面的n位不受到截斷誤差的影響.。本文選取了18位CORDIC來說明擴展位和迭代次數對精度的影響。對于圓坐標旋轉模式,輸入向量為x=2.0,y=0.0,z∈[-π,π]的所有可能輸入向量;對于圓坐標向量模式,輸入為模為2.0,夾角屬于[-π,π]的所有可能輸入向量;對于雙曲旋轉模式,輸入為 x=2.0,y=0.0,z∈[-1.0,1.0]的所有可能輸入向量;對于雙曲向量模式,輸入為 x=a+1,y=a-1,α∈[0.25,1.0]的所有可能輸入。各個函數的誤差與所選取的迭代次數和擴展位的關系如表4所示。

表4 幾種配置下的各函數的計算誤差,單位為ulpTab.4 Error of functions under several configurations
其中ulp為最低位單位。對于20位定點數有1ulp=2-18。配置為擴展位數和迭代次數。通過實驗可以發現,對擴展后的向量進行截斷會帶來顯著的誤差增長。并且只會在使用較長擴展位時對結果進行舍入才有明顯的效果,因為擴展位較少的情況下舍入反而可能會讓結果向錯誤的方向舍入。最后值得注意的是雙曲模式下的函數誤差高于圓坐標模式,由(8)可知這是由于雙曲模式向量角度收斂性更差而導致的。
本文提出的CORDIC協處理器可以容易地集成進目前常見的航天級FPGA中,為CPU提供更強的三角函數和超越函數運算能力。在中端的V4系列FPGA中實現萬分之1精度的三角函數和超越函數只需要不到10%的資源,并可運行于較高的系統時鐘下。CORDIC協處理器提供了并行運算幾種常用三角和超遠函數的功能,不僅適用于導航計算機,也可以被用于其他有大量實時性計算需求的嵌入式系統中。
[1]閆捷,徐曉蘇,李瑤,等.基于DSP與FPGA的嵌入式組合導航計算機系統設計 [J].測控技術,2013,32(12):61-64 YAN Jie,XU Xiao-su,LI Yao,et al.Design of embedded integrated navigation system based on DSP and FPGA[J].Measurement&Control Technology,2013,32(12):61-64.
[2]Texas Instrument.TMS320C54x DSP Library Programmer’s Reference[EB/OL].(2013).http://www.ti.com.cn/lit/ug/spru518d/spru518d.pdf.
[3]Volder JE.The birth of CORDIC[J].Journal of VLSI Signal Processing Systems,2000,25(2):101-105.
[4]Walther J S.The story of unified CORDIC[J].Journal of VLSI Signal Processing Systems,2000,25(2):107-112.
[5]Jean-Michel Muller.Elementary Functions Algorithms and Implementation Second Edition[M].Boston:Birkh?user,2005.
[6]王思聰,文治平,于立新.空間用CORDIC處理器的結構級設計方法[J].微電子學與計算機,2006,23(8):57-60.WANG Si-cong,WEN ZHi-ping,YU Li-xin.The architecture desing method of CORDIC processor used in space[J].Microelectronics&Computer,2006,23(8):57-60.
[7]Munoz DM,Sanchez DF,Llanos CH,et al.FPGA based floating-point library for CORDIC algorithms[A].Programmable Logic Conference (SPL) 2010 VI Southern[C]//Ipojuca,2005:55-60.
[8]Sarbishei O,Radecka K.On the Fixed-Point Accuracy Analysis and Optimization of FFT Units with CORDICMultipliers[A].Computer Arithmetic (ARITH),2011 20th IEEE Symposium[C]//Tubingen,2011:62-69.