文獻標識碼： A
CORDIC算法(Coordinate Rotation Digital Computer)(即“坐標旋轉數字計算機"),最早是由VOLDER J E.于1959年提出[1]，當時是為研制B-58轟炸機的導航系統而設計的。但該算法僅僅用到了移位和加減運算，因此其硬件實現非常簡單。后經眾多學者加以研究和發展，使其成為數字電路中計算各種超越函數的一種簡捷有效的算法[2]。CORDIC算法在FPGA中應用甚廣[2-3]，在DSP上也有所應用[4]，但鮮有在單片機上的應用報道。
1 CORDIC算法簡介
CORDIC算法的基本原理來自二維坐標的旋轉變換。當平面直角坐標系繞原點旋轉一個角度θ時，新舊坐標間的變換關系為：
以上介紹的CORDIC算法屬于最基本的“圓周型”CORDIC算法。對迭代步驟稍作改動，即可將其變形為“直線型”或“雙曲型”CORDIC算法，從而可用于計算乘除法、雙曲函數、指數函數和對數函數[2]。
CORDIC算法中的移位運算很適合于在FPGA中用硬件實現，不但計算速度快，而且相對于查表法明顯節省了硬件資源。而在DSP中，由于缺少專門的硬件支持，CORDIC算法相對于冪級數展開法或查表法則缺乏優勢[2-4]。對于經典的PIC、8051等單片機，因為每執行1條移位指令只能將數據左移或右移1位，而執行CORDIC算法時會在移位運算上消耗大量的時間,從而影響CORDIC算法的效率。
2 羅盤載體姿態及其解算
式(7)、式(8)中的Bxh與Byh即是地磁場在水平面內的2個正交分量，分別沿xh軸和yh軸方向。式(7)所描述的正是將載體坐標系中測得的地磁分量變換到地平坐標系的過程。
參考文獻[5]給出了3種求解羅盤載體姿態的算法，上文介紹的就是其中的第1種，而第2種算法原理與其完全相同,只是利用三角函數的性質略微簡化了計算。應當指出,由于坐標系及角度定義上的差異，式(5)、式(7)、式(8)與參考文獻[5]中的公式略有不同，但并不影響問題的本質。
CORDIC算法本身可以計算三角函數與反三角函數，因此可以直接將其用于上述算法中[6]，但這樣并沒有完全發揮CORDIC算法的優勢。如前所述，CORDIC算法的原理來自于坐標旋轉變換,因此,可以直接利用CORDIC算法完成式⑺中的變換并由此獲得載體的姿態，而無需計算三角與反三角函數。參考文獻[5]所給出的第3種算法就是基于這一思想，其步驟為：
(1)用矢量模式將重力加速度矢量在yOz平面內的投影旋轉到與z軸重合，求得橫滾角。
(2)再用矢量模式將重力加速度矢量（此時該矢量已在xOz平面內）旋轉到與z軸重合，求得俯仰角。
(3)用旋轉模式補償橫滾角，計算補償后的地磁場分量。
(4)再用旋轉模式補償俯仰角，計算補償后的地磁場分量。
(5)用矢量模式，根據地磁場水平分量計算航向角。
3 CORDIC算法在單片機上的實現
3.1 用CORDIC算法計算三角函數
表1給出了在VRS51L3074單片機上用CORDIC算法計算反正切函數的結果,并與Keil C51編譯器（V8.12）所提供的浮點庫函數進行了對比。VRS51L3074是美國Ramtron公司的產品，基于增強型8051內核，在兼容傳統8051指令集的基礎上提升了指令執行速度。同時，該單片機還擁有“增強型算術單元”的外圍模塊，該模塊中有32位桶形移位器。表1也列出了在CORDIC算法中使用該移位器的結果。
此處所用CORDIC算法的輸入、輸出均采用16位數據格式，角度以度為單位，小數點設在兩字節之間，即角度分辨率為1/256度，迭代次數為14次(不含第1次45°的旋轉），這是由角度的分辨率決定的（第15次及其后的迭代所轉角度已小于1/256度）。為便于比較，所有計算結果均已換算為十進制格式(單位為度，保留4位小數)。表1還列出了由單片機的定時器T0所記錄的函數執行時間，T0時鐘為系統時鐘(不分頻)。因此,計時結果即為函數從調用到返回所經歷的系統時鐘周期數。
從表1的結果可看出，CORDIC算法具有一定的速度優勢，但其精度無法與浮點函數相比。另外，使用桶形移位器代替單片機的移位指令，可明顯提升CORDIC算法的執行速度。這也證明了CORDIC算法在單片機上運行時，在移位操作上消耗了大量的時間。CORDIC算法的主要優點在于： (1)代碼較為簡潔；(2)可自定義輸出角度單位及格式(浮點庫函數所給出的角度單位總是弧度)； (3)在計算y/x的反正切時不必計算除法，并且該算法在±90°時也是收斂的。
在單片機上運行CORDIC算法還需注意，輸入矢量的模長既不能過大也不能過小。如模長過大，則由于CORDIC算法在將矢量旋轉的同時還會將其“拉長”(這也正是CORDIC算法需要模校正因子的原因)，將會導致計算過程出現溢出。而若模長過小，則會導致精度降低[7]。表2為輸入矢量的模長對CORDIC算法的影響：在4個輸入矢量的輻角均相等的情況下，第1個矢量由于模長過大導致CORDIC計算過程溢出，第2個矢量計算結果正常，第3和第4個矢量由于模長過小導致CORDIC計算精度降低。
移位操作耗時過多以及計算精度受輸入矢量模長影響這兩個缺點，使得用CORDIC算法在單片機上計算三角函數的實用價值并不理想。
3.2 CORDIC算法用于姿態解算
如上所述，CORDIC算法并不適合于在單片機上直接計算三角函數。但將其用于羅盤載體的姿態解算仍然是可行的，這是基于以下兩點：(1)由于傳感器特性的不一致性，其采集的原始數據均需經過零點、增益等校準后才可用于姿態解算，這就相當于為CORDIC算法的輸入矢量模長做了歸一化處理；(2)如前所述,利用CORDIC算法計算俯仰角與橫滾角并對其加以補償時，并不需要具體計算三角與反三角函數，只需完成坐標變換即可。
參考文獻[5]不僅給出了用CORDIC算法求解航向角的思路，并且在FPGA上進行了驗證。將這一算法移植到單片機時，還可再稍作簡化，具體計算步驟為：
(1)將重力加速度矢量在yOz平面內的投影旋轉到與z軸重合，同時地磁場矢量亦作同步的旋轉，從而得到橫滾角并同時對其加以補償。
(2)將重力加速度矢量(此時該矢量已在xOz平面內)旋轉到與z軸重合，同時地磁場矢量亦作同步的旋轉，從而得到俯仰角并同時對其加以補償。
(3)根據地磁場水平分量計算航向角。
以上計算過程相當于把參考文獻[5]所述算法步驟的第1、第3步和第2、第4步分別合并。與CORDIC算法的矢量模式和旋轉模式在FPGA上的實現方式不同，在單片機上則不必拘泥于這兩種模式的區分，從而可以使程序更加緊湊。
表3給出了在VRS51L3074上使用增強算術單元實現CORDIC求解載體姿態的實驗結果，并同時給出按式(5)~式(8)利用浮點運算所得結果進行對比。磁場矢量的模長已被限制到4 096(相當于十六進制的0x1000)，重力加速度矢量的模長則限制為8 192(相當于十六進制的0x2000)。
電子羅盤的俯仰角、橫滾角和航向角要求達到0.1°的分辨率。由表3可見，CORDIC算法的計算精度能夠滿足這一要求，并且該算法完成1次姿態解算所需的時間僅相當于浮點運算所需時間的1/8～1/9，能很好地滿足電子羅盤對姿態解算的實時性要求。
CORDIC算法無需乘除運算和大規模的查找表即可實現三角函數、反三角函數以及其他超越函數的運算。盡管CORDIC算法通常被用于FPGA，實踐結果證明,這一算法也可有效地在8051內核單片機上運行，并且具有優于C語言浮點庫函數的運行速度和一定的精度。
在單片機上，制約CORDIC算法效率的主要因素是移位操作耗時過多。若單片機擁有桶形移位器，如VRS51L3074單片機的增強型算術單元，則可改善移位操作耗時過多的問題，從而提升該算法的效率。此外，適當調節輸入矢量的模長也是應用該算法時需要注意的問題。
采用合理安排的計算步驟，CORDIC算法可以準確、高效地完成電子羅盤中的傾角補償與航向計算，同時也為電子羅盤中的單片機節約了程序空間和時間。
參考文獻
[1] VOLDER J E. The birth of CORDIC[J]. Journal of VLSI Signal Processing, 2000,25:101-105.
[2] ANDRAKA R. A survey of CORDIC algorithm for FPGA based computers[A]. Proceedings of the 1998 ACM/SIGDA International Symposium on Field Programmable Gate Arrays[C]. ACM, 1998:191-200.
[3] 李全，陳石平，李曉歡，等. 正交三角函數的CORDIC實現[J]. 微計算機信息，2008，24(12-3)：268-269，252.
[4] 馬士超，王貞松. 基于DSP的三角函數快速計算[J]. 計算機工程，2005，31(22)：12-14.
[5] LAULAINEN E, KOSKINEN L, KOSUNEN M, et al. Compass tilt compensation algorithm using CORDIC[A]. IEEE International Symposium on Circuits and Systems, 2008[C]. ISCAS, 2008:1188-1191.
[6] 崔曉松，胡建萍，李陬. CORDIC 算法在導航解算系統中的應用[J].杭州電子科技大學學報，2007，27(6)：5-8.
[7] KOTA K, CAVALLARO J R. Numerical accuracy and hardware tradeoffs for CORDIC arithmetic special-purpose processors[J]. IEEE Transactions on Computers, 1993, 42(7):769-779.