CORDIC算法在Verilog中计算三角函数的理论基础与MATLAB实现

一、数学原理

如图所示，(x2,y2)是(x1,y1)旋转θ角度得到，因此可以有：

经过变换得到：

提出cos(θ)便可以得到伪旋转方程：

伪旋转仅实现了正确的角度旋转，但向量模值变为原来的1/cosθ。
tanθ = 2^(-i)，i为自然数，则上式变为：

只需要提前制作一个表格将N次迭代的数据保存，在Verilog中就可以直接调用：

而计算角度和cos、sin值采取的方法不同，计算角度值采取向量模式的CORDIC算法，后者则采用旋转模式的CORDIC算法，两者的使用和区别见：https://wenku.baidu.com/view/ef46b5040740be1e650e9a4b.html?rec_flag=default&sxts=1568622848231（仅学习使用，原作者：何宾）
上述算法通俗的来讲：
1.计算角度：将任意一个向量先旋转45°，然后减半旋转，再减半旋转，向量模式最终是使yi趋向于零，大于零则逆时针旋转，小于0则顺时针旋转，指导yi无线趋向于零最终在正确的值附近摇摆，迭代次数越高，越精确。
2.计算正弦、余弦：将任意一个向量先旋转45°，然后减半旋转，再减半旋转，旋转模式最终是使zi趋向于零，大于零则逆时针旋转，小于0则顺时针旋转，指导zi无线趋向于零最终在正确的值附近摇摆，迭代次数越高，越精确。

二、MATLAB实现

向量模式计算角度（反正切）：

%% ***********************************************************************************
%     已知坐标，用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;
% 初始化----------------------------------------
N = 16;  %迭代次数
tan_table = 2.^-(0 : N-1);
angle_LUT = atan(tan_table);

K = 1;
for k = 0 : N-1
    K = K*(1/sqrt(1 + 2^(-2*k)));
end
 
x = 3;
y = sqrt(3);
angle_accumulate = 0;
 
% cordic算法计算-------------------------------
if (x==0 && y==0) 
    radian_out = 0;
    amplitude_out = 0;
else  % 先做象限判断，得到相位补偿值
    if (x > 0)
        phase_shift = 0;
    elseif (y < 0)
        phase_shift = -pi;
    else
        phase_shift = pi;
    end
  
    for k = 0 : N-1   % 迭代开始
        x_temp = x;
        y_temp = y;
        angle_accumulate_temp=angle_accumulate;
        if (y < 0)  % d(k)=1，逆时针旋转
            x = x_temp - y_temp*2^(-k);
            y = y_temp + x_temp*2^(-k);
            angle_accumulate = angle_accumulate_temp - angle_LUT(k+1);
        else          % d(k)=-1，顺时针旋转
            x = x_temp + y_temp*2^(-k);
            y = y_temp - x_temp*2^(-k);
            angle_accumulate = angle_accumulate_temp + angle_LUT(k+1);
        end     
        radian_out = angle_accumulate + phase_shift; %弧度输出
    end
    
    amplitude_out = x*K;  %幅值输出
end
    
angle_out = radian_out*180/pi;  %相角输出

旋转模式（计算正弦、余弦）：

%% ***********************************************************************************
%     已知相角theta，计算其正弦和余弦值。基本公式如下：
%     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;
% 初始化----------------------------------------
N = 16;  %迭代次数
tan_table = 2.^-(0 : N-1);
angle_LUT = atan(tan_table);
 
K = 1;
for k = 0 : N-1
    K = K*(1/sqrt(1 + 2^(-2*k)));
end
 
theta = 60;
x = 1;
y = 0;
phase_accumulate = theta/180*pi;  %转化为弧度
 
% cordic算法计算-------------------------------
if (phase_accumulate > pi/2)  % 先做象限判断，得到相位补偿值
    phase_accumulate = phase_accumulate - pi;
    sign_x = -1;
    sign_y = -1;
elseif (phase_accumulate < -pi/2)
    phase_accumulate = phase_accumulate + pi;
    sign_x = -1;
    sign_y = -1;
else
    sign_x = 1;
    sign_y = 1;
end
     
 for k = 0 : N-1   % 迭代开始
        x_temp = x;
        y_temp = y;
        phase_accumulate_temp=phase_accumulate;
        if (phase_accumulate > 0)  % d(k)=1，逆时针旋转
            x = x_temp - y_temp*2^(-k);
            y = y_temp + x_temp*2^(-k);
            phase_accumulate = phase_accumulate_temp - angle_LUT(k+1);
        else                                     % d(k)=-1，顺时针旋转
            x = x_temp + y_temp*2^(-k);
            y = y_temp - x_temp*2^(-k);
            phase_accumulate = phase_accumulate_temp + angle_LUT(k+1);
        end
end
    
cos_out = sign_x*x*K;  %余弦输出
sin_out = sign_y*y*K;   %正弦输出

上述代码原作者：https://blog.csdn.net/longxuekun1992/article/details/52435024