CORDIC演算法解釋及verilog HDL實現(圓座標系)

Handat發表於2024-08-18

CORDIC演算法原理闡述

CORDIC(Coordinate Rotation Digital Computer)演算法,即座標旋轉數字計算方法,是J.D.Volder1於1959年首次提出,主要用於三角函式、雙曲線、指數、對數的計算。

偽旋轉

在笛卡爾座標平面(下方左圖)由 \(({x_1},{y_1})\) 旋轉 θ 角度至 \(({x_2},{y_2})\) 得到:\(({\hat x_2},{\hat y_2})\)

提出因數 \(\cos \theta\) ,方程轉化為:\(\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.\)

待去除 \(\cos \theta\) 項,得到“偽旋轉”公式\(\left\{ {\matrix{ {{{\hat x}_2} = {x_1} - {y_1}\tan \theta } \cr {{{\hat y}_2} = {y_1} + {x_1}\tan \theta } \cr } } \right.\)

經“偽旋轉”後,向量 R 模值將增加 $1/\cos \theta $ 倍(角度保持一致)。
image

角度累加器

為便於FPGA硬體實現(正切項需改為移位操作):以 $ \tan {\theta ^i} = {2^{ - i}}$ 設定旋轉角度 θ ;

故方程可轉換為\(\left\{ {\matrix{ {{{\hat x}_{_2}} = {x_1} - {y_1}{2^{ - i}}} \cr {{{\hat y}_{_2}} = {y_1} + {x_1}{2^{ - i}}} \cr } } \right.\)\(\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]\)

其中矩陣 \(\left[ {\matrix{ 1 & { - {2^{ - i}}} \cr {{2^{ - i}}} & 1 \cr } } \right]\) 可進行拆分為多個類似矩陣乘積,即旋轉角度 θ ,可拆分為多次小的旋轉(下圖為對應的反正切角度表)。
image

由於旋轉角度 θ 可為任意值,故將旋轉變換採用迭代演算法實現,即多次角度迭代關係無限趨近於目標θ角度(以 θ 旋轉角度限制)。以55°度旋轉角為例逼近55° = 45.0° + 26.6° -14.0°- 7.1° + 3.6° + 1.8° - 0.9°。

旋轉過程需引入一個判決因子 \({d_i}\) ,用於確定角度旋轉的方向。

根據判決因子 \({d_i}\) 來設定一個角度累加器:$\eqalign{
& {z^{(i + 1)}} = {z^{(i)}} - {d_i}{\theta ^{(i)}} \cr
& where:{d_i} = \pm 1 \cr} $,其中z(旋轉角度差)無限趨近於0。

並且偽旋轉可表示為\(\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.\)

象限預處理

當然,每次旋轉的方向都影響到最終要旋轉的累積角度,角度範圍大致為: $ - 99.7 \le \theta \le 99.7$。對於範圍外的角度,需要使用三角恆等式轉化進行“預處理”,即象限判斷。

image

因此,原始演算法規整為使用向量的偽旋轉來表示迭代移位-相加演算法,即:\(\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.\)

前面提到了,在進行“偽旋轉”操作時,每次迭代運算都忽略了\(\cos \theta\)項,最終得到的 \({x^{(n)}},{y^{(n)}}\) 被伸縮了 \({k_n}\)

${k_n} = \prod\limits_n {({1 \over {\cos {\theta ^{(i)}}}})} = \prod\limits_n {(\sqrt {1 + {2^{( - 2i)}}} )} $ (伸縮因子)。

\({k_n}\) 求無限積,${k_n} = \prod\limits_n {(\sqrt {1 + {2^{( - 2i)}}} )} \to 1.6476,as:n \to \infty $( \(1/{k_n} = 0.6073\)

若已知執行的迭代次數,便可直接求得 \({k_n}\) 最終值。


關於圓座標系下,CORDIC演算法應用包括旋轉模式和向量模式兩種:

旋轉模式

應用場景:已知相角angle,用Cordic演算法計算其正弦和餘弦值。

具體過程:判決因子\({d_i}{\rm{ = sign}}({z^{(i)}}) \Rightarrow {z^{(i)}} \to 0\),N次迭代後得到\(\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)}}\) = θ)透過設定 \({x^{(0)}} = {1 \over {{k_n}}}{{\rm{y}}^{(0)}} = 0\),可最終求到 $\cos \theta、 \sin \theta $ 。

向量模式

應用場景:已知座標,用cordic演算法計算相角和幅值。

具體過程:直角座標系轉換的極座標系,迭代過程變化為\(\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.\)

其中判決因子 \({d_i}{\rm{ = - sign}}({x^{(i)}}{y^{(i)}}) \Rightarrow {y^{(i)}} \to 0\),N次迭代得到:\(\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.\)

透過設定\({x^{(0)}} = 1,{z^{(0)}} = 0\),可最終求得 \({\tan ^{ - 1}}{y^{(0)}}\)

Verilog HDL實現CORDIC

針對\(\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.\) ,每次迭代計算需要2次移位 \(({x^{(i)}{,y^{(i)}}})\) 、1次查詢表\({\theta ^{(i)}}\)、3次加法(x、y、z累加)。

對應的CORDIC硬體結構如下:

image

在Cordic—旋轉模式下,Matlab程式碼實現:

點選檢視程式碼
%% ***********************************************************************************
%     圓座標系下:Cordic—旋轉模式
%     已知相角angle,計算其正弦和餘弦值。基本公式如下:
%     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;    %設定旋轉角度

% 初始化-------------------------------
N = 16;  %迭代次數
tan_table = 2.^-(0 : N-1);
angle_LUT = atan(tan_table);    %建立arctan&angle查詢表

An = 1;
for k = 0 : N-1
    An = An*(1/sqrt(1 + 2^(-2*k)));  
end
Kn = 1/An;%計算歸一化伸縮因子引數:Kn = 1.6476,1/Kn = 0.6073

Xn = 1/Kn; %相對於X軸上開始旋轉
Yn = 0;

Zi = angle/180*pi;  %角度轉化為弧度

% cordic演算法計算-------------------------------
if (Zi > pi/2)  % 先做象限判斷,得到相位補償值
    Zi = Zi - pi;
    sign_x = -1;
    sign_y = -1;
elseif (Zi < -pi/2)
    Zi = Zi + pi;
    sign_x = -1;
    sign_y = -1;
else
    sign_x = 1;
    sign_y = 1;
end

 for k = 0 : N-1   % 迭代開始
        Di = sign(Zi);
     
        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

cos_out = sign_x*Xn;  %餘弦輸出
sin_out = sign_y*Yn;  %正弦輸出


Verilog HDL在旋轉模式下,程式:

點選檢視程式碼
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
);

//旋轉角度查詢表
wire   [31:0]rot[15:0];

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


設定多種角度值,模擬如下圖:

image

在Cordic—向量模式下,Matlab程式碼實現:

點選檢視程式碼
%% ***********************************************************************************
%     圓座標系下:Cordic—向量模式
%     已知座標,用cordic演算法計算相角和幅值。基本公式如下:
%     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;
% 初始化----------------------------------------
Xn = -1;
Yn = sqrt(3);

Zi = 0;
Di = 0;

N = 16;  %迭代次數
tan_table = 2.^-(0 : N-1);
angle_LUT = atan(tan_table);

An = 1;
for k = 0 : N-1
    An = An*(1/sqrt(1 + 2^(-2*k)));  
end
Kn = 1/An;%計算歸一化伸縮因子引數:Kn = 1.6476,1/Kn = 0.6073

% cordic演算法計算-------------------------------
if (Xn==0 && Yn==0)     %移至原點,未旋轉角度
    radian_out = 0;
    amplitude_out = 0;
else  % 先做象限判斷,得到相位補償值
    if (Xn > 0)         %第一、四象限:(-pi/2,0)/(0,pi/2)-->Zn
        phase_shift = 0;
    elseif (Yn < 0)     %第三象限:(-pi,-pi/2)-->預旋轉-pi,Zn+pi/2
        phase_shift = -pi;
    else                %第二象限:(pi/2,pi)-->預旋轉pi,Zn-pi/2
        phase_shift = pi;
    end
  
    for k = 0 : N-1   % 迭代開始
        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
    radian_out = Zi + phase_shift; %弧度輸出
    amplitude_out = abs(Xn)/Kn;  %幅值輸出
end

angle_out = radian_out*180/pi;  %相角輸出:角度(度)=角度(弧度)x pi/180


Verilog HDL在向量模式下,程式:

點選檢視程式碼
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
);


//旋轉角度查詢表
wire   [31:0]rot[15:0];

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


設定三種不同x,y值,模擬如下圖:

image


本篇文章中使用的Verilog程式模組,若有需見網頁左欄Gitee倉庫連結:https://gitee.com/silly-big-head/little-mouse-funnyhouse/tree/FPGA-Verilog/

相關文章