怎樣使用CORDIC算法求解角度正余弦呢？

1、算法簡(jiǎn)介

CORDIC（Coordinate Rotation Digital Computer）算法即坐標(biāo)旋轉(zhuǎn)數(shù)字計(jì)算方法，是J.D.Volder1于1959年首次提出，主要用于三角函數(shù)、雙曲線、指數(shù)、對(duì)數(shù)的計(jì)算。該算法通過(guò)基本的加和移位運(yùn)算代替乘法運(yùn)算，使得矢量的旋轉(zhuǎn)和定向的計(jì)算不再需要三角函數(shù)、乘法、開(kāi)方、反三角、指數(shù)等函數(shù)，計(jì)算向量長(zhǎng)度并能把直角坐標(biāo)系轉(zhuǎn)換為極坐標(biāo)系。因?yàn)镃ordic 算法只用了移位和加法，很容易用純硬件來(lái)實(shí)現(xiàn)，非常適合FPGA實(shí)現(xiàn)。

CORDIC算法是天平稱重思想在數(shù)值運(yùn)算領(lǐng)域的杰出范例。核心的思想是把非線性的問(wèn)題變成了線性的迭代問(wèn)題【4】。

CORDIC算法完成坐標(biāo)或向量的平面旋轉(zhuǎn)（下圖以逆時(shí)針旋轉(zhuǎn)為例）。

旋轉(zhuǎn)后，可得如下向量：

旋轉(zhuǎn)的角度θ經(jīng)過(guò)多次旋轉(zhuǎn)得到的（步步逼近，接近二分查找法），每次旋轉(zhuǎn)一小角度。單步旋轉(zhuǎn)定義如下公式：

公式（2）提取cosθ，可修改為：

修改后的公式把乘法次數(shù)從4次改為3次，剩下的乘法運(yùn)算可以通過(guò)選擇每次旋轉(zhuǎn)的角度去除，將每一步的正切值選為2的指數(shù)（二分查找法），除以2的指數(shù)可以通過(guò)右移操作完成（verilog）。

每次旋轉(zhuǎn)的角度可以表示為：

所有迭代角度累加值等于最終需要的旋轉(zhuǎn)角度θ：

這里Sn為1或者-1，根據(jù)旋轉(zhuǎn)方向確定（后面有確定方法，公式（15）），順時(shí)針為-1，逆時(shí)針為1。

可以得到如下公式：

結(jié)合公式（3）和（7），得到公式（8）：

到這里，除了余弦值這個(gè)系數(shù)，算法只要通過(guò)簡(jiǎn)單的移位和加法操作完成。而這個(gè)系數(shù)可以通過(guò)預(yù)先計(jì)算最終值消掉。首先重新重寫(xiě)這個(gè)系數(shù)如下：

第二步計(jì)算所有的余弦值并相乘，這個(gè)值K稱為增益系數(shù)。

由于K值是常量，我們可以先忽略它。

到這里我們發(fā)現(xiàn)，算法只剩下移位和加減法，這就非常適合硬件實(shí)現(xiàn)了，為硬件快速計(jì)算三角函數(shù)提供了一種新的算法。在進(jìn)行迭代運(yùn)算時(shí)，需要引入一個(gè)新的變量Z，表示需要旋轉(zhuǎn)的角度θ中還沒(méi)有旋轉(zhuǎn)的角度。

這里，我們可以把前面提到確定旋轉(zhuǎn)方向的方法介紹了，就是通過(guò)這個(gè)變量Z的符號(hào)確定。

通過(guò)公式（5）和（15），將未旋轉(zhuǎn)的角度變?yōu)?。

一個(gè)類編程風(fēng)格的結(jié)構(gòu)如下，反正切值是預(yù)先計(jì)算好的。

1.1 旋轉(zhuǎn)模式

旋轉(zhuǎn)模式下，CORDIC算法驅(qū)動(dòng)Z變?yōu)?，結(jié)合公式（13）和（16），算法的核心計(jì)算如下：

一種特殊情況是，另初始值如下：

因此，旋轉(zhuǎn)模式下CORDIC算法可以計(jì)算一個(gè)輸入角度的正弦值和余弦值。

1.2 向量模式

向量模式下，有兩種特例：

因此，向量模式下CORDIC算法可以用來(lái)計(jì)算輸入向量的模和反正切，也能開(kāi)方計(jì)算，并可以將直角坐標(biāo)轉(zhuǎn)換為極坐標(biāo)。

2、硬件算法實(shí)現(xiàn)

根據(jù)【5】，可以看到CORDIC迭代算法的一種直接實(shí)現(xiàn)方式是反饋結(jié)構(gòu)，此結(jié)構(gòu)只設(shè)計(jì)一級(jí)CORDIC運(yùn)算迭代單元、然后在系統(tǒng)時(shí)鐘的驅(qū)動(dòng)下，將本級(jí)的輸出作為本級(jí)的輸入，通過(guò)同一級(jí)迭代完成運(yùn)算。這種方法硬件開(kāi)銷小、但控制相對(duì)復(fù)雜。

所以根據(jù)【1】、【2】，使用流水線結(jié)構(gòu)實(shí)現(xiàn)了CORDIC迭代算法、求取了角度的正弦和余弦值。

下面分段介紹下各部分代碼：

首先是角度的表示，進(jìn)行了宏定義，360讀用16位二進(jìn)制表示2^16，每一度為2^16/360。

//360°--2^16,phase_in = 16bits (input [15:0] phase_in)
//1°--2^16/360
`define rot0  16'h2000    //45
`define rot1  16'h12e4    //26.5651
`define rot2  16'h09fb    //14.0362
`define rot3  16'h0511    //7.1250
`define rot4  16'h028b    //3.5763
`define rot5  16'h0145    //1.7899
`define rot6  16'h00a3    //0.8952
`define rot7  16'h0051    //0.4476
`define rot8  16'h0028    //0.2238
`define rot9  16'h0014    //0.1119
`define rot10 16'h000a    //0.0560
`define rot11 16'h0005    //0.0280
`define rot12 16'h0003    //0.0140
`define rot13 16'h0002    //0.0070
`define rot14 16'h0001    //0.0035
`define rot15 16'h0000    //0.0018

然后是流水線級(jí)數(shù)定義、增益放大倍數(shù)以及中間結(jié)果位寬定義。流水線級(jí)數(shù)16，為了滿足精度要求，有文獻(xiàn)指出流水線級(jí)數(shù)必須大于等于角度位寬16（針對(duì)正弦余弦計(jì)算的CORDIC算法優(yōu)化及其FPGA實(shí)現(xiàn)）。增益放大2^16，為了避免溢出狀況中間結(jié)果（x，y，z）定義為17為，最高位作為符號(hào)位判斷，1為負(fù)數(shù)，0為正數(shù)。

module cordic
(
    input clk,
    
    input [15:0] phase_in,
    output reg signed [16:0] eps,
    output reg signed [16:0] sin,
    output reg signed [16:0] cos
);
parameter PIPELINE = 16;

parameter K = 16'h9b74;
//gian k=0.607253*2^16,9b74,n means the number pipeline
//pipeline 16-level    //maybe overflow,matlab result not overflow
//MSB is signed bit,transform the sin and cos according to phase_in[15:14]
reg signed [16:0] x0=0,y0=0,z0=0;
reg signed [16:0] x1=0,y1=0,z1=0;
reg signed [16:0] x2=0,y2=0,z2=0;
reg signed [16:0] x3=0,y3=0,z3=0;
reg signed [16:0] x4=0,y4=0,z4=0;
reg signed [16:0] x5=0,y5=0,z5=0;
reg signed [16:0] x6=0,y6=0,z6=0;
reg signed [16:0] x7=0,y7=0,z7=0;
reg signed [16:0] x8=0,y8=0,z8=0;
reg signed [16:0] x9=0,y9=0,z9=0;
reg signed [16:0] x10=0,y10=0,z10=0;
reg signed [16:0] x11=0,y11=0,z11=0;
reg signed [16:0] x12=0,y12=0,z12=0;
reg signed [16:0] x13=0,y13=0,z13=0;
reg signed [16:0] x14=0,y14=0,z14=0;
reg signed [16:0] x15=0,y15=0,z15=0;
reg signed [16:0] x16=0,y16=0,z16=0;

還需要定義memory型寄存器數(shù)組并初始化為0，用于寄存輸入角度高2位的值。

reg [1:0] quadrant [PIPELINE:0];
integer i;
initial
begin
    for(i=0;i<=PIPELINE;i=i+1)
    quadrant[i] = 2'b0;
end

接著，是對(duì)輸入角度象限處理，將角度都轉(zhuǎn)換到第一象限，方便處理。輸入角度值最高兩位賦值0，即轉(zhuǎn)移到第一象限[0°,90°]。此外，完成x0，y0和z0的初始化，并增加一位符號(hào)位。

always @ (posedge clk)//stage 0,not pipeline
begin
    x0 <= {1'b0,K}; //add one signed bit,0 means positive
    y0 <= 17'b0;
    z0 <= {3'b0,phase_in[13:0]};//control the phase_in to the range[0-Pi/2]
end

接下來(lái)根據(jù)剩余待旋轉(zhuǎn)角度z的符號(hào)位進(jìn)行16次迭代處理，即完成16級(jí)流水線處理。

View Code

其中使用到了算數(shù)右移（>>>）運(yùn)算、所以在之前聲明時(shí)將相應(yīng)的reg/wire均聲明為signed類型。這一點(diǎn)在【1】的最后也有說(shuō)明。

需要注意的是這里的算數(shù)移位運(yùn)算（這一運(yùn)算的詳細(xì)過(guò)程在【3】中進(jìn)行了說(shuō)明），與之區(qū)分的是邏輯移位運(yùn)算。

二者規(guī)則為：

邏輯左移和右移，空出的位均補(bǔ)零。

算數(shù)左移與邏輯左移相同，都在低位補(bǔ)零；算數(shù)右移、移出的高位比特使用符號(hào)位填充（0正1負(fù)）

舉例說(shuō)明，對(duì)8'b_1011_0111進(jìn)行邏輯、算數(shù)移位的結(jié)果如下圖所示：

比如這里的原數(shù)值是8'b10110111、為負(fù)數(shù)（補(bǔ)碼形式）、換算成十進(jìn)制為-73.

算數(shù)右移一位之后的結(jié)果是8'b11011011、由補(bǔ)碼換算成原碼再換算為十進(jìn)制為-37.

由于進(jìn)行了象限的轉(zhuǎn)換，最終流水結(jié)果需要根據(jù)象限進(jìn)行轉(zhuǎn)換為正確的值。這里寄存17次高2位角度輸入值，配合流水線結(jié)果用于象限判斷，并完成轉(zhuǎn)換。

//according to the pipeline,register phase_in[15:14]
always @ (posedge clk)
begin
quadrant[0] <= phase_in[15:14];
quadrant[1] <= quadrant[0];
quadrant[2] <= quadrant[1];
quadrant[3] <= quadrant[2];
quadrant[4] <= quadrant[3];
quadrant[5] <= quadrant[4];
quadrant[6] <= quadrant[5];
quadrant[7] <= quadrant[6];
quadrant[8] <= quadrant[7];
quadrant[9] <= quadrant[8];
quadrant[10] <= quadrant[9];
quadrant[11] <= quadrant[10];
quadrant[12] <= quadrant[11];
quadrant[13] <= quadrant[12];
quadrant[14] <= quadrant[13];
quadrant[15] <= quadrant[14];
quadrant[16] <= quadrant[15];
end

最后，根據(jù)寄存的高2位角度輸入值，利用三角函數(shù)關(guān)系，得出最后的結(jié)果，其中負(fù)數(shù)進(jìn)行了補(bǔ)碼操作。

//alter register, according to quadrant[16] to transform the result to the right result
always @ (posedge clk)
    eps <= z16;

always @ (posedge clk) begin
case(quadrant[16]) //or 15
2'b00:begin //if the phase is in first quadrant,the sin(X)=sin(A),cos(X)=cos(A)
        cos <= x16;
        sin <= y16;
        end
2'b01:begin //if the phase is in second quadrant,the sin(X)=sin(A+90)=cosA,cos(X)=cos(A+90)=-sinA
        cos <= ~(y16) + 1'b1;//-sin
        sin <= x16;//cos
        end
2'b10:begin //if the phase is in third quadrant,the sin(X)=sin(A+180)=-sinA,cos(X)=cos(A+180)=-cosA
        cos <= ~(x16) + 1'b1;//-cos
        sin <= ~(y16) + 1'b1;//-sin
        end
2'b11:begin //if the phase is in forth quadrant,the sin(X)=sin(A+270)=-cosA,cos(X)=cos(A+270)=sinA
        cos <= y16;//sin
        sin <= ~(x16) + 1'b1;//-cos
        end
endcase
end

完整代碼：

Whole Code

testbench測(cè)試代碼：

Testbench

3、Modelsim仿真結(jié)果

仿真結(jié)果的補(bǔ)充說(shuō)明：

（1）程序全程未使用復(fù)位信號(hào)，testbench中第一個(gè)計(jì)算的角度為16'h2000也就是45度，如果以圖示時(shí)刻為0時(shí)刻、仿真結(jié)果對(duì)應(yīng)的波形即分別為sin(x+π/4)和cos(x+π/4)的波形。作為參考，0.5*√2*65535≈46340.

（2）關(guān)于運(yùn)算過(guò)程中的位數(shù)溢出

根據(jù)仿真結(jié)果，本測(cè)試?yán)?，x4出現(xiàn)過(guò)16位位數(shù)溢出。

（3）關(guān)于流水線設(shè)計(jì)的理解

前文提到過(guò)實(shí)現(xiàn)CORDIC迭代算法時(shí)可以使用反饋結(jié)構(gòu)（只使用一級(jí)）、也可以使用流水線結(jié)構(gòu)（多級(jí)），如果任務(wù)是只單獨(dú)計(jì)算一個(gè)角度的正弦或者余弦值，那么所需要的迭代次數(shù)或者說(shuō)消耗的時(shí)鐘周期數(shù)量其實(shí)是相同的，本設(shè)計(jì)中為16個(gè)時(shí)鐘。

流水線結(jié)構(gòu)的威力是在連續(xù)、源源不斷地計(jì)算一組多個(gè)角度的正余弦值的時(shí)候才展現(xiàn)出來(lái)，當(dāng)初始流水線被填滿之后，每經(jīng)過(guò)一個(gè)時(shí)鐘周期、都會(huì)在輸出上獲得一個(gè)更新的角度的正余弦結(jié)果值，上圖仿真結(jié)果圖中黃色cursor左側(cè)的時(shí)間段內(nèi)、流水線即被逐步填滿。

換句話說(shuō)，如果現(xiàn)在的任務(wù)是要計(jì)算n個(gè)角度的正余弦值、計(jì)算一個(gè)角度需要的迭代次數(shù)為x，反饋結(jié)構(gòu)需要的時(shí)長(zhǎng)為（n*x）個(gè)時(shí)鐘周期，流水線結(jié)構(gòu)只需要（n+x-1）個(gè)時(shí)鐘周期。

審核編輯：劉清