龙空技术网

使用CORDIC算法求解角度正余弦及Verilog实现

梦Love天下 128

前言:

当前姐妹们对“传统角度算法”大体比较讲究,各位老铁们都需要分析一些“传统角度算法”的相关知识。那么小编在网上汇集了一些有关“传统角度算法””的相关资讯,希望我们能喜欢,朋友们一起来学习一下吧!

#记录我的2023#使用CORDIC算法求解角度正余弦及Verilog实现

本文是用于记录在了解和学习CORDIC算法期间的收获,以供日后自己及他人参考;并且附上了使用Verilog实现CORDIC算法求解角度的正弦和余弦的代码、简单的testbench测试代码、以及在Modelsim下的仿真结果

本文主要参考了:

【1】 (cordic算法的verilog实现及modelsim仿真)

【2】(CORDIC算法--流水线结构)

【3】([文档].艾米电子 - 算术运算电路,Verilog)

【4】(如何理解CORDIC算法)

【5】(针对正弦余弦计算的CORDIC算法优化及其FPGA实现)

感谢!

-------------------------------------------------------------------------------------------------------------

1、算法简介

CORDIC(Coordinate Rotation Digital Computer)算法即坐标旋转数字计算方法,是J.D.Volder1于1959年首次提出,主要用于三角函数、双曲线、指数、对数的计算。该算法通过基本的加和移位运算代替乘法运算,使得矢量的旋转和定向的计算不再需要三角函数、乘法、开方、反三角、指数等函数,计算向量长度并能把直角坐标系转换为极坐标系。因为Cordic 算法只用了移位和加法,很容易用纯硬件来实现,非常适合FPGA实现。

CORDIC算法是天平称重思想在数值运算领域的杰出范例。核心的思想是把非线性的问题变成了线性的迭代问题【4】。

CORDIC算法完成坐标或向量的平面旋转(下图以逆时针旋转为例)。

图片

旋转后,可得如下向量:

图片

旋转的角度θ经过多次旋转得到的(步步逼近,接近二分查找法),每次旋转一小角度。单步旋转定义如下公式:

图片

公式(2)提取cosθ,可修改为:

图片

修改后的公式把乘法次数从4次改为3次,剩下的乘法运算可以通过选择每次旋转的角度去除,将每一步的正切值选为2的指数(二分查找法),除以2的指数可以通过右移操作完成(verilog)。

每次旋转的角度可以表示为:

image

所有迭代角度累加值等于最终需要的旋转角度θ:

image

这里Sn为1或者-1,根据旋转方向确定(后面有确定方法,公式(15)),顺时针为-1,逆时针为1。

image

可以得到如下公式:

image

结合公式(3)和(7),得到公式(8):

image

到这里,除了余弦值这个系数,算法只要通过简单的移位和加法操作完成。而这个系数可以通过预先计算最终值消掉。首先重新重写这个系数如下:

图片

第二步计算所有的余弦值并相乘,这个值K称为增益系数。

图片

由于K值是常量,我们可以先忽略它。

图片

图片

到这里我们发现,算法只剩下移位和加减法,这就非常适合硬件实现了,为硬件快速计算三角函数提供了一种新的算法。在进行迭代运算时,需要引入一个新的变量Z,表示需要旋转的角度θ中还没有旋转的角度。

图片

这里,我们可以把前面提到确定旋转方向的方法介绍了,就是通过这个变量Z的符号确定。

图片

图片

通过公式(5)和(15),将未旋转的角度变为0。

一个类编程风格的结构如下,反正切值是预先计算好的。

图片

1.1 旋转模式

旋转模式下,CORDIC算法驱动Z变为0,结合公式(13)和(16),算法的核心计算如下:

图片

一种特殊情况是,另初始值如下:

图片

因此,旋转模式下CORDIC算法可以计算一个输入角度的正弦值和余弦值。

1.2 向量模式

向量模式下,有两种特例:

图片

因此,向量模式下CORDIC算法可以用来计算输入向量的模和反正切,也能开方计算,并可以将直角坐标转换为极坐标。

算法介绍:,

2、硬件算法实现

根据【5】,可以看到CORDIC迭代算法的一种直接实现方式是反馈结构,此结构只设计一级CORDIC运算迭代单元、然后在系统时钟的驱动下,将本级的输出作为本级的输入,通过同一级迭代完成运算。这种方法硬件开销小、但控制相对复杂。

所以根据【1】、【2】,使用流水线结构实现了CORDIC迭代算法、求取了角度的正弦和余弦值。

下面分段介绍下各部分代码:

首先是角度的表示,进行了宏定义,360读用16位二进制表示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

图片

然后是流水线级数定义、增益放大倍数以及中间结果位宽定义。流水线级数16,为了满足精度要求,有文献指出流水线级数必须大于等于角度位宽16(针对正弦余弦计算的CORDIC算法优化及其FPGA实现)。增益放大2^16,为了避免溢出状况中间结果(x,y,z)定义为17为,最高位作为符号位判断,1为负数,0为正数。

图片

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型寄存器数组并初始化为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

图片

接着,是对输入角度象限处理,将角度都转换到第一象限,方便处理。输入角度值最高两位赋值0,即转移到第一象限[0°,90°]。此外,完成x0,y0和z0的初始化,并增加一位符号位。

图片

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

图片

接下来根据剩余待旋转角度z的符号位进行16次迭代处理,即完成16级流水线处理。

图片View Code

其中使用到了算数右移(>>>)运算、所以在之前声明时将相应的reg/wire均声明为signed类型。这一点在【1】的最后也有说明。

需要注意的是这里的算数移位运算(这一运算的详细过程在【3】中进行了说明),与之区分的是逻辑移位运算。

二者规则为:

逻辑左移和右移,空出的位均补零。

算数左移与逻辑左移相同,都在低位补零;算数右移、移出的高位比特使用符号位填充(0正1负)

举例说明,对8'b_1011_0111进行逻辑、算数移位的结果如下图所示:

图片

比如这里的原数值是8'b10110111、为负数(补码形式)、换算成十进制为-73.

算数右移一位之后的结果是8'b11011011、由补码换算成原码再换算为十进制为-37.

由于进行了象限的转换,最终流水结果需要根据象限进行转换为正确的值。这里寄存17次高2位角度输入值,配合流水线结果用于象限判断,并完成转换。

图片

//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

图片

最后,根据寄存的高2位角度输入值,利用三角函数关系,得出最后的结果,其中负数进行了补码操作。

图片

//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测试代码:

图片Testbench

3、Modelsim仿真结果

图片

仿真结果的补充说明:

(1)程序全程未使用复位信号,testbench中第一个计算的角度为16'h2000也就是45度,如果以图示时刻为0时刻、仿真结果对应的波形即分别为sin(x+π/4)和cos(x+π/4)的波形。作为参考,0.5*√2*65535≈46340.

(2)关于运算过程中的位数溢出

根据仿真结果,本测试例下,x4出现过16位位数溢出。

(3)关于流水线设计的理解

前文提到过实现CORDIC迭代算法时可以使用反馈结构(只使用一级)、也可以使用流水线结构(多级),如果任务是只单独计算一个角度的正弦或者余弦值,那么所需要的迭代次数或者说消耗的时钟周期数量其实是相同的,本设计中为16个时钟。

流水线结构的威力是在连续、源源不断地计算一组多个角度的正余弦值的时候才展现出来,当初始流水线被填满

标签: #传统角度算法