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 $ 倍（角度保持同等）。
角度累加器
为便于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]\) 可进行拆分为多个类似矩阵乘积，即旋转角度 θ ，可拆分为多次小的旋转（下图为对应的反正切角度表）。
由于旋转角度 θ 可为恣意值，故将旋转变更采用迭代算法实现，即多次角度迭代关系无限趋近于目标θ角度（以 θ 旋转角度限制）。以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$。对于范围外的角度，需要利用三角恒等式转化进行“预处理惩罚”，即象限判断。
因此，原始算法规整为利用向量的伪旋转来表示迭代移位-相加算法，即：\(\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硬件结构如下：
点击查看代码- %% ***********************************************************************************
- % 圆坐标系下：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; %正弦输出
点击查看代码[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);//旋转角度查找表wire [31:0]rot[15:0];assign rot[0] = 32'd2949120 ; //45.0000度*2^16assign rot[1] = 32'd1740992 ; //26.5651度*2^16assign rot[2] = 32'd919872 ; //14.0362度*2^16assign rot[3] = 32'd466944 ; //7.1250度*2^16assign rot[4] = 32'd234368 ; //3.5763度*2^16assign rot[5] = 32'd117312 ; //1.7899度*2^16assign rot[6] = 32'd58688 ; //0.8952度*2^16assign rot[7] = 32'd29312 ; //0.4476度*2^16assign rot[8] = 32'd14656 ; //0.2238度*2^16assign rot[9] = 32'd7360 ; //0.1119度*2^16assign rot[10] = 32'd3648 ; //0.0560度*2^16assign rot[11] = 32'd1856 ; //0.0280度*2^16assign rot[12] = 32'd896 ; //0.0140度*2^16assign rot[13] = 32'd448 ; //0.0070度*2^16assign rot[14] = 32'd256 ; //0.0035度*2^16assign rot[15] = 32'd128 ; //0.0018度*2^16//FSM_parameterlocalparam 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 16)*0.6073; end default :; endcaseenalways @(posedge sys_clk or negedge sys_rst) begin if(!sys_rst) cnt |
