1、CORDIC算法
CORDIC（Coordinate Rotation DIgital Computer）是一種經過迭代實現快速平面旋轉的算法，經過變形擴展，它能夠對多種數學函數求值，例如求解三角/反三角函數、雙曲函數等，除此以外還能夠作特殊的向量運算甚至算術運算，這些內容不在本文討論的重點之中。git
在CORDIC被髮明以前，要對特殊函數求值，最天然的方法即是級數展開，例如利用泰勒展開來逼近目標函數，只要階數取得足夠大，就能夠無限逼近目標函數。級數展開在數學上是完美的，但運用到計算機時，咱們很快就會發現問題：級數展開本質是用多項式函數來近似目標函數，這其中包括大量複雜浮點運算，對於沒有硬件浮點運算單元的平臺，只能經過軟件浮點實現。github
CORDIC的出現解決了這個問題。該算法利用迭代逼近的方法，僅僅經過加/減和移位操做便可求解，極大的方便了計算機實現。算法
本文所作的工做：
函數
1、從CORDIC算法正向角度分析三角函數（sin, cos, tan）對任意角的求值；工具
2、從CORDIC算法逆向角度分析反三角函數（arcsin、arccos、arctan）對任意角的求值、向量模的求值；優化
3、分析定點運算的實現及軟件模型的創建spa
4、經過Verilog HDL設計硬件。設計
本文代碼倉庫詳見：https://github.com/cassuto/CORDIC-all-in-one-verilog 3d
2、核心思想
正如該算法的名字所說，CORDIC最初是爲一種用來進行座標軸旋轉的專用計算機開發的（其原型硬件於1959年成功應用於實時導航）。追根溯源，CORDIC的核心就是座標軸旋轉，算法最先的出發點是解決旋轉問題。它與函數求值究竟有什麼關係？早期研究者是如何想到這個算法的呢？這正是本文試圖說明的問題，固然並不完善，歡迎批評指正。調試
一、正向角度（已知旋轉角，求點旋轉後的座標）
假設如今咱們要將一個直角座標系中的點P(x0, y0)繞原點逆時針旋轉z0角度，則變換後的點P1(x1, y1)座標以下：
咱們發現這個轉換式中涉及三角函數，但如今假設機器還不能求任意三角函數值，那麼可否改進？
能夠考慮查表實現，但查表法要求數據必須是離散的，這樣旋轉角只能取有限個值。如何對任意的旋轉角求解？
簡單，利用迭代逼近，把目標角度分紅若干個小角度，每次迭代只旋轉一個特定角度並靠近目標，經過若干次迭代後，點就被旋轉到了近似的位置上。關鍵點是，經過特定的分割，使得每次旋轉的角度都是特定角，若是將些特定角的三角函數值固化到列表中，就能夠利用查表繞開三角函數的計算。對計算機而言雖然須要迭代屢次，但每次都是簡單運算，整體速度快於複雜的三角求值。通常該算法的結果只是近似值，但經過設置迭代次數就能夠控制精度，迭代次數越多，精度越高，並最終趨於穩定。
例如：將點P0旋轉z = 50°到P‘，求P’座標。
爲了便於製做三角函數表，首先假設：對全部問題來講，z的取值都是[0,90°]，所以假設每次旋轉的角度都是90°的n分之一。
第一步正向旋轉45°變換到P1，發現沒有達到目標角度；所以又繼續將P1正向旋轉22.5°變換到P2，發現超出了目標角度；因而繼續逆向旋轉11.25°變換到P3，仍然超出目標角度；再逆向旋轉5.625°，到達P4，這時若P4已經在設定精度內趨近答案，則該精度內能夠用P4近似表示P’座標。
這個算法基本上解決了旋轉變換式中三角函數求值的問題，但這並非咱們想要的：每次迭代都須要進行4次浮點乘法運算，旋轉角z的正餘弦都涉及浮點數，例如cos 22.5° ≈ 0.923879。
繼續觀察旋轉變換式，能夠變造成：
上式仍然涉及兩個三角函數，可是cos被放到了一邊，因而重點研究tan。將角度z微分爲n個小角度，若是像下面這樣特殊選取z，使tan z剛好只與2的冪有關，則本來複雜的乘tan zn運算就能夠經過右移n位實現（xn，yn都是整數，相關問題在後面討論），這對二進制計算機是很是天然的。有人指出0.618黃金分割是更優的，但並不適合二進制計算機實現。
結合上面的變換式子，並把cosz隱藏到係數P中，咱們就有了遞推關係（其中n表示當前迭代次數，n+1則表示下一次迭代。這裏只寫出了逆時針旋轉的狀況，順時針旋轉改變1/2n項的符號便可。）
因爲n是離散的，而且z = atan(1/2n)能夠提早計算，因此可用查表的方法快速得出z的值，這樣係數Pn也能夠經過查表求出，但係數Pn還是浮點數。
到這裏，原先算法的4次浮點運算，被成功減小到了只有2次，咱們離真正的CORDIC算法已經很近了。
但問題尚未結束，若是迭代n次的話，仍須要進行2n次浮點運算，有沒有優化的餘地呢？引發2n次浮點運算的緣由是每次都須要與係數Pn相乘，這樣作真的有必要嗎？
構造輔助三角形，能夠得出：
所以Pn的取值只與n有關，而與別的變量沒有關係。Pn徹底能夠在全部迭代都完成後單獨計算，在最後將結果乘P=ΠPn便可。
對於根據迭代較少的狀況，能夠將P的不一樣取值固化到表格中，經過查錶快速求出。
而對於迭代次數較多的狀況，能夠將P看做常數近似處理，這是由於隨着n的增長，P逐漸穩定下來：
經過數值統計，咱們能夠看出具體多少次迭代之後開始能夠將P看做常量。從圖中能夠看出從n=8開始比較合適。當n<3時P的取值變化較大，但事實上小於3次的迭代實用價值不大。
到這裏整個算法只需在最後進行一次浮點運算，而每次迭代都只涉及簡單的加/減和位移運算。
雖然有了這個快速旋轉變換算法，但咱們的最終目標並不在於此，而是求三角函數的值，如何經過旋轉求三角函數？
簡單。藉助單位圓這個工具，在直角座標系中將點(0, 1)繞原點旋轉至α角上，則旋轉後的點的橫縱座標對應的正是cosα、sinα的值，tanα也能夠利用比值求出，這正是利用CORDIC求三角函數的核心思路。下圖顯示了CORDIC算法的迭代過程。
上述算法可用僞代碼描述以下：
1 For n=0 to i-1 2 If (Z(n) >= 0) Then 3 X(n + 1) := X(n) – (Yn >> n) 4 Y(n + 1) := Y(n) + (Xn >> n) 5 6 Z(n + 1) := Z(n) - atan(1/2^n) 7 Else 8 X(n + 1) := X(n) + (Yn >> n) 9 Y(n + 1) := Y(n) – (Xn >> n) 10 11 Z(n + 1) := Z(n) + atan(1/2^n) 12 End if 13 End for
二、逆向角度（已知點旋轉後的點座標，求旋轉角）
從正向分析或者逆向分析，角度不一樣，解決的問題也不一樣。
仍然在單位圓內，假設已知正弦或餘弦中的一個值，如何求對應的反三角函數？這裏以求反正弦函數爲例說明，已知反正弦函數自變量爲S，咱們將點(1, 0)逆時針旋轉，若是轉到某個角度該點的縱座標恰等於S，說明這個角度就是反正弦函數的值；同理反餘弦函數也可用用相似的方法算出。
接下來繼續分析反正切的狀況，已知正切值，則能夠得出對應的單位圓上的點P0(x0, y0)。若是將點P0旋轉一個角度，使旋轉後的點縱座標變成0，那麼這個角度正是反正切函數的值。
以上即是利用CORDIC求反三角函數的全部思路。
如今留下的問題是到底旋轉多少角度才合適。CORDIC一樣以迭代逼近的方法解決這個問題。
以求反正切函數爲例：這裏只討論z在0到90°範圍內（即x0, y0都爲正）的狀況。咱們能夠先將原始點旋轉atan(-1/2)角度進行試探，旋轉完成後發現點縱座標仍然大於0，因而繼續旋轉atan(-1/4)、atan(-1/8)到達圖中草綠色終點，這個時候發現縱座標小於0，因而再反方向旋轉atan(-1/16)角度到達藍色終點，若此時點的縱座標在規定精度內接近0則迭代結束，能夠經過累加全部迭代中旋轉的角度增量求出α的具體數值，角度α正是atan的值。
經過第一小節的討論能夠看出，遞推式中係數P在幾何上的效果是縮放了旋轉半徑（點到座標原點的距離）。因爲求反正切函數過程當中只用到旋轉角，而旋轉角與旋轉半徑與無關，因此不須要考慮係數P。
上述算法可用僞代碼描述以下：
1 A := 0 2 For n=0 to i-1 3 If (Y(n) >= 0) Then 4 X(n + 1) := X(n) + (Yn >> n) 5 Y(n + 1) := Y(n) - (Xn >> n) 6 7 A := A + atan(1/2^n) 8 Else 9 X(n + 1) := X(n) - (Yn >> n) 10 Y(n + 1) := Y(n) + (Xn >> n) 11 12 A := A - atan(1/2^n) 13 End if 14 End for
理解了求反正切函數的思路後，咱們就能夠着手考慮反正弦和反餘弦函數的求值。比較可知，旋轉過程當中，對反正切函數要求旋轉點的縱座標趨於0，而對反正弦函數要求旋轉點縱座標趨於特定值S（函數自變量的取值）。能夠發現，前者所涉及的的旋轉是後者的一種特殊狀況，這樣只需簡單修改前者就能夠得出適用後者的旋轉算法。
1 limP := 0.607253 2 A := 0 3 B := S / limP 4 For n=0 to i-1 5 If (Y(n) >= B) Then 6 X(n + 1) := X(n) + (Yn >> n) 7 Y(n + 1) := Y(n) - (Xn >> n) 8 9 A := A + atan(1/2^n) 10 Else 11 X(n + 1) := X(n) - (Yn >> n) 12 Y(n + 1) := Y(n) + (Xn >> n) 13 14 A := A - atan(1/2^n) 15 End if 16 End for
須要解釋的是B := S / limP，這裏爲何須要將S放大1/limP倍呢？緣由很簡單，因爲在每次迭代時都忽略了係數Pn<1，這會致使Y(n)比未忽略係數時大出1/Pn倍。爲了使Y(n)與S可以進行比較，須要將S同比例放大1/Pn倍，這樣每次迭代又將引入浮點運算。直接近似Pn = limP雖然會引入細微的精度損失，但避免了大量浮點運算。
同理，也能夠得出旋轉點橫座標趨於特定值的算法。
1 limP := 0.607253 2 A := 0 3 B := S / limP 4 For n=0 to i-1 5 If (X(n) >= B) Then 6 X(n + 1) := X(n) - (Yn >> n) 7 Y(n + 1) := Y(n) + (Xn >> n) 8 9 A := A - atan(1/2^n) 10 Else 11 X(n + 1) := X(n) + (Yn >> n) 12 Y(n + 1) := Y(n) - (Xn >> n) 13 14 A := A + atan(1/2^n) 15 End if 16 End for
至此全部關於利用CORDIC求三角和反三角函數的分析就結束了。
借鑑上面求反正切角函數的思路，咱們也能夠看出向量求模的運算方法。直角座標系中，將目標向量V的起點平移到座標原點，終點用座標(x0, y0)表示。若是把這個點旋轉到x軸（或y軸）上，則旋轉後的點對應的橫座標（或縱座標）正是該向量的模長|V|。
該算法的核心與求反正切函數基本相同，這裏再也不贅述。
3、定點算法的軟件模型創建
若是能用整數表示浮點數，則能夠經過整數運算完成實數運算。在必定精度內，定點數與浮點數能夠經過比例放縮互相轉換。
首先肯定整數位寬，這裏以16位爲例，
一、量化角度：在直角範圍內，能夠將正整數的取值區間均分爲coeff = 2^14 / 90 = 182個單位，則一個單位對應角度的數量級爲10^-2。任意給定浮點角度Z°，其對應的定點角度爲floor(Z * coeff)；
二、量化座標(x, y)：將原始座標放大coeff倍並取整便可得出定點座標。但若座標取值較大則存在溢出風險，需另選系數。
上述兩點措施完成了算法輸入的定點量化；而對於算法的輸出，只需將定點數縮小coeff倍便可獲得最終的浮點結果。
在進行硬件設計以前，有必要先考慮軟件模型。由於軟件語言的抽象程度高於硬件描述語言，能夠更加緊湊地對算法進行描述，同時能快速地進行調試。
軟件模型主要涉及兩個核心算法：旋轉和反旋轉。
旋轉算法的C語言模型：
https://github.com/cassuto/CORDIC-all-in-one-verilog/blob/master/CORDIC-rotate-fixed-point.c
反旋轉算法：
https://github.com/cassuto/CORDIC-all-in-one-verilog/blob/master/CORDIC-anti-rotate-fixed-point.c
4、FPGA實現
這一小節主要解決的問題是硬件實現的相關細節，例如象限轉換，流水化處理等。採用流水線的目的在於提升時鐘頻率，例如在DDS（直接數字頻率合成）應用中，CORDIC算法能夠代替採樣表生成波形數據。爲了實現較高的合成頻率，流水化是有必要的。
關於Verilog實現有幾點須要說明的地方，首先利用verilog標準中規定的generate語句，能夠實現任意深度流水線的綜合；
其次在此以前對算法的全部討論都只涉及第Ⅰ象限，如今加入象限的處理。電路直接取相位輸入高2位判斷象限，並按函數值在4個象限的符號關係對輸出進行處理。
完整的Verilog模塊以下：該模塊經過端口phase_i輸入相位。經過sin_o，cos_o輸出三角函數值，經過err_o輸出迭代引發的相位偏差。從流水線爲空開始，需等待PIPE_DEPTH+2個時鐘週期才能得出結果；而流水線填滿後，對於後續相位輸入，模塊均可以在2個時鐘週期內更新輸出。
1 module cordic_dds # ( 2 parameter DW = 16, /* Data width */ 3 parameter PIPE_DEPTH = 14, /* Pipeline depth */ 4 parameter limP = 16'h4dba /* P = 0.607253 * 2^15 */ 5 ) 6 (/*AUTOARG*/ 7 // Outputs 8 sin_o, cos_o, err_o, 9 // Inputs 10 clk, phase_i 11 ); 12 13 input clk; 14 input [DW-1:0] phase_i; /* Phase */ 15 output [DW:0] sin_o, cos_o; /* Function value output */ 16 output [DW:0] err_o; /* Phase Error output */ 17 18 reg [DW:0] cos_r=0, sin_o_r=0; 19 reg [DW:0] x[PIPE_DEPTH:0]; 20 reg [DW:0] y[PIPE_DEPTH:0]; 21 reg [DW:0] z[PIPE_DEPTH:0]; 22 23 reg [DW:0] atan_rom[PIPE_DEPTH:0]; 24 25 reg [1:0] quadrant [PIPE_DEPTH:0]; 26 27 integer i; 28 initial begin 29 for(i=0; i<=PIPE_DEPTH; i=i+1) begin 30 x[i] = 0; y[i] = 0; z[i] = 0; 31 quadrant[i] = 2'b0; 32 end 33 end 34 35 initial begin 36 atan_rom[0] <= 8189; 37 atan_rom[1] <= 4834; 38 atan_rom[2] <= 2554; 39 atan_rom[3] <= 1296; 40 atan_rom[4] <= 650; 41 atan_rom[5] <= 325; 42 atan_rom[6] <= 162; 43 atan_rom[7] <= 81; 44 atan_rom[8] <= 40; 45 atan_rom[9] <= 20; 46 atan_rom[10] <= 10; 47 atan_rom[11] <= 5; 48 atan_rom[12] <= 2; 49 atan_rom[13] <= 1; 50 end 51 52 53 // ================= // 54 // Pipeline stages // 55 // ================= // 56 always @ (posedge clk) begin // stage 0 57 x[0] <= {1'b0, limP}; 58 y[0] <= 0; 59 z[0] <= {3'b0, phase_i[DW-1-2:0]}; // control the phase_i to the range[0-Pi/2] 60 end 61 62 always @ (posedge clk) begin // stage 1 63 x[1] <= x[0] - y[0]; 64 y[1] <= x[0] + y[0]; 65 z[1] <= z[0] - atan_rom[0]; // reversal 45deg 66 end 67 68 generate 69 genvar k; 70 for(k=1; k<PIPE_DEPTH; k=k+1) begin 71 always @ (posedge clk) begin 72 if (z[k][DW]) begin /* the diff is negative on clockwise */ 73 x[k+1] <= x[k] + {{k{y[k][DW]}},y[k][DW:k]}; /* >> k */ 74 y[k+1] <= y[k] - {{k{x[k][DW]}},x[k][DW:k]}; /* >> k */ 75 z[k+1] <= z[k] + atan_rom[k]; 76 end else begin 77 x[k+1] <= x[k] - {{k{y[k][DW]}},y[k][DW:k]}; 78 y[k+1] <= y[k] + {{k{x[k][DW]}},x[k][DW:k]}; 79 z[k+1] <= z[k] - atan_rom[k]; 80 end 81 end 82 end 83 endgenerate 84 85 // ================= // 86 // Count quadrant // 87 // ================= // 88 always @ (posedge clk) begin 89 quadrant[0] <= phase_i[DW-1:DW-2]; 90 end 91 generate 92 genvar j; 93 for(j=0; j<PIPE_DEPTH; j=j+1) begin 94 always @ (posedge clk) begin 95 quadrant[j+1] <= quadrant[j]; 96 end 97 end 98 endgenerate 99 100 // ================= // 101 // Adjust quadrant // 102 // ================= // 103 always @ (posedge clk) 104 case(quadrant[PIPE_DEPTH]) 105 2'b00: begin 106 cos_r <= x[PIPE_DEPTH]; /* cos */ 107 sin_o_r <= y[PIPE_DEPTH]; /* sin */ 108 end 109 2'b01: begin 110 cos_r <= ~(y[PIPE_DEPTH]) + 1'b1; /* -sin */ 111 sin_o_r <= x[PIPE_DEPTH]; /* cos */ 112 end 113 2'b10: begin 114 cos_r <= ~(x[PIPE_DEPTH]) + 1'b1; /* -cos */ 115 sin_o_r <= ~(y[PIPE_DEPTH]) + 1'b1; /* -sin */ 116 end 117 2'b11: begin 118 cos_r <= y[PIPE_DEPTH]; /* sin */ 119 sin_o_r <= ~(x[PIPE_DEPTH]) + 1'b1; /* -cos */ 120 end 121 endcase 122 123 assign cos_o = cos_r; 124 assign sin_o = sin_o_r; 125 assign err_o = z[PIPE_DEPTH]; 126 127 endmodule
經過仿真得出波形以下：
本文僅對CORDIC算法進行了一些淺顯的分析，基本覆蓋到了算法的核心思路和簡單應用，但從工業應用的角度出發，還有不少值得討論的問題未在文章中分析。因爲筆者時間有限，不能將本文作得很是全面，疏漏之處有待往後逐步完善，還望各位讀者海涵。