一種降低CORDIC算法截斷誤差的方法

曲世雋，王翾

（中國傳媒大學信息與通信工程學院，北京 100024）

1 引言

坐標旋轉數字計算機（Coordinate Rotation Digital Computer,CORDIC）算法是一種用于高效計算基本函數的迭代算法。該算法通過一般的加減和移位操作實現了乘除法的相關計算，使得向量的的旋轉和定向的計算不再依賴于三角函數、乘法、開方、反三角、指數等復雜函數。CORDIC 算法的硬件實現只需要移位器和加法器，簡單的硬件就能有效的運算復雜函數。

在1959年，Volder開發了一類計算三角函數、雙曲函數的算法，其中包括指數運算和對數運算，這就是CORDIC算法的雛形。在提出這個算法之后，Volder將其應用在導航領域之中。但是在當時這個算法并沒有被廣泛的討論和應用。直到1980 年，Haviland 和Tuszynski[1]給出了一種全模式的CORDIC算法的處理芯片。在這之后，CORDIC算法才引起人們的關注，被應用于各種領域。Hu Y H[2]將CORDIC算法應用于基于VLSI的數字信號處理領域。應用包括離散傅里葉變換的計算、矩陣運算和三角函數發生器。在2009 年，Vachhani L,Sridharan K,Meher P K[3]通過FPGA實現高效CORDIC 算法并將其應用在機器人探測上。當然，CORDIC算法也可以應用于其他領域。例如，計算機制圖中求點到線的距離、直角坐標與極坐標的相互變換及求多維向量的歐幾里得范數等，可以預見，CORDIC算法將在未來有更廣泛的應用。

目前人們對于CORDIC算法的關注點主要集中在硬件資源消耗和計算精度兩個問題上。想要獲得高計算精度就必然要消耗更多的硬件資源，因此如何平衡兩者關系就成了一個重要問題。在文獻[4]中，提出了一種在硬件資源和計算精度之間尋求折中的方法。文中介紹了兩種解決CORDIC算法精度問題的方法，第一種建立在定點數CORDIC算法的運算基礎上，文中認為輸入變量的不規范會使誤差變大，因此需要對輸入的變量進行統一、規范的處理，以此提高CORDIC算法的精度，第二種方法是針對于浮點數的CORDIC算法，通過使用混合體系結構降低復雜度，以此提高精度。在實際應用中，為了硬件結構的簡單和運算的高效，我們更加偏向于使用定點數進行CORDIC算法。但是文中介紹的提高定點數精度的CORDIC算法在硬件結構上還是過于復雜，因此需要提出一種硬件結構簡單，更容易實現的方法來降低誤差。

為了提出一種更好的方法實現硬件資源和計算精度的平衡，就需要對CORDIC 算法的誤差來源進行詳細的分析。CORDIC 算法的誤差可以分為近似誤差和截斷誤差，在任何坐標系的任何模式中，近似誤差是由有限的迭代次數產生的，而截斷誤差是由有限的數據位位寬產生的[5]。Tze-Yun Sung,Yi-Hsun Sung[6]對所有坐標系下CORDIC算法的誤差都進行了分析。包括圓坐標系，雙曲坐標系，線性坐標系的旋轉模式和矢量模式，用近似的方式給出了一個誤差的最大值。但是文獻主要是基于數學分析的角度，給出了近似誤差和截斷誤差的分析結果。在實際CORDIC 算法的應用中，需要有硬件結構的設計和反正切、正余弦函數的計算結果，這些在文章中沒有提到。

本文提出了一種降低CORDIC算法中截斷誤差的方法，將CORDIC算法的流程分為預處理、迭代移位、后處理三個部分，對每個部分進行不同的處理，旨在使用較少的硬件資源同時保持較高的計算精度。文中給出了基于此種方法的硬件架構設計，并做了基于FPGA的仿真實驗，驗證了結果的可行性和正確性。

2 CORDIC算法原理

CORDIC算法的基本公式如下所示[5]：

其中s(m,i)是非遞減的移位序列，滿足

是第i次旋轉的角度。m= +1,- 1,0確定坐標系是圓，雙曲或是線性。δ(i)確定旋轉角的方向，+1和-1分別代表不同的方向。由于CORDIC 算法的實現方式很多，以圓坐標系為例，分旋轉模式和矢量模式進行原理分析。

通常在旋轉模式下，設定x(0) = 1/k1,其中km=，n 為迭代次數。y(0) = 0。z(0)為輸入待計算的角度θ。旋轉方向δ(i)由z(i)的正負決定，若z(i)為正，則δ(i)為+1，若z(i)為負，則δ(i)為-1。根據公式(1)(2)(3)(5)可以得到角度θ的正弦值y(n)和余弦值x(n)。

通常在矢量模式下，設定輸入的值x 和y，通過迭代和旋轉，使y(i)逐漸逼近于0。旋轉方向δ(i)由y(i)的正負決定，若y(i)為正，則δ(i)為-1，若y(i)為負，則δ(i)為+1。同樣根據公式(1)(2)(3)(5)，可以得到反正切函數z(n) = arctan(y/x)，得到向量(x,y)的長度d=x(n)*(1/k1)。

3 CORDIC算法的誤差

在研究CORDIC算法的時候，我們需要在硬件資源的使用和計算精度上做出平衡，這就需要對CORDIC算法的誤差進行分析。CORDIC算法的誤差由迭代次數n，旋轉角度小數位位數b，x和y的小數位位數bv確定。在本文中b和bv采用相同的大小。CORDIC算法的誤差可以分為兩部分：近似誤差和截斷誤差。EA表示近似誤差最大值，由有限的迭代次數造成，ER表示截斷誤差最大值，由有限的數據位位寬造成。

根據文獻[6]，圓坐標系下旋轉模式的近似誤差為

截斷誤差為

矢量模式下的近似誤差為

截斷誤差為

基于以上的理論推導，可以確定迭代次數和小數位位數對于CORDIC算法的誤差的影響，如圖1和圖2所示。

圖1 圓坐標系旋轉模式誤差

圖2 圓坐標系矢量模式誤差

4 降低截斷誤差的方法

4.1 方法簡述

從前文的分析中可以得知，增加數據位位寬和迭代次數可以顯著的提高CORDIC 算法的精度，但是同時也會消耗更多的硬件資源。降低近似誤差，只需要增加迭代的次數，但是降低截斷誤差卻不能只是簡單的增加數據位位寬。過于長的數據位位寬不僅會消耗更多的硬件資源，同時也會對后續的數據處理造成影響。因此，需要有一個合適的方法來降低截斷誤差。

通過此種方式實現CORDIC 算法，完整的步驟可以分為預處理、迭代移位、后處理三個部分。三個步驟如圖3所示。

圖3 實現改進CORDIC算法完整步驟

為了有更高的精度，我們選擇在預處理部分進行數據的規范化處理和擴充數據位位寬。在迭代移位的過程中如果改變數據位位寬，很可能會產生溢出問題，因此選擇在后處理部分對迭代移位后的數據進行截斷。這樣就可以保證在迭代移位部分提高了計算精度，同時輸出位寬又可以根據實際的需要進行截取。

4.2 預處理

預處理部分可以分為對輸入數據進行規范化處理和對輸入數據的位寬進行擴充兩個模塊。假設輸入數據x_in和y_in為16bit 定點數，最高位為符號位。對輸入數據進行規范化處理要計算除符號位之外的先導零的個數。比較x_in和y_in的先導零個數，j為兩者先導零個數的較小值，根據公式

可以得到j的值。其中m為輸入數據位寬，ε= 2-b為精度。在硬件設計中，可以通過一個15-4優先編碼器和一個逐位比較的或門來得到j的值。x_in和y_in并行輸入逐位比較的或門，可以得到兩者之間較高位的1所在的位置。再通過15-4優先編碼器得到較高位的1之前先導零的個數。在得到j之后，統一將x_in和y_in左移j位。為了防止迭代過程中發生溢出，還要在符號位后增加兩個0 位。根據計算，可以得到k1的值在1.414 和1.646 之間，因此擴充兩個0 位完全可以避免溢出。15-4優先編碼器的真值表如表1所示。

表1 15-4優先編碼器真值表

對輸入數據的位寬進行擴充就是在輸入數據末尾補零。我們將輸入數據擴充為32bit。也就是在進行完數據的規范化處理之后需要在末尾補14 位的0。預處理硬件設計如下圖4所示。

圖4 預處理硬件設計

4.3 迭代移位

CORDIC 算法的單層迭代移位結構如圖5 所示。它需要1 個多路復用器，3 個加/減法器，2 個移位器。旋轉方向由selx,sely,selz定義，多路復用器根據旋轉模式或矢量模式來進行選擇。當模式為旋轉模式時，旋轉方向由sign(zi)確定。當模式為矢量模式時，旋轉方向由sign(yi)確定。要實現并行流水線的CORDIC 算法，只需要對單層的結構進行疊加就可以實現。疊加的次數取決于迭代次數N。

圖5 單層迭代結構

4.4 后處理

后處理的部分實際上就是對迭代移位后輸出的數據進行位寬的縮減，本質上也是對前導零的去除。縮減的原則是，首先從符號位后最高位開始判斷，如果最高位為0，則向左移1 位，繼續判斷下一位，若還為0，則繼續向左移1位，判斷下一位，直到遇到為1的位數，則停止判斷。對0的位數舍去完之后，整體保留16bit。設后處理部分總共左移的位數為num。后處理部分流程圖如圖6所示。

圖6 后處理流程圖

通過流程圖可知，實現后處理部分需要通過循環結構來進行判斷，為了節省硬件資源，可以使用計數器來實現循環結構。因此，后處理部分需要1 個累加器和1 個計數器，1 個寄存器。將需要處理的數據放入寄存器，在每個計數器+1時進行逐位判斷。后處理結構如圖7所示。

圖7 后處理結構

最后可以總觀整個流程中對數據的左移，可以得到數據在預處理部分左移的位數為14 +j,在后處理部分左移的位數為num- 16。所以整個流程對數據的放大為2j+num-2倍。也就是說得到的數據精度為ε=2j+num-2,在改進后的CORDIC 算法中，可以在預處理和后處理部分去除無意義的前導零，以此提升精度。

5 仿真與測試

在本設計中采用Spartan 3E 開發板，用Verilog 對CORDIC 算法在圓坐標系下的矢量模式進行了描述，采用流水線結構。對于100000 組輸入范圍在(0,10000)的獨立均勻分布的16bit 數據進行基于FPGA的圓坐標系矢量模式下的CORDIC 算法計算。選擇迭代次數為16 次，在預處理時將其擴充為32bit，在后處理時保留16bit。為了驗證其正確性，對其實現過程用ModelSim10.5進行仿真，仿真結果如圖8所示。

圖8 ModelSim仿真結果

從仿真結果中可以看出輸入和輸出結果之間存在延時，這跟預處理和迭代移位過程有關。

提取仿真過程中激勵文件隨機生成的100000 組輸入值，對比理想輸出相位值、幅度值和實際輸出相位值、幅度值之間的絕對誤差。其絕對誤差的分布如圖9所示。

圖9 使用截斷方法的CORDIC算法的絕對誤差分布

計算得到幅度值的平均絕對誤差為1.67*10-2，方差為1.3956*10-4。相位值的平均絕對誤差為1.1*10-3，方差為2.0246*10-6。對比使用相同數據位位寬的傳統CORDIC算法，在精度上有了一定的提升。

6 結束語

該文簡要介紹了CORDIC 算法的原理，從硬件實際應用角度出發，提出了一個降低CORDIC 算法截斷誤差的方法。實現了對復雜函數的高精度運算，并給出了圓坐標系矢量模式下的仿真結果。在實驗過程中，可以通過改變迭代的次數或是數據位位寬來滿足不同的精度要求。從仿真結果來看，提出的方法可以很好的在精度要求和硬件資源消耗上達到平衡。這個方法在本文中只給出了在圓坐標系矢量模式的應用分析，在以后的研究中，可以將這個方法推廣到圓坐標系、線性坐標系、雙曲坐標系的不同模式中。為CORDIC 算法計算開方、反三角等復雜函數問題提供了一個新的處理方法，具有一定的借鑒意義。