Explanation and FPGA implementation of CORDIC algorithm (circular coordinate system)

CORDIC Algorithm Explanation and Verilog Simulation (Circular Coordinate System)

Explanation of the principle of CORDIC algorithm

The CORDIC (Coordinate Rotation Digital Computer) algorithm, or Coordinate Rotation Digital Computing method, was first proposed by .Volder1 in 1959, and is mainly used for trigonometric, hyperbolic, exponential, and logarithmic calculations.

pseudorotation

In the Cartesian coordinate plane (left below) is given by$({x_1},{y_1})$ Rotate θ angle to$({x_2},{y_2})$ Get:$({\hat x_2},{\hat y_2})$ ；

propose a factor$\cos \theta$ , the equations are transformed:$\left\{ {\matrix{ {{x_2} = \cos \theta ({x_1} - {y_1}\tan \theta )} \cr {{y_2} = \cos \theta ({y_1} + {x_1}\tan \theta )} \cr } } \right.$；

To be removed$\cos \theta$ term to obtain the "pseudo-rotation" equation$\left\{ {\matrix{ {{{\hat x}_2} = {x_1} - {y_1}\tan \theta } \cr {{{\hat y}_2} = {y_1} + {x_1}\tan \theta } \cr } } \right.$。

After "pseudo-rotation", the modulus of the vector R is increased by a factor of $1/\cos \theta $ (the angle remains the same).

Angle Totalizer

For FPGA hardware implementation (the tangent term needs to be changed to a shift operation): set the rotation angle θ with $ \tan {\theta ^i} = {2^{ - i}}$;

Therefore, the equation can be converted to$\left\{ {\matrix{ {{{\hat x}_{_2}} = {x_1} - {y_1}{2^{ - i}}} \cr {{{\hat y}_{_2}} = {y_1} + {x_1}{2^{ - i}}} \cr } } \right.$ maybe$\left[ {\matrix{ {{{\hat x}_{_2}}} \cr {{{\hat y}_{_2}}} \cr } } \right] = \left[ {\matrix{ 1 & { - {2^{ - i}}} \cr {{2^{ - i}}} & 1 \cr } } \right]\left[ {\matrix{ {{x_1}} \cr {{y_1}} \cr } } \right]$。

where the matrix$\left[ {\matrix{ 1 & { - {2^{ - i}}} \cr {{2^{ - i}}} & 1 \cr } } \right]$ The splitting into multiple similar matrix products, i.e., the rotation angle θ , can be performed into multiple small rotations (the corresponding table of inverse tangent angles is shown below).

Since the rotation angle θ can be an arbitrary value, the rotation transformation is realized using an iterative algorithm, i.e., multiple angle iterations relate to an infinite convergence to the target θ angle (limited by the θ rotation angle). Take the 55° degree rotation angle as an example to approximate 55° = 45.0° + 26.6° -14.0°- 7.1° + 3.6° + 1.8° - 0.9°.

The rotation process requires the introduction of a judgment factor${d_i}$ , which is used to determine the direction of angular rotation.

Based on the judgment factor${d_i}$ to set an angle accumulator: $\eqalign{
& {z^{(i + 1)}} = {z^{(i)}} - {d_i}{\theta ^{(i)}} \cr
& where:{d_i} = \pm 1 \cr} $, where z (the difference in rotation angles) tends infinitely to zero.

and the pseudo-rotation can be expressed as$\left\{ {\matrix{ {{x^{(i + 1)}} = {x^{(i)}} - {d_i}({2^{ - i}}{y^{(i)}})} \cr {{y^{(i + 1)}} = {y^{(i)}} + {d_i}({2^{ - i}}{x^{(i)}})} \cr } } \right.$。

quadrant preprocessing

Of course, the direction of each rotation affects the final cumulative angle to be rotated, which is roughly in the range: $ - 99.7 \le \theta \le 99.7$. For angles outside of this range, it is necessary to "preprocess" using the trigonometric constant transformation, i.e., quadrant determination.

Thus, the original algorithm regularizes to represent the iterative shift-add algorithm using pseudo-rotations of vectors, i.e:$\left\{ {\matrix{ {{x^{(i + 1)}} = {x^{(i)}} - {d_i}({2^{ - i}}{y^{(i)}})} \cr {{y^{(i + 1)}} = {y^{(i)}} + {d_i}({2^{ - i}}{x^{(i)}})} \cr {{z^{(i + 1)}} = {z^{(i)}} - {d_i}{\theta ^{(i)}}} \cr } } \right.$，

As mentioned earlier, the "pseudo-rotation" operation ignores each iteration of the$\cos \theta$Items that end up getting${x^{(n)}},{y^{(n)}}$ It's been stretched.${k_n}$times (multiplier)

${k_n} = \prod\limits_n {({1 \over {\cos {\theta ^{(i)}}}})} = \prod\limits_n {(\sqrt {1 + {2^{( - 2i)}}} )} $ (Stretch factor).

treat (sb a certain way)${k_n}$ Find the infinite product, ${k_n} = \prod\limits_n {(\sqrt {1 + {2^{( - 2i)}}} )} \to 1.6476,as:n \to \infty $($1/{k_n} = 0.6073$ ）

If the number of iterations performed is known, it is straightforward to obtain the${k_n}$ Final value.

Regarding the circular coordinate system, CORDIC algorithm applications include both rotational and vector modes:

rotary mode

Application scenario: the phase angle angle is known and its sine and cosine values are computed using the Cordic algorithm.

Specific process: judgment factor${d_i}{\rm{ = sign}}({z^{(i)}}) \Rightarrow {z^{(i)}} \to 0$which is obtained after N iterations$\left\{ {\matrix{ {{x^{(n)}} = {k_n}({x^{(0)}}\cos {z^{(0)}} - {y^{(0)}}\sin {z^{(0)}})} \cr {{y^{(n)}} = {k_n}({y^{(0)}}\cos {z^{(0)}} + {x^{(0)}}\sin {z^{(0)}})} \cr {{z^{(n)}} = 0} \cr } } \right.$（ ${z^{(0)}}$ = θ) by setting${x^{(0)}} = {1 \over {{k_n}}}{{\rm{y}}^{(0)}} = 0$The final result is $ \cos \theta, \sin \theta $ .

vector pattern

Application scenario: coordinates are known and phase angle and magnitude are calculated using the cordic algorithm.

Specific process: the polar coordinate system transformed by the right-angle coordinate system, the iterative process changes to$\left\{ {\matrix{ {{x^{(i + 1)}} = {x^{(i)}} - {d_i}({2^{ - i}}{y^{(i)}})} \cr {{y^{(i + 1)}} = {y^{(i)}} + {d_i}({2^{ - i}}{x^{(i)}})} \cr {{z^{(i + 1)}} = {z^{(i)}} - {d_i}{\theta ^{(i)}}} \cr } } \right.$，

where the judgment factor${d_i}{\rm{ = - sign}}({x^{(i)}}{y^{(i)}}) \Rightarrow {y^{(i)}} \to 0$, N iterations are obtained:$\left\{ {\matrix{ {{x^{(n)}} = {k^{(n)}}\sqrt {x_0^2 + y_0^2} } \cr {{y^{(n)}} = 0} \cr {{z^{(n)}} = {z^{(0)}} + {{\tan }^{ - 1}}({y_0}/{x_0})} \cr {{k^{(n)}} = \prod\limits_n {\sqrt {1 + {2^{ - 2i}}} } } \cr } } \right.$，

By setting the${x^{(0)}} = 1，{z^{(0)}} = 0$which can be finally solved for${\tan ^{ - 1}}{y^{(0)}}$。

Verilog HDL Implementation of CORDIC

be directed against$\left\{ {\matrix{ {{x^{(i + 1)}} = {x^{(i)}} - {d_i}({2^{ - i}}{y^{(i)}})} \cr {{y^{(i + 1)}} = {y^{(i)}} + {d_i}({2^{ - i}}{x^{(i)}})} \cr {{z^{(i + 1)}} = {z^{(i)}} - {d_i}{\theta ^{(i)}}} \cr } } \right.$ The calculation requires 2 shifts per iteration.$({x^{(i)}{,y^{(i)}}})$ 1 look-up table${\theta ^{(i)}}$, 3 additions (x, y, z cumulative).

The corresponding CORDIC hardware structure is as follows:

Matlab code implementation in Cordic-rotation mode:

Click to view code

%% ***********************************************************************************
% In the circular coordinate system: Cordic-rotation mode
% The phase angle angle is known and its sine and cosine values are calculated. The basic formula is as follows:
% x(k+1) = x(k) - d(k)*y(k)*2^(-k)
% y(k+1) = y(k) + d(k)*x(k)*2^(-k)
% z(k) = z(k) - d(k)*actan(2^(-k))
%% ***********************************************************************************
clear;close all;clc.

angle = 30; %Set the rotation angle

% Initialize -------------------------------
N = 16; %Number of iterations
tan_table = 2.^-(0 : N-1); % set the rotation angle.
angle_LUT = atan(tan_table); % create arctan&angle lookup table

An = 1; for k = 0 : N-1
for k = 0 : N-1
    An = An*(1/sqrt(1 + 2^(-2*k))); %Build arctan&angle lookup table.
end
Kn = 1/An;%Calculate the normalized scaling factor parameter: Kn = 1.6476,1/Kn = 0.6073

Xn = 1/Kn; %Start rotation with respect to the X-axis
Yn = 0;.

Zi = angle/180*pi; %Angle converted to radians

% cordic algorithm calculated by -------------------------------
if (Zi > pi/2) % first do quadrant judgment to get phase compensation value
    Zi = Zi - pi.
    sign_x = -1; sign_y = -1; sign_x = -1
    sign_y = -1; elseif (Zi < pi/2)
elseif (Zi < -pi/2)
    Zi = Zi + pi; sign_x = -1; sign_y = -1; elseif (Zi < -pi/2)
    sign_x = -1; sign_y = -1; elseif (Zi < -pi/2)
    sign_y = -1; elseif (Zi < -pi/2)
sign_y = -1; elseif (Zi < -pi/2)
    sign_x = 1; sign_y = 1; else
    sign_x = 1; sign_y = 1; else
sign_x = 1; sign_y = 1; end

 for k = 0 : N-1 % start of iteration
        Di = sign(Zi).

        x_temp = Xn; Xn = x_temp - Di*Yn*2^(-k); x_temp
        Xn = x_temp - Di*Yn*2^(-k);
        Yn = Yn + Di*x_temp*2^(-k);
        Zi = Zi - Di*angle_LUT(k+1); end
Zi = Zi - Di*angle_LUT(k+1); end

cos_out = sign_x*Xn; %cosine output
sin_out = sign_y*Yn; %sine output

Verilog HDL inrotary modeUnder the program:

Click to view code

module Cordic_rotate_mode(
    input sys_clk ,
    input sys_rst ,

    input signed [31:0] angle ,

    output reg [31:0] cosout ,
    output reg [31:0] sinout
);

//Angle of rotation lookup table
wire [31:0]rot[15:0];

assign rot[0] = 32'd2949120 ; //45.0000degree (angles, temperature etc)*2^16
assign rot[1] = 32'd1740992 ; //26.5651degree (angles, temperature etc)*2^16
assign rot[2] = 32'd919872 ; //14.0362degree (angles, temperature etc)*2^16
assign rot[3] = 32'd466944 ; //7.1250degree (angles, temperature etc)*2^16
assign rot[4] = 32'd234368 ; //3.5763degree (angles, temperature etc)*2^16
assign rot[5] = 32'd117312 ; //1.7899degree (angles, temperature etc)*2^16
assign rot[6] = 32'd58688 ; //0.8952degree (angles, temperature etc)*2^16
assign rot[7] = 32'd29312 ; //0.4476degree (angles, temperature etc)*2^16
assign rot[8] = 32'd14656 ; //0.2238degree (angles, temperature etc)*2^16
assign rot[9] = 32'd7360 ; //0.1119degree (angles, temperature etc)*2^16
assign rot[10] = 32'd3648 ; //0.0560degree (angles, temperature etc)*2^16
assign rot[11] = 32'd1856 ; //0.0280degree (angles, temperature etc)*2^16
assign rot[12] = 32'd896 ; //0.0140degree (angles, temperature etc)*2^16
assign rot[13] = 32'd448 ; //0.0070degree (angles, temperature etc)*2^16
assign rot[14] = 32'd256 ; //0.0035degree (angles, temperature etc)*2^16
assign rot[15] = 32'd128 ; //0.0018degree (angles, temperature etc)*2^16

//FSM_parameter
localparam IDLE = 2'd0;
localparam WORK = 2'd1;
localparam ENDO = 2'd2;

reg [1:0] state ;
reg [1:0] next_state ;
reg [3:0] cnt;


always @(posedge sys_clk or negedge sys_rst)begin
    if(!sys_rst)
        next_state <= IDLE;
    else begin
        state <= next_state;
        case(state)
            IDLE:next_state <= WORK;
            WORK:next_state <= cnt == 15 ? ENDO:WORK;
            ENDO:next_state <= IDLE;
            default:next_state <= IDLE;
        endcase
    end
end


reg signed [31:0] x_shift;
reg signed [31:0] y_shift;
reg signed [31:0] z_rot;

wire D_sign;
assign D_sign= z_rot[31];

always @(posedge sys_clk) begin
    case(state)
    IDLE:
        begin
            x_shift <= 32'd39800;
            y_shift <= 32'd0;
            z_rot <= (angle<<16);
        end
        
    WORK:
        if(D_sign)begin
            x_shift <= x_shift + (y_shift>>>cnt);
            y_shift <= y_shift - (x_shift>>>cnt);
            z_rot <= z_rot + rot[cnt];
        end
        else begin
            x_shift <= x_shift - (y_shift>>>cnt);
            y_shift <= y_shift + (x_shift>>>cnt);
            z_rot <= z_rot - rot[cnt];
        end
        
    ENDO:
        begin
            cosout <= x_shift;
            sinout <= y_shift;
        end
        
    default :;
    endcase
end

always @(posedge sys_clk or negedge sys_rst) begin
    if(!sys_rst)
        cnt <= 4'd0;
    else if(state == IDLE && next_state == WORK)
        cnt <= 4'd0;
    else if(state==WORK)begin
        if(cnt<4'd15)
            cnt <= cnt + 1'b1;
        else
            cnt <= cnt;
    end
    else
        cnt <= 4'd0;
end

endmodule

Set a variety of angle values and the simulation is shown below:

In Cordic-vector patternunder the Matlab code implementation:

Click to view code

%% ***********************************************************************************
% In circular coordinate system: Cordic-vector mode
% The coordinates are known and the phase angle and magnitude are calculated using the cordic algorithm. The basic formula is as follows:
% x(k+1) = x(k) - d(k)*y(k)*2^(-k)
% y(k+1) = y(k) + d(k)*x(k)*2^(-k)
% z(k) = z(k) - d(k)*actan(2^(-k))
%% ***********************************************************************************
clear;close all;clc.
% Initialize ----------------------------------------
Xn = -1.
Yn = sqrt(3);

Zi = 0; Di = 0;; % Initialize
Di = 0.

N = 16; % number of iterations
tan_table = 2.^-(0 : N-1);
angle_LUT = atan(tan_table);

An = 1; for k = 0 : N-1
for k = 0 : N-1
    An = An*(1/sqrt(1 + 2^(-2*k))); end
An = An*(1/sqrt(1 + 2^(-2*k)); end
Kn = 1/An;%Calculate the normalized scaling factor parameter: Kn = 1.6476,1/Kn = 0.6073

% cordic algorithm calculation -------------------------------
if (Xn==0 && Yn==0) %move to origin, not rotated by angle
    radian_out = 0;
    amplitude_out = 0.
else % first do quadrant judgment, get phase compensation value
    if (Xn > 0) %first quadrant: (-pi/2,0)/(0,pi/2)-->Zn
        phase_shift = 0;
    elseif (Yn < 0) %Quadrant 3: (-pi,-pi/2)-->pre-rotate -pi,Zn+pi/2
        phase_shift = -pi.
    else %Quadrant 2: (pi/2,pi)-->pre-rotate pi,Zn-pi/2
        phase_shift = pi; end
    phase_shift = pi; end

    for k = 0 : N-1 % start iteration
        Di = -sign(Xn*Yn);

        x_temp = Xn.
        Xn = x_temp - Di*Yn*2^(-k);
        Yn = Yn + Di*x_temp*2^(-k);
        Zi = Zi - Di*angle_LUT(k+1); end
    Zi = Zi - Di*angle_LUT(k+1); end
    radian_out = Zi + phase_shift; % radian out
    amplitude_out = abs(Xn)/Kn; %amplitude_out
amplitude_out = abs(Xn)/Kn; %amplitude_out

angle_out = radian_out*180/pi; %phase_out: angle in degrees = angle in radians x pi/180

Verilog HDL in vector mode, program:

Click to view code

module Cordic_vector_mode(
    input sys_clk ,
    input sys_rst ,

    input signed [31:0] x ,
    input signed [31:0] y ,

    output reg [31:0] phase ,
    output reg [31:0] mo_value
);


//Angle of rotation lookup table
wire [31:0]rot[15:0];

assign rot[0] = 32'd2949120 ; //45.0000degree (angles, temperature etc)*2^16
assign rot[1] = 32'd1740992 ; //26.5651degree (angles, temperature etc)*2^16
assign rot[2] = 32'd919872 ; //14.0362degree (angles, temperature etc)*2^16
assign rot[3] = 32'd466944 ; //7.1250degree (angles, temperature etc)*2^16
assign rot[4] = 32'd234368 ; //3.5763degree (angles, temperature etc)*2^16
assign rot[5] = 32'd117312 ; //1.7899degree (angles, temperature etc)*2^16
assign rot[6] = 32'd58688 ; //0.8952degree (angles, temperature etc)*2^16
assign rot[7] = 32'd29312 ; //0.4476degree (angles, temperature etc)*2^16
assign rot[8] = 32'd14656 ; //0.2238degree (angles, temperature etc)*2^16
assign rot[9] = 32'd7360 ; //0.1119degree (angles, temperature etc)*2^16
assign rot[10] = 32'd3648 ; //0.0560degree (angles, temperature etc)*2^16
assign rot[11] = 32'd1856 ; //0.0280degree (angles, temperature etc)*2^16
assign rot[12] = 32'd896 ; //0.0140degree (angles, temperature etc)*2^16
assign rot[13] = 32'd448 ; //0.0070degree (angles, temperature etc)*2^16
assign rot[14] = 32'd256 ; //0.0035degree (angles, temperature etc)*2^16
assign rot[15] = 32'd128 ; //0.0018degree (angles, temperature etc)*2^16

//FSM_parameter
localparam IDLE = 2'd0;
localparam WORK = 2'd1;
localparam ENDO = 2'd2;

reg [1:0] state ;
reg [1:0] next_state ;
reg [3:0] cnt;

reg signed [31:0] x_shift;
reg signed [31:0] y_shift;
reg signed [31:0] z_rot;


always @(posedge sys_clk or negedge sys_rst)begin
    if(!sys_rst)
        next_state <= IDLE;
    else begin
        state <= next_state;
        case(state)
            IDLE:next_state <= WORK;
            WORK:next_state <= cnt == 15 ? ENDO:WORK;
            ENDO:next_state <= IDLE;
            default:next_state <= IDLE;
        endcase
    end
end

wire D_sign;
assign D_sign=~y_shift[31];


always @(posedge sys_clk) begin
    case(state)
    IDLE:
        begin
            x_shift <= x;
            y_shift <= y;
            z_rot <= 0;
        end
        
    WORK:
        if(D_sign)begin
            x_shift <= x_shift + (y_shift>>>cnt);
            y_shift <= y_shift - (x_shift>>>cnt);
            z_rot <= z_rot + rot[cnt];
        end
        else begin
            x_shift <= x_shift - (y_shift>>>cnt);
            y_shift <= y_shift + (x_shift>>>cnt);
            z_rot <= z_rot - rot[cnt];
        end
        
    ENDO:
        begin
            phase <= z_rot>>>16;
            mo_value <= (x_shift>>>16)*0.6073;
        end
        
    default :;
    endcase
en


always @(posedge sys_clk or negedge sys_rst) begin
    if(!sys_rst)
        cnt <= 4'd0;
    else if(state == IDLE && next_state == WORK)
        cnt <= 4'd0;
    else if(state==WORK)begin
        if(cnt<4'd15)
            cnt <= cnt + 1'b1;
        else
            cnt <= cnt;
    end
    else
        cnt <= 4'd0;
end

endmodule

Setting three different x, y values, the simulation is shown below:

The Verilog program modules used in this post can be found in the Gitee repository link in the left column of this page:/silly-big-head/little-mouse-funnyhouse/tree/FPGA-Verilog/