当前位置：首页 > 后端开发 > 正文内容

CORDIC算法解说及verilog HDL完成（圆坐标系）

邻居的猫1个月前 (12-09)后端开发1622

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/

扫描二维码推送至手机访问。

版权声明：本文由51Blog发布，如需转载请注明出处。

本文链接：https://www.51blog.vip/?id=254

标签: FPGA-Verilog类
分享给朋友：
返回列表

上一篇：FM的正交解调法

下一篇：FPGA驱动adc128s052的几个问题

“CORDIC算法解说及verilog HDL完成（圆坐标系）” 的相关文章

Go Vue3 CMS办理后台(前后端别离形式)

Go Vue3 CMS办理后台(前后端别离形式)

本后台运用前后端别离形式开发，前端UI为Vue3+Ant Design Vue，后端Api为Go+Gin，解耦前后端逻辑，使开发更专心 技能栈 前端：Vue3，Ant Design Vue，Axios，分页，OTP动态码登录 后端：Go，Gin，Gorm，Mysql，Docker，JWT，跨域，...

【日记】咱们行发工资真的便是 Black Box……（577 字）

【日记】咱们行发工资真的便是 Black Box……（577 字）

正文 今日头好油…… 昨日应付完了真实太晚，就没洗澡。现在的头几乎无法看…… 回想了一下，今日如同什么都没干。字面意义上的。今日新行长下来，带了一堆东西。去帮了忙。他看见我还一愣。估量是头太油了……. 发工资了。市分行的搭档问我怎样比跟我同一批进来的人高那么多。你问我我也不知道啊…… 人力也不发个工...

【日记】自己心里戏很多（笑（968 字）

【日记】自己心里戏很多（笑（968 字）

正文 　　本来想手写来着，成果找了快一个小时的图。没找到。抛弃了。时间也不大够用了，就不手写了。 　　找图首要是由于一件事——今日遽然告诉要拍证件照。 　　我特别疑惑，之前不是拍过了吗，并且也没怎样用到，这东西。如同必需要从头拍，不知道为什么。并且正午才告诉。 　　还必需要打领带…… 　　谁没事儿系...

shell （3）脚本参数传递与数学运算

shell （3）脚本参数传递与数学运算

🌟声明🌟 红客全栈教程 学习视频来自UP 泷羽sec，如涉及侵权马上删除文章 以下只涉及学习内容，其他都与本人无关，切莫逾越法律红线，否则后果自负。 星河飞雪网络安全人才培养计划，绝对零区，公益免费教学！没有网络安全，就没有国家安全！ 脚本参数如何传递？ echo 执行的文件名是:$0 echo...

go人体,基因本体与生物信息学的桥梁

go人体,基因本体与生物信息学的桥梁

您好，关于“go人体”的搜索结果中，大部分内容与围棋相关，并未找到直接与“人体”相关的信息。请问您是否需要了解有关围棋的内容，例如围棋的基本规则、历史背景、艺术价值等？如果您有其他具体需求，请告诉我，我会尽力为您提供帮助。探索GO人体：基因本体与生物信息学的桥梁随着生物信息学的发展，基因本体（Gen...

java面试宝典,java官网

java面试宝典,java官网

1. JavaGuide 这是一个全面的Java学习与面试指南，涵盖了Java基础、集合、IO、并发、JVM、新特性等多方面的知识。非常适合准备Java面试的朋友使用。 2. 2024最全Java面试八股文 这篇文章分享了一套详细的Java面试手册，涵盖了MyBatis、Zooke...

Copyright © 2024-2028 51blog All Rights Reserved

XML地图