CORDIC算法原理及硬件實現（FPGA）

一、CORDIC算法

　　CORDIC（Coordinate Rotation DIgital Computer）是一種通過迭代實現快速平面旋轉的算法，通過變形擴展，它可以對多種超越函數求值，例如三角/反三角函數、雙曲函數等。

　　對超越函數求值，常見方法為用多項式近似，例如利用泰勒展開來逼近目標函數，只要階數取得足夠大，就可以無限逼近目標函數。當我們把這個方法運用到某些特殊機器時，很快會發現問題：泰勒展開包括大量復雜浮點運算，對於沒有硬件浮點運算單元的機器，只能通過軟件浮點實現。

　　CORDIC的出現解決了這個問題。該算法利用迭代逼近的方法，僅僅通過加/減和移位操作即可計算，極大的方便了機器實現。

 

本文的：

　　一、解釋CORDIC旋轉的原理。

　　二、從旋轉模式出發，解釋三角函數（sin, cos, tan）的計算；

　　三、從旋轉模式出發，解釋反三角函數（arcsin、arccos）的計算；

　　四、從向量模式出發，解釋向量模的計算；

　　五、介紹定點運算的實現及軟件模型的建立

　　六、通過Verilog HDL設計硬件。

 

本文代碼倉庫詳見：https://github.com/cassuto/CORDIC-all-in-one-verilog 

 

二、基本思想

　　正如該算法的名字所說，CORDIC最初是為一種用來進行坐標軸旋轉的專用計算機開發的（其原型硬件於1959年成功應用於實時導航）。追根溯源，CORDIC的核心就是坐標軸旋轉。其后又由不同的旋轉方式導出了各類函數的計算方法。

　　CORDIC包括旋轉模式和向量模式，下面逐一解釋。

 

1、旋轉模式

1.1 偽旋轉

　　假設現要將一個直角坐標系中的點P(x₀, y₀)繞原點逆時針旋轉z₀角度，則變換后的點P¹(x₁, y₁)坐標如下：

　　我們發現這個轉換式中涉及三角函數，但現在假設機器還不能求任意三角函數值，那么能否改進？

　　繼續觀察旋轉變換式，可以等價得到：

　　上式仍然涉及兩個三角函數，但是cos被放到了一邊，於是重點研究tan(z₀)。將角度z₀微分為n個小角度，每一步旋轉θ_n，通過若干步的旋轉最終逼近目標_。如果像下面這樣特殊選取θ_n，使tan θ_n恰好只與2的冪有關，則原本復雜的乘tan θ_n運算就可以通過二進制右移n位實現（x_n，y_n都是整數，相關問題在后面討論），這對二進制計算機非常自然。

$$ tan \theta_n = \frac{1}{2^n} $$

　　結合上面的變換公式，令cos θ_n=Pn，得到如下迭代格式：

　　分析計算量。θ_n = atan(1/2ⁿ)可查表得出。同理系數P_n也可以通過查表求出，但它們仍是浮點數。

　　到這里，算法原先的4次浮點運算，被成功減少到了2次，我們離CORDIC偽旋轉算法已經很近了。

　　但問題還沒有結束，如果迭代n次的話，仍需要進行2n次浮點運算。引起2n次浮點運算的原因是每次都需要與系數Pn相乘，這樣做真的有必要嗎？

　　構造輔助三角形，可以得出：

　　因此P_n的取值只與n有關。完全可以在所有迭代都完成后，再在最終結果上乘P=ΠP_n。

　　對於根據迭代較少的情況，可以將P_n的不同取值固化到表格中，通過查表求出。

　　但這仍需查表。進一步分析發現，對於迭代次數較多的情況，可以將P_n看作常數近似處理，因為隨着n的增加，P_n逐漸穩定下來：

　　即迭代多次時，可取P_n ≈ lim(n→∞)P_n，這成功去除了Pn查找表。不過引入了截斷誤差。誤差分析至關重要，本文有待完善。

　　通過數值方法，可以形象地看出Pn隨n增大的變化趨勢。

　　

　

　　同理可得順時針方向的旋轉算法。

對兩旋轉方向的公式進行歸納：定義符號函數S(n)控制旋轉方向，逆時針旋轉時S(n)取1；順時針旋轉時S(n)取-1。可以寫出更一般的偽旋轉公式：

（式1.1-1）

　　至此我們得到了一個完整的旋轉算法，一些文獻中也稱這個過程為“偽旋轉”。

 

1.2 迭代逼近

　　偽旋轉算法雖然避免了復雜的計算，但結果並不一定精確，原因在於每次旋轉的角度θ_n是嚴格受限的。例如目標角度Z₀取30°，第一次旋轉θ₀=atan(1)=45°，顯然已經遠遠超過了目標角，因此偽旋轉算法不一定收斂於精確值。

　　為了確保結果收斂於精確值，需要引入迭代逼近。思路很直觀：上例中“超過了目標角”，就以更小的角度向相反方向旋轉；若未“超過”，則繼續同向旋轉，如此反復，經過無窮次迭代后，點（X_n，Y_n）將最終收斂於目標點。

為了實現迭代逼近，首先引入角度累加器Zn，以確定相對於目標角度旋轉了多少。Z₀取目標角度，規定每步迭代滿足：

（式1.2-2）

　　則Zn正是目標角度與累計旋轉角度的差值，稱為相位誤差。如果已經轉過的角度大於目標角，則Zn<0，反之Zn>=0。

回到逼近的原理：如果轉過的角度大於目標角，則向相反方向旋轉，反之則繼續同向旋轉。發現，這正是Zn+1和S(n)的關系。

　　假設目標角Z₀>0表示逆時針旋轉，則S(n)與Zn+1可總結如下：

（式1.2-3）

　　聯立式1.1-1、1.2-2、1.2-3，便可得到完整的旋轉算法。可以證明，當n→∞時，Zn→0，點(X_n, Y_n)將收斂於目標點。

　　稱這種CORDIC模式為旋轉模式。

上述算法可用偽代碼描述如下：

For n=0 to i-1
    If (Z(n) >= 0) Then
        X(n + 1) := X(n) – (Yn >> n)
        Y(n + 1) := Y(n) + (Xn >> n)
        
        Z(n + 1) := Z(n) - atan(1/2^n)
    Else
        X(n + 1) := X(n) + (Yn >> n)
        Y(n + 1) := Y(n) – (Xn >> n)
        
        Z(n + 1) := Z(n) + atan(1/2^n)
    End if
End for
14 X(i) *= 0.607253 // limP

15 Y(i) *= 0.607253 // limP

 

1.3 通過旋轉模式計算三角函數

　　借助單位圓工具，在直角坐標系中，將點(1, 0)以原點為中心逆時針旋轉α角，旋轉后點的坐標(x',y')正是(cosα, sinα)。

下圖顯示了CORDIC算法計算三角函數的迭代過程。

 

 

 

 

2、向量模式

　　向量模式是為了計算向量模和反三角函數而引出的。

1.1 向量模的計算

　　直角坐標系中，設目標向量V(x₀, y₀)，向量模的定義如下：

　　向量模的計算：將點P₀(x₀, y₀)作為初始坐標。如果將點P₀旋轉某個角度，使縱坐標y'=0，則此時橫坐標x'即為向量V(x₀, y₀)的模。

　　與CORDIC旋轉模式類似，關鍵在於逼近的目標不同。取S(n)：

 

 

　　則旋轉過程中，Zn→0，向量將向Y_n=0的方向逼近，其過程如下（注意圖中向量模長比例並不精確）：

算法用偽代碼描述如下：

limP := 0.607253
For n=0 to i-1
    If (Y(n) >= 0) Then
        X(n + 1) := X(n) + (Yn >> n)
        Y(n + 1) := Y(n) - (Xn >> n)
 6     Else
        X(n + 1) := X(n) - (Yn >> n)
        Y(n + 1) := Y(n) + (Xn >> n)
 9     End if
End for
11 X(i) *= limP

 

1.2 反正弦函數的計算

　　已知正弦函數值S，在直角坐標系中，若當點(1, 0)逆時針旋轉某個角度A時縱坐標恰等於S，則A=asinS 即為反正弦函數的值。可見這只是正弦函數計算的逆過程。

　　仍為向量模式，只是逼近的目標發生了變化。取Zn = Y(n) - S，Zn為縱坐標相對於S的誤差。

　　則旋轉過程中，Zn向量將不斷向豎直分量Y_n=S的方向逼近。

 

算法用偽代碼描述如下：

A := 0
 2 S /= limP
For n=0 to i-1
    If (Y(n) >= S) Then
        X(n + 1) := X(n) + (Yn >> n)
        Y(n + 1) := Y(n) - (Xn >> n)
        
        A := A + atan(1/2^n)
    Else
        X(n + 1) := X(n) - (Yn >> n)
        Y(n + 1) := Y(n) + (Xn >> n)
        
        A := A - atan(1/2^n)
    End if
End for

 

同理可得反余弦函數的計算。

 

總結，CORDIC旋轉模式可完成sinA、cosA、asinA、acosA的計算。

 

三、定點算法的軟件模型建立

定點數

　　定點數可通過整數來表示。

　　1、量化角度：在直角范圍內，將b位無符號整數的取值區間均分為coeff₁ = 2^b / 90 個單位。

則對應任意給定的浮點角度Z°，其對應的整數可表示為：floor(Z * coeff₁)；

　　2、轉換坐標(x, y)：將原始坐標放大coeff₂倍並取整即可得出定點坐標。coeff₂應根據實際問題適當選取，盡量減小舍入誤差。

坐標用整數可表示為：（floor(X * coeff₂)、floor(Y * coeff₂)）。

　　3、對於算法的整數輸出，只需將整數除以相應coeff即可得到對應的浮點數。

　　浮點數到定點數的轉換過程將引入舍入誤差。

 

軟件模型

       軟件的抽象程度高於硬件，可以更加緊湊地對算法進行描述，同時能快速地進行調試。

　　通過軟件實現CORDIC算法：

　　　 旋轉模式的C語言模型：

　　　　　　https://github.com/cassuto/CORDIC-all-in-one-verilog/blob/master/CORDIC-rotate-fixed-point.c

　　　 向量模式的C語言模型：

　　　　　　https://github.com/cassuto/CORDIC-all-in-one-verilog/blob/master/CORDIC-anti-rotate-fixed-point.c

 　　完整軟件實現詳見：https://github.com/cassuto/CORDIC-all-in-one-verilog 

 

四、FPGA實現

　　CORDIC在硬件中的典型的應用包括DDS（直接數字頻率合成）等等。

　　本小節主要介紹硬件實現的相關細節，包括象限轉換，流水化處理。

 

象限轉換

　　根據三角函數的對稱性，只需實現一個象限，即可得出其余所有象限的情況。

　　由於相位取值范圍被象限四等分，可取相位輸入的高2位決定象限。得出象限后，即可根據三角函數值在4個象限的符號關系得到正確的輸出。

 

流水化處理

　　流水化提高了時鍾頻率，“以面積換取速度”。

　　從流水線為空開始，下游模塊需等待PIPE_DEPTH+2個時鍾周期才能從流水線得出結果；

　　而當流水線填滿后，對於相位輸入，流水線都可以在單時鍾周期內更新輸出。

 

Verilog設計

輸入：相位phase_i。

輸出：三角函數值sin_o，cos_o，相位誤差err_o

module cordic_dds # (
   parameter DW = 16,               /* Data width */
   parameter PIPE_DEPTH = 14,       /* Pipeline depth */
   parameter limP = 16'h4dba        /* P = 0.607253 * 2^15 */
)
(/*AUTOARG*/
   // Outputs
   sin_o, cos_o, err_o,
   // Inputs
   clk, phase_i
);

   input             clk;
   input [DW-1:0]    phase_i;       /* Phase */
   output [DW:0]     sin_o, cos_o;  /* Function value output */
   output [DW:0]     err_o;         /* Phase Error output */

   reg [DW:0] cos_r=0, sin_o_r=0;
   reg [DW:0] x[PIPE_DEPTH:0];
   reg [DW:0] y[PIPE_DEPTH:0];
   reg [DW:0] z[PIPE_DEPTH:0];

   reg [DW:0] atan_rom[PIPE_DEPTH:0];
   
   reg [1:0] quadrant [PIPE_DEPTH:0];

   integer i;
   initial begin
      for(i=0; i<=PIPE_DEPTH; i=i+1) begin
         x[i] = 0; y[i] = 0; z[i] = 0;
         quadrant[i] = 2'b0;
      end
   end
   
   initial begin
      atan_rom[0] <= 8189;
      atan_rom[1] <= 4834;
      atan_rom[2] <= 2554;
      atan_rom[3] <= 1296;
      atan_rom[4] <= 650;
      atan_rom[5] <= 325;
      atan_rom[6] <= 162;
      atan_rom[7] <= 81;
      atan_rom[8] <= 40;
      atan_rom[9] <= 20;
      atan_rom[10] <= 10;
      atan_rom[11] <= 5;
      atan_rom[12] <= 2;
      atan_rom[13] <= 1;
   end
   
   
   // ================= //
   // Pipeline stages   //
   // ================= //
   always @ (posedge clk) begin // stage 0
      x[0] <= {1'b0, limP};
      y[0] <= 0;
      z[0] <= {3'b0, phase_i[DW-1-2:0]}; // control the phase_i to the range[0-Pi/2]
   end

   always @ (posedge clk) begin // stage 1
      x[1] <= x[0] - y[0];
      y[1] <= x[0] + y[0];
      z[1] <= z[0] - atan_rom[0]; // reversal 45deg
   end

   generate
      genvar k;
      for(k=1; k<PIPE_DEPTH; k=k+1) begin
         always @ (posedge clk) begin
            if (z[k][DW]) begin /* the diff is negative on clockwise */
               x[k+1] <= x[k] + {{k{y[k][DW]}},y[k][DW:k]}; /* >> k */
               y[k+1] <= y[k] - {{k{x[k][DW]}},x[k][DW:k]}; /* >> k */
               z[k+1] <= z[k] + atan_rom[k];
            end else begin
               x[k+1] <= x[k] - {{k{y[k][DW]}},y[k][DW:k]};
               y[k+1] <= y[k] + {{k{x[k][DW]}},x[k][DW:k]};
               z[k+1] <= z[k] - atan_rom[k];
            end
         end
      end
   endgenerate

   // ================= //
   // Count quadrant    //
   // ================= //
   always @ (posedge clk) begin
      quadrant[0] <= phase_i[DW-1:DW-2];
   end
   generate
      genvar j;
      for(j=0; j<PIPE_DEPTH; j=j+1) begin
         always @ (posedge clk) begin
            quadrant[j+1] <= quadrant[j];
         end
      end
   endgenerate

   // ================= //
   // Adjust quadrant   //
   // ================= //
   always @ (posedge clk)
      case(quadrant[PIPE_DEPTH])
         2'b00: begin
            cos_r <= x[PIPE_DEPTH]; /* cos */
            sin_o_r <= y[PIPE_DEPTH]; /* sin */
         end
         2'b01: begin
            cos_r <= ~(y[PIPE_DEPTH]) + 1'b1; /* -sin */
            sin_o_r <= x[PIPE_DEPTH]; /* cos */
         end
         2'b10: begin
            cos_r <= ~(x[PIPE_DEPTH]) + 1'b1; /* -cos */
            sin_o_r <= ~(y[PIPE_DEPTH]) + 1'b1; /* -sin */
         end
         2'b11: begin
            cos_r <= y[PIPE_DEPTH]; /* sin */
            sin_o_r <= ~(x[PIPE_DEPTH]) + 1'b1; /* -cos */
         end
      endcase

   assign cos_o = cos_r;
   assign sin_o = sin_o_r;
   assign err_o = z[PIPE_DEPTH];

endmodule

 

 

　　通過仿真得出波形如下：