三角函數的計算是個復雜的主題，有計算機之前，人們通常通過查找三角函數表來計算任意角度的三角函數的值。這種表格在人們剛剛產生三角函數的概念的時候就已經有了，它們通常是通過從已知值（比如sin(π/2)=1）開始並重復應用半角和和差公式而生成。
現在有了計算機，三角函數表便推出了歷史的舞台。但是像我這樣的喜歡刨根問底的人，不禁要問計算機又是如何計算三角函數值的呢。最容易想到的辦法就是利用級數展開，比如泰勒級數來逼近三角函數，只要項數取得足夠多就能以任意的精度來逼近函數值。除了泰勒級數逼近之外，還有其他許多的逼近方法，比如切比雪夫逼近、最佳一致逼近和Padé逼近等。所有這些逼近方法本質上都是用多項式函數來近似我們要計算的三角函數，計算過程中必然要涉及到大量的浮點運算。
在缺乏硬件乘法器的簡單設備上（比如沒有浮點運算單元的單片機或FPGA），用這些方法來計算三角函數會非常的費時。為了解決這個問題，J.Volder於1959年提出了一種快速算法，稱之為CORDIC(Coordinate Rotation Digital Computer) 算法，這個算法只利用移位和加減運算，就能計算常用三角函數值，如Sin，Cos，Sinh，Cosh等函數。
CORDIC可以應用於圓周系統、線性系統和雙曲系統等，在不同的系統內解決不同的復雜計算問題。圓周系統中解決三角函數計算問題，線性系統中解決乘法和除法問題。
1. 基本思路
通常，旋轉模式常用來解決三角函數的問題，體現的應用主要是極坐標向直角坐標轉換，即已知一點的極坐標(α,r)，求其直角坐標(x,y)，實際上是求sinα、cosα或者tanα。
為了便於理解，開頭先拋出cordic算法的實質，即類似於二分法加旋轉的概念，下面來解釋。
如圖1.1所示，在單位圓上，向量OP與X軸的正半軸夾角為α，故P點的坐標可表示為：
將向量OP順時針旋轉θ角至向量OQ，此時OQ與X軸正半軸的夾角為α-θ，舉個例子來簡單說明一下思想，若P點對應的極坐標為(51°,1)，要轉換為直角坐標系，我們需要知道cos51°和sin51°，假設我們不能借助其他的手段來直接計算該值，那么我們逐次的旋轉一定的角度來逼近51°，而旋轉的這些角度值的選擇有一定的特殊性，一會再解釋。假設我們從x軸開始逆時針旋轉（或者從51°開始順時針旋轉），旋轉45°，發現沒到51°，再逆時針旋轉45°/2=22.5°，則一共旋轉了67.5°，超了，再順時針旋轉22.5°/2=11.25°……逐次逼近51°，這個過程用公式描述為：
C語言描述過程如下，以及模擬旋轉結果：
#include <stdio.h>
#include <stdlib.h>
double my_tan2(double z);
int main(viod)
{
double x = my_tan2(51.0);
return 0;
}
double my_tan2(double z)
{
const double theta[] = { 45.0, 26.56505118, 14.03624347, 7.125016349,
3.576334375, 1.789910608, 0.8951737102, 0.4476141709,
0.2238105004, 0.1119056771, 0.05595289189, 0.02797645262,
0.01398822714, 0.006994113675, 0.003497056851, 0.001748528427
};
int i = 0;
double angle_new = 0;
double angle_remain = z;
char detection;
for( i=0; i<16;i++)
{
if(angle_remain > 0)
{
angle_remain = angle_remain - theta[i];
detection = '+';
}
else
{
angle_remain = angle_remain + theta[i];
detection = '-';
}
printf(" 旋轉次數 = %-8d 旋轉角度 = %-12f 旋轉方向：%-8c 剩余角度 = %-8f\n", i, theta[i],detection, angle_remain);
}
return angle_remain;
}
再精度滿足要求的前提下，旋轉一定從次數，則可以求得相應點的直角坐標。有了舉例鋪墊，再從理論上分析過程，以及選擇旋轉角度的依據。
2. 減少乘法運算
圖中，Q點的坐標可表示為：
這里定義θ為目標旋轉角度。根據三角函數公式可將上式展開為：
現在已經有點 Cordic 算法的樣子了，但是我們看到每次旋轉都要計算 4 次浮點數的乘法運算，運算量還是太大了。還需要進一步的改進，改進的切入點當然還是坐標變換的過程。
將式(1.1)代入到式(1.3)中可得
提取cosθ，上式可重新表示為：
從上式中可看出，cosθ只是改變了目標向量OQ的模長，每次旋轉后的新坐標點到原點的距離都變長了，如果去掉cosθ，縮放系數就為1/cosθ，圖1.2所示，此時OP旋轉θ角之后到達OR，這種旋轉稱為偽旋轉。
不難看出，OQ與OR幅度上的差異，且可證明向量OP和PR是正交的，此時R點的坐標可表示為:
我們求的是θ，不關心模長r的改變，這樣的變形非常的簡單，最終可以通過對偽旋轉的輸出加以補償來獲取真實的旋轉結果。
對於偽旋轉，可以將θ分解為一系列滿足要求的小角度和，即
這樣，將一次旋轉分解成了一系列的微旋轉。那么第i+1次旋轉后的結果為
但是每次循環的運算量一下從4次乘法降到了2次乘法了。同理，在編程過程中，可以利用查表法，事先將如45°、22.5°等角度的正切值存入表中，在計算過程中，對照旋轉值，最終解得極角的三角函數，從而得到某極坐標的直角坐標。
3. 消除乘法運算
這一步至關重要，其實就是cordic算法能在FPGA上應用的關鍵所在。上述分析中，我們已經成功的將乘法的次數減少了一半，但是畢竟還是有乘法運算，以為這在硬件上實現還是會消耗很多資源，那么還有沒有可能進一步降低運算量呢？還要從計算式入手。
第一次循環時，tan45°=1，所以第一次循環實際上是不需要乘法運算的。第二次運算呢？
tan22.5°=0.4142135623731，第二次循環乘數是小數，選擇22.5°是因為二分查找法的查找效率最高。如果選用個在22.5到45度之間的值，通過犧牲查找效率而減少乘法，這是可以考慮的。
如果令
d_i代表方向，取+1時代表順時針，取-1代表逆時針。(1.8)式，可重新改為
從上式不難看出，每次微旋轉只需加法、減法和移位操作（乘以2的-i次冪，即右移i位）即可完成。
我們發現tan26.565051177078°=0.5，如果我們第二次旋轉采用26.565051177078°，那么乘數變為0.5，如果我們采用定點數運算的話（沒有浮點協處理器時為了加速計算我們會大量的采用定點數算法）乘以0.5就相當於將乘數右移一位。右移運算是很快的，這樣第二次循環中的乘法運算也被消除了。類似的方法，第三次循環中不用11.25°，而采用 14.0362434679265 °。tan14.0362434679265°=1/4，乘數右移兩位就可以了。剩下的都以此類推。
為了確定d_i的值，引入一個新的變量z，定義
因為旋轉可以從兩個方向進行，可以從0°旋轉來逼近目標極角α，也可以從目標極角α旋轉來逼近0°角。若設z初始化為0°，即z_0 = 0。根據條件對z_i執行加或者減tan^(-1) 2^(-i)，使得z的最終值逼近為α。z的迭代過程，是將z與α的差值收斂於0的過程，也正是將旋轉角度θ分解為一系列θ_i的過程，過z_i可認為是第i次旋轉之后的終點角度。至此，cordic算法之圓周系統的旋轉模式迭代過程可表示為：
旋轉角度θ的最大值和最小值，可計算得
即表明，算法只能支持{-99.8829°，99.8829°}范圍內的角度處理。在很多場合中需要目標旋轉角度可以覆蓋[-180°,180°]，這就需要預處理。只有當目標角度|α|>π/2時需要預旋轉π/2。
上述討論都是基於偽旋轉，即每次旋轉模長都有變化，省掉了cosθ_i。以k_i表示第i次微旋轉模長補償因子，故第i次微旋轉真實的旋轉結果應為
模長補償因子可以建立一個查找表，根據旋轉次數定，此時不得不增加乘法，但是若旋轉次數為n，則總的模長補償因子k可表示為：
當n趨於無窮大時，K逼近0.607252935。所以當旋轉次數定了之后，可以在最后的結果只增加一次乘法，來補償模長，從而得到真實的直角坐標，至此就完成了極坐標到直角坐標的轉換。滿足要求的角度和模長補償因子可以用matlab生成，附程序：
matlab生成旋轉角度和模長補償因子
clear all
clc
i = 15; %旋轉階數-1
for n = 0:i;
theta(n+1) = atand(1/(2^n)); %生成符合要求的角度
k(n+1) = 1/sqrt(1+(2^(-2*n))); %生成模長補償因子
end
vpa(theta,10); %控制精度
vpa(k,10);
round(theta); %取整
附C程序來理解過程，注意我采用的方法是從0°（x軸）逐次旋轉逼近目標角度，注意代碼還是以浮點數運算，稍后會根據FPGA等暫時適合做定點處理的平台進行微改。
#include <stdio.h>
#include <stdlib.h>
double cordic_c(double a,double r);
double x = 1, y = 0; //以X軸為旋轉起始點
int main(viod)
{
double remain = cordic_c(31.3,15.0); //極坐標值（極角，極徑）
printf("旋轉角度誤差：%f, 直角坐標：x = %f, y = %f",remain,x,y);
return 0;
}
double cordic_c(double a,double r)
{
const double theta[] = { 45.0, 26.56505118, 14.03624347, 7.125016349,
3.576334375, 1.789910608, 0.8951737102, 0.4476141709,
0.2238105004, 0.1119056771, 0.05595289189, 0.02797645262,
0.01398822714, 0.006994113675, 0.003497056851, 0.001748528427
}; //旋轉角度
const double k[] = { 0.7071067812, 0.894427191, 0.9701425001, 0.9922778767,
0.9980525785, 0.9995120761, 0.999877952, 0.9999694838, 0.9999923707,
0.9999980927, 0.9999995232, 0.9999998808, 0.9999999702, 0.9999999925,
0.9999999981, 0.9999999995
};//模長補償因子
int i = 0;
double x_temp = 1.0, y_temp = 0.0;
double angle_new = 0.0; //旋轉后終止角度
double angle_remain = a; //旋轉后，剩余角度
char detection; //旋轉方向
for( i=0; i<16;i++)
{
if( angle_remain == 0)
{
break;
}
if(angle_remain > 0)
{
angle_new = angle_new + theta[i];
angle_remain = a - angle_new;
x_temp = (x - y / (1<<i))*k[i];
y_temp = (y + x /(1<<i))*k[i];
x = x_temp;
y = y_temp;
detection = '+';
}
else
{
angle_new = angle_new - theta[i];
angle_remain = a - angle_new;
x_temp = (x + y /(1<<i))*k[i];
y_temp = (y - x /(1<<i))*k[i];
x = x_temp;
y = y_temp;
detection = '-';
}
printf(" 旋轉次數 = %-8d 旋轉角度 = %-12f 旋轉方向：%-8c 終點角度 = %-8f\n", i+1, theta[i],detection,angle_new);
}
x = r*x;
y = r*y;
return angle_remain;
}
從結果中，看到旋轉精度能接受。如用查表法的話，其實是空間換時間，速度最快，但是想要精度高的話會占用大量的存儲空間（你不可能每個角度都占用一個空間）。 cordic 屬於比較均衡的算法，空間占用很少，計算時間也挺快。
值得注意的是，理論分析過程是針對第一象限的旋轉，如果極坐標在其他象限，則首先有個預處理過程，將極坐標先轉到第一象限后，再進行cordic算法。處理完之后再做相應的后處理。
參考
-
《基於FPGA的數字信號處理（第2版）》——高亞軍著