1.CORDIC算法简介和基本原理
关于CORDIC算法的介绍和原理部分我就不重复介绍了，读者可参阅下面的连接，也可以自行查找其它介绍，或者相关文献：
http://mp.ofweek.com/it/a656714611187上面的连接中介绍了CORDIC算法的基本原理。实际上最初提出的CODIC算法是基于圆坐标的旋转，后来的学者们根据CORDIC算法的基本原理发展了很多其它的变种形式用来计算其它函数，如上面链接中所提及的线性旋转模式、双曲线旋
转模式等，也有相应的优化形式，其优化出发点主要是针对迭代收敛速度的、减少所需寄存器的等，这里不过多介绍。
尽管CORDIC算法也可以实现如乘除法等运算，但它的限制较大，实际应用CORDIC算法的场景主要是计算三角函数的圆坐标旋转模式，包含旋转模式和向量模式两种模式，分别用来计算cos、sin值和actan值。
2.CORDIC算法的MATLAB实现
2.0获得参考基
n=16;
tanQ=zeros(n,1);
tanhQ=zeros(n,1);
thetad=zeros(n,1);
thetah=zeros(n,1);
cosQ=zeros(n,1);
coshQ=zeros(n,1);
K_m1=zeros(n,1);
K_h=zeros(n,1);
max=0;
max2=0;
for i=1:n
tanQ(i)=2^(-i+1);
tanhQ(i)=2^(-i);
thetad(i)=atand(tanQ(i));% 以度为单位
cosQ(i)=cosd(thetad(i));% 以度为单位
thetah(i)=atanh(tanhQ(i));
thetahd=rad2deg(thetah);% 以度为单位
coshQ(i)=cosh(thetah(i));
% 计算伸缩系数K_n的倒数
if(i==1)
K_m1(i)=1*cosQ(i); % 第一个值为1与当前cos的乘积
else
K_m1(i)=K_m1(i-1)*cosQ(i);% 第二个以后的值为前一值与当前cos的乘积，
end
% 计算双曲线旋转的伸缩因子
if(i==1)
K_h(i)=1*sqrt(1-2^(-2*(i))); % 第一个值为1与当前cos的乘积
else
K_h(i)=K_h(i-1)*sqrt(1-2^(-2*(i)));% 第二个以后的值为前一值与当前cos的乘积，
end
max=max+thetad(i);% 计算可计算的角度边界
end
K=1./K_m1;
theta=deg2rad(thetad);
output=table(thetad,theta,tanQ,cosQ,K_m1,K,thetahd,thetah,tanhQ,coshQ,K_h)
writetable(output,'CORDIC_table.xlsx') % 输出结果到Excel表格以备查用
max % 输出圆坐标旋转模式可计算的最大值角度
上面程序在MATLAB的.mlx文件中输出结果如上图。
这里的思想就是通过一些列的tan=2^(-i)值调用反三角函数算出相对应的角度（本文中MATLAB程序采用单位为弧度，单位为度的结果也有输出如有需要可自行调取），运行上面的代码可以得到一个包含所需参考角度的Excel表格，下面计算中可以直接读取和调用表格中的数据。
2.1旋转模式（Rotation Mode）
n=16;% 迭代次数
x=zeros(n+1,1);
y=zeros(n+1,1);
z=zeros(n+1,1);
theta_r=zeros(n+1,1);
d=zeros(n+1,1);
x(1)=K_m1(16);% x初始值为1/K_n
y(1)=0; % y初始值为0
Q=pi/4;
z(1)=Q;% z初始值为待求角度
for i=1:n
if z(i)>=0
d(i)=1;
else
d(i)=-1;
end
x(i+1)=x(i)-d(i)*y(i)*2^(-i+1);
y(i+1)=y(i)+d(i)*x(i)*2^(-i+1);
z(i+1)=z(i)-d(i)*output.theta(i);
theta_r(i+1)=z(i+1);% theta result in rad during process
end
m=10;% 结果保留m位有效数字
cosQ=vpa(x(17),m);% 计算得出的cos值
ecosQ=vpa(cos(Q)-x(17),m);
sinQ=vpa(y(17),m);%计算得出的sin值
esinQ=vpa(sin(Q)-y(17),m);
c=vpa(z(17),m);
theta_d=rad2deg(theta_r);
output_c=table(x,y,z,d,theta_r,theta_d)
上面程序在MATLAB的.mlx文件中输出结果如上图。
程序中给定了角度为pi/4它的正弦和余弦都为，程序仅通过16次迭代即得到了相对精确的结果。
2.2向量模式（Vector Mode）
n=16;% 迭代次数
x=zeros(n+1,1);
y=zeros(n+1,1);
z=zeros(n+1,1);
d=zeros(n+1,1);
theta_r=zeros(n+1,1);
x(1)=1;% x初始值为1/K_n
y(1)=1; % y初始值为0
z(1)=0;% z初始值为待求角度
for i=1:n
if y(i)>=0
d(i)=-1;
else
d(i)=1;
end
x(i+1)=x(i)-d(i)*y(i)*2^(-i+1);
y(i+1)=y(i)+d(i)*x(i)*2^(-i+1);
z(i+1)=z(i)-d(i)*output.theta(i);
theta_r(i+1)=z(i+1);% theta result in rad during process
end
m=10;% 结果保留m位有效数字
k_n=1.64676025786545e+00;
x_n=vpa(x(17),m);
y_n=vpa(y(17),m);
y_0=sqrt((x_n/k_n)^2-x(1)^2);
theta_d=rad2deg(theta_r);
theta_f=vpa(rad2deg(z(17)),m)% 计算得出的角度值
output_v=table(x,y,z,d,theta_r,theta_d)
上面程序在MATLAB的.mlx文件中输出结果如上图。
程序中给定了，结果输出角度为45.00088972，十分接近精确结果了。如果想获得更高的精度，只需要增加迭代旋转的次数，所能达到的精度在参考链接中有所讨论，不再赘述。
3.CORDIC算法的Verilog实现
注：为了避免在Verilog中的浮点运算，实现过程中所有的参数都做了放大处理*2^16，角度也放大后映射到了相应的寄存器数值上。
3.0旋转子的实现
该算法在桌面系统基础上的软件实现智能提供验证指导，CORDIC算法主要是在快速硬件电路中使用的，下面按照CORDIC的基本原理，用Verilog进行硬件电路的搭建。
CORDIC算法的核心就是一套迭代公式，也可以说是一个程序操作，程序操作在硬件中可以抽象为一个器件，然后不同的模式重复调用该器件进行最终算法的搭建。旋转模式和向量模式之间主体的操作是一致，区别仅在于初始化输入和迭代判断标志不同。
对于圆坐标模式
。
的值，如果进行旋转模式则根据z的符号判断；如果进行向量模式则根据xy的符号判断。
旋转子的Verilog实现如下：
`timescale 1ns / 1ps
//
// Company: SUSTech
// Engineer: Jack PENG
//
// Create Date: 2022/09/14 17:25:45
// Design Name:
// Module Name: CORDIC_Roter
// Project Name:
// Target Devices:
// Tool Versions:
// Description:
// Basic operation for CORDIC circular;
// MODE=0, rotaiton mode; MODE=1, vector mode;
// Dependencies:
// NONE
// Revision:
// Revision 0.01 - File Created
// Additional Comments:
//
//
module CORDIC_Roter(
SYS_CLK,
RST_N,
Input_x_n_1,
Input_y_n_1,
Input_z_n_1,
Input_sign_n_1,
Input_rote_base,
Output_x_n,
Output_y_n,
Output_z_n,
Output_sign_n
);
parameter ROTE_BASE = 0;
parameter SHIFT_BASE =0;
parameter MODE = 0;
input SYS_CLK;
input RST_N;
input wire signed [31:0] Input_x_n_1;
input wire signed [31:0] Input_y_n_1;
input wire signed [31:0] Input_z_n_1;
input wire [0:0] Input_sign_n_1;
input wire signed [31:0] Input_rote_base;
output reg signed [31:0] Output_x_n;
output reg signed [31:0] Output_y_n;
output reg signed [31:0] Output_z_n;
output reg [0:0] Output_sign_n;
always @ (posedge SYS_CLK or negedge RST_N)
begin
if (!RST_N)
begin
Output_x_n<=1'b0;
Output_y_n<=1'b0;
Output_z_n<=1'b0;
Output_sign_n<=1'b0;
end
else
begin
if (!MODE)
begin // Rotation Mode
if (Input_z_n_1[31])
begin
Output_x_n<=Input_x_n_1+(Input_y_n_1>>>SHIFT_BASE);
Output_y_n<=Input_y_n_1-(Input_x_n_1>>>SHIFT_BASE);
Output_z_n<=Input_z_n_1+Input_rote_base;
Output_sign_n<=Input_sign_n_1;
end
else
begin
Output_x_n<=Input_x_n_1-(Input_y_n_1>>>SHIFT_BASE);
Output_y_n<=Input_y_n_1+(Input_x_n_1>>>SHIFT_BASE);
Output_z_n<=Input_z_n_1-Input_rote_base;
Output_sign_n<=Input_sign_n_1;
end
end
else
begin // Vector Mode
if (!Input_y_n_1[31])
begin
Output_x_n<=Input_x_n_1+(Input_y_n_1>>>SHIFT_BASE);
Output_y_n<=Input_y_n_1-(Input_x_n_1>>>SHIFT_BASE);
Output_z_n<=Input_z_n_1+Input_rote_base;
Output_sign_n<=Input_sign_n_1;
end
else
begin
Output_x_n<=Input_x_n_1-(Input_y_n_1>>>SHIFT_BASE);
Output_y_n<=Input_y_n_1+(Input_x_n_1>>>SHIFT_BASE);
Output_z_n<=Input_z_n_1-Input_rote_base;
Output_sign_n<=Input_sign_n_1;
end
end
end
end
endmodule该模块为待调用的基础模块，没有相应的testbench文件。
3.1旋转模式（Rotation Mode）
3.1.1旋转模式模块文件
上面的CORDIC_Roter模块内置了两种模式，通过MODE进行选择，对于旋转模式例化时MODE=0，模块中默认使用该模式。
`timescale 1ns / 1ps
//
// Company: SUSTech
// Engineer: Jack Peng
//
// Create Date: 2022/09/14 18:33:01
// Design Name: Cordic Algorithm
// Module Name: CORDIC_Rotation
// Project Name:
// Target Devices:
// Tool Versions: 220914
// Description:
// This program is an application of CORDIC algorithm with circlar rotation mode.
// Dependencies: CORDIC_Roter
//
// Revision:
// Revision 0.01 - File Created
// Additional Comments:
//
//
module CORDIC_Rotation(
SYS_CLK,
RST_N,
Phase,
Sin,
Cos,
Error
);
input SYS_CLK;
input RST_N;
input [31:0] Phase;
output [31:0] Sin;
output [31:0] Cos;
output [31:0] Error;
//--- Define the rotation base angle ---
// the rotation is processed in degrees; the values are magnified 2^16 times during processing; totally 16 rotations.
wire signed [31:0] base [15:0];
assign base[00]=32'd2949120; //45deg*2^16
assign base[01]=32'd1740992; //26.5651deg*2^16
assign base[02]=32'd919872; //14.0362deg*2^16
assign base[03]=32'd466944; //7.1250deg*2^16
assign base[04]=32'd234368; //3.5763deg*2^16
assign base[05]=32'd117312; //1.7899deg*2^16
assign base[06]=32'd58688; //0.8952deg*2^16
assign base[07]=32'd29312; //0.4476deg*2^16
assign base[08]=32'd14656; //0.2238deg*2^16
assign base[09]=32'd7360; //0.1119deg*2^16
assign base[10]=32'd3648; //0.0560deg*2^16
assign base[11]=32'd1856; //0.0280deg*2^16
assign base[12]=32'd896; //0.0140deg*2^16
assign base[13]=32'd448; //0.0070deg*2^16
assign base[14]=32'd256; //0.0035deg*2^16
assign base[15]=32'd128; //0.0018deg*2^16
parameter Pipeline = 16; // use a pipline of 16 stage
parameter K=32'h09b74; //K=0.607253*2^16,32'h09b74
parameter PI_half = 5898240;// 90*2^16
parameter PI = 11796480;// 180*2^16
reg signed [31:0] Sin;
reg signed [31:0] Cos;
reg signed [31:0] Error;
reg signed [31:0] x_00,y_00,z_00;
wire signed [31:0] x [16:1];
wire signed [31:0] y [16:1];
wire signed [31:0] z [16:1];
reg signal;
wire [16:1] sign;
// initialize the parameters for cosine and sine compute
/*
x_00=K; is the reciprocal of K_n
y_00=0
z_00=Phase; is the aimed angles
*/
always @ (posedge SYS_CLK or negedge RST_N)
begin
if (!RST_N)
begin
x_00<=1'b0;
y_00<=1'b0;
z_00<=1'b0;
signal<=1'b0;
end
else
begin
x_00<=K;
y_00<=32'd0;
z_00<=Phase;
signal<=1'b0;
if ((Phase>$signed(PI_half))&(Phase<=$signed(PI)))
begin
z_00<=Phase-PI;
signal<=1'b1;
end
else if((Phase<$signed(-PI_half))&(Phase>=$signed(-PI)))
begin
z_00<=Phase+PI;
signal<=1'b1;
end
end
end
//--- generate operation pipeline ---
generate
genvar i;
for (i = 0; i < 16 ; i = i + 1)
begin: roter
if (i==0)
begin
CORDIC_Roter #(.SHIFT_BASE(i),.MODE(0))
rote00 (.SYS_CLK(SYS_CLK),.RST_N(RST_N),.Input_x_n_1(x_00),.Input_y_n_1(y_00),.Input_z_n_1(z_00),.Input_sign_n_1(signal),.Input_rote_base(base[i]),
.Output_x_n(x[i+1]),.Output_y_n(y[i+1]),.Output_z_n(z[i+1]),.Output_sign_n(sign[i+1]));
end
else
begin
CORDIC_Roter #(.SHIFT_BASE(i),.MODE(0))
rote00 (.SYS_CLK(SYS_CLK),.RST_N(RST_N),.Input_x_n_1(x[i]),.Input_y_n_1(y[i]),.Input_z_n_1(z[i]),.Input_sign_n_1(sign[i]),.Input_rote_base(base[i]),
.Output_x_n(x[i+1]),.Output_y_n(y[i+1]),.Output_z_n(z[i+1]),.Output_sign_n(sign[i+1]));
end
end
endgenerate
//--- get the result ---
always @ (posedge SYS_CLK or negedge RST_N)
begin
if (!RST_N)
begin
Cos<=1'b0;
Sin<=1'b0;
Error<=1'b0;
end
else
begin
if (sign[16])
begin
Cos<=~(x[16])+1'b1;
Sin<=~(y[16])+1'b1;
Error<=~(z[16])+1'b1;
end
else
begin
Cos<=x[16];
Sin<=y[16];
Error<=z[16];
end
end
end
endmodule3.1.2旋转模式tb文件
testbench文件如下
`timescale 1ns / 1ps
//
// Company:
// Engineer:
//
// Create Date: 2022/09/14 16:36:00
// Design Name:
// Module Name: tb_CORDIC_Rotation
// Project Name:
// Target Devices:
// Tool Versions:
// Description:
//
// Dependencies:
//
// Revision:
// Revision 0.01 - File Created
// Additional Comments:
//
//
module tb_CORDIC_Rotation;
// Inputs
reg SYS_CLK;
reg RST_N;
reg signed [31:0] cnt;
reg signed [31:0] cnt_n;
reg signed [31:0] Phase;
// reg signed [31:0] Phase_n;
wire [31:0] Sin;
wire [31:0] Cos;
wire [31:0] Error;
real sin_compute;
real cos_compute;
real theta;
// Instantiate the Unit Under Test (UUT)
CORDIC_Rotation uut
(
.SYS_CLK (SYS_CLK ),
.RST_N (RST_N ),
.Phase (Phase ),
.Sin (Sin ),
.Cos (Cos ),
.Error (Error )
);
initial
begin
#0 SYS_CLK = 1'b0;
#1000 RST_N = 1'b0;
#1000 RST_N = 1'b1;
#2000000 $stop;
end
always #1000
begin
SYS_CLK = ~SYS_CLK;
end
always @ (posedge SYS_CLK or negedge RST_N)
begin
if(!RST_N)
cnt <= 1'd0;
else
cnt <= cnt_n;
theta<=(cnt_n-32'd65536*180-32'd65536*34)*3.1415926535897932384626/65536/180;
end
always @ (*)
begin
if(cnt >= 32'd65536*360-1)
cnt_n = 32'd0;
else
cnt_n = cnt + 32'd65536;
end
//生成相位0-359度,Phase[17:16]为相位的象限,Phase[15:10]为相位的值
always @ (posedge SYS_CLK or negedge RST_N)
begin
if(!RST_N)
Phase <= 1'b0;
else
Phase <= cnt-32'd65536*180;
sin_compute<= $sin(theta)*65535;
cos_compute<= $cos(theta)*65535;
end
endmodule
3.1.2仿真波形结果
3.2向量模式（Vector Mode）
3.2.1 向量模式模块文件
对于向量模式例化时MODE=1。
`timescale 1ns / 1ps
//
// Company: SUSTech
// Engineer: Jack PENG
//
// Create Date: 2022/09/14 16:32:07
// Design Name:
// Module Name: CORDIC_Vector
// Project Name:
// Target Devices:
// Tool Versions:
// Description:
// This program is an application of CORDIC algorithm with circlar vector mode.
//
// Dependencies: CORDIC_Roter; set MODE=1;
//
// Revision:
// Revision 0.01 - File Created
// Additional Comments:
//
//
module CORDIC_Vector(
SYS_CLK,
RST_N,
Input_x0,
Input_y0,
Output_xn,
Phase,
Error
);
input SYS_CLK;
input RST_N;
input [31:0] Input_x0;
input [31:0] Input_y0;
output [31:0] Output_xn;
output [31:0] Phase;
output [31:0] Error;
//--- Define the rotation base angle ---
// the rotation is processed in degrees; the values are magnified 2^16 times during processing; totally 16 rotations.
wire signed [31:0] base [15:0];
assign base[00]=32'd2949120; //45deg*2^16
assign base[01]=32'd1740992; //26.5651deg*2^16
assign base[02]=32'd919872; //14.0362deg*2^16
assign base[03]=32'd466944; //7.1250deg*2^16
assign base[04]=32'd234368; //3.5763deg*2^16
assign base[05]=32'd117312; //1.7899deg*2^16
assign base[06]=32'd58688; //0.8952deg*2^16
assign base[07]=32'd29312; //0.4476deg*2^16
assign base[08]=32'd14656; //0.2238deg*2^16
assign base[09]=32'd7360; //0.1119deg*2^16
assign base[10]=32'd3648; //0.0560deg*2^16
assign base[11]=32'd1856; //0.0280deg*2^16
assign base[12]=32'd896; //0.0140deg*2^16
assign base[13]=32'd448; //0.0070deg*2^16
assign base[14]=32'd256; //0.0035deg*2^16
assign base[15]=32'd128; //0.0018deg*2^16
parameter Pipeline = 16; // use a pipline of 16 stage
parameter K=32'h09b74; //K=0.607253*2^16,32'h09b74
parameter PI_half = 5898240;// 90*2^16
parameter PI = 11796480;// 180*2^16
reg signed [31:0] Output_xn;
reg signed [31:0] Phase;
reg signed [31:0] Error;
reg signed [31:0] x_00,y_00,z_00;
wire signed [31:0] x [16:1];
wire signed [31:0] y [16:1];
wire signed [31:0] z [16:1];
reg signal;
wire [16:1] sign;
// initialize the parameters for cosine and sine compute
/*
x_00=K; is the reciprocal of K_n
y_00=0
z_00=Input_x0; is the aimed angles
*/
always @ (posedge SYS_CLK or negedge RST_N)
begin
if (!RST_N)
begin
x_00<=1'b0;
y_00<=1'b0;
z_00<=1'b0;
end
else
begin
x_00<=Input_x0;
y_00<=Input_y0;
z_00<=32'd0;
signal<=1'b0;
end
end
//--- generate operation pipeline ---
generate
genvar i;
for (i = 0; i < 16 ; i = i + 1)
begin: roter
if (i==0)
begin
CORDIC_Roter #(.SHIFT_BASE(i),.MODE(1))
rote00 (.SYS_CLK(SYS_CLK),.RST_N(RST_N),.Input_x_n_1(x_00),.Input_y_n_1(y_00),.Input_z_n_1(z_00),.Input_sign_n_1(signal),.Input_rote_base(base[i]),
.Output_x_n(x[i+1]),.Output_y_n(y[i+1]),.Output_z_n(z[i+1]),.Output_sign_n(sign[i+1]));
end
else
begin
CORDIC_Roter #(.SHIFT_BASE(i),.MODE(1))
rote00 (.SYS_CLK(SYS_CLK),.RST_N(RST_N),.Input_x_n_1(x[i]),.Input_y_n_1(y[i]),.Input_z_n_1(z[i]),.Input_sign_n_1(sign[i]),.Input_rote_base(base[i]),
.Output_x_n(x[i+1]),.Output_y_n(y[i+1]),.Output_z_n(z[i+1]),.Output_sign_n(sign[i+1]));
end
end
endgenerate
always @ (posedge SYS_CLK or negedge RST_N)
begin
if (!RST_N)
begin
Phase<=1'b0;
Output_xn<=1'b0;
Error<=1'b0;
end
else
begin
if (sign[16])
begin
Output_xn<=~(x[16])+1'b1;
Phase<=~(z[16])+1'b1;
Error<=~(y[16])+1'b1;
end
else
begin
Output_xn<=x[16];
Phase<=z[16];
Error<=y[16];
end
end
end
endmodule3.2.2向量模式tb文件
testbench文件如下：
`timescale 1ns / 1ps
//
// Company:
// Engineer:
//
// Create Date: 2022/09/14 16:33:28
// Design Name:
// Module Name: tb_CORDIC_Vector
// Project Name:
// Target Devices:
// Tool Versions:
// Description:
//
// Dependencies:
//
// Revision:
// Revision 0.01 - File Created
// Additional Comments:
//
//
module tb_CORDIC_Vector;
// Inputs
reg SYS_CLK;
reg RST_N;
reg signed [31:0] cnt;
reg signed [31:0] cnt_n;
reg signed [31:0] Input_x0;
reg signed [31:0] Input_y0;
reg signed [31:0] Input_yy;
wire [31:0] Output_xn;
wire [31:0] Phase;
wire [31:0] Error;
real atan;
// Instantiate the Unit Under Test (UUT)
CORDIC_Vector uut
(
.SYS_CLK (SYS_CLK ),
.RST_N (RST_N ),
.Input_x0 (Input_x0 ),
.Input_y0 (Input_y0 ),
.Output_xn (Output_xn ),
.Phase (Phase ),
.Error (Error )
);
parameter K=0.607253; //K=0.607253*2^16,32'h09b74
initial
begin
#0 SYS_CLK = 1'b0;
#10000 RST_N = 1'b0;
#10000 RST_N = 1'b1;
#100000000 $stop;
end
always #10000
begin
SYS_CLK = ~SYS_CLK;
end
always @ (posedge SYS_CLK or negedge RST_N)
begin
if(!RST_N)
cnt <= 1'd0;
else
cnt <= cnt_n;
end
always @ (*)
begin
if(cnt >= 32'd65536*10)
cnt_n = 32'd0;
else
cnt_n = cnt + 32'd655;
end
//生成相位0-359度,Phase[17:16]为相位的象限,Phase[15:10]为相位的值
always @ (posedge SYS_CLK or negedge RST_N)
begin
if(!RST_N)
begin
Input_x0 <= 1'b1;
Input_y0 <= 1'b1;
end
else
begin
Input_x0 <= 32'd65536;
Input_y0 <= (cnt-32'd65536*5);
Input_yy<=(cnt-32'd655*17-32'd65536*5);
atan<=$atan2(Input_yy,Input_x0)*65536*180/3.1415926535897932384626;
//tan<=Input_y0>>>16;
end
end
endmodule
3.2.3仿真波形结果
