【算法】Cordic算法的原理及matlab/verilog应用

一、前言

单片机或者FPGA等计算能力弱的嵌入式设备进行加减运算还是容易实现，但是想要计算三角函数（sin、cos、tan），甚至双曲线、指数、对数这样复杂的函数，那就需要费些力了。通常这些函数的计算需要通者查找表或近似计算（如泰勒级数逼近）等技术来转换为硬件易于实现的方式。

CORDIC(Coordinate Rotation Digital Computer, 坐标旋转数字计算方法)算法就是一种化繁为简的算法，通过基本的加减和移位运算代替乘法运算，逐渐逼近目标值，得出函数的数值解。

二、Cordic算法理论推导

理论推导参考：CORDIC算法理论详解_cordic算法详解-CSDN博客，这篇博客的推导仔细而全面。

Cordic算法的基石在于一个规律：从tan45°开始，角度每减半，tan值也接近减半。这一规律直接将三角函数运算变成2的幂运算，而这在数字电路中可直接用移位运算来实现。

三、Cordic算法 matlab实现

由于Cordic旋转角度有限（-99.88°，99.88°），因此对于在第二三象限的目标角度，需要做预处理。一般常用的两种方式：1，第二象限逆时针旋转90°到第一象限，第三象限顺时针旋转90°到第四象限，然后按cordic算法参与计算；2，第二象限顺时针旋转180°到第四象限，第三象限顺时针旋转180°到第一象限，然后按cordic算法参与计算。后面分别用这两种方法进行了matlab实现。

3.1 已知相位（角度）求坐标（正弦余弦）

function [sin_theta,cos_theta] = cordic_sincos(theta,n)
% n:iterations
% theta: -180~180
tan_table = 2.^-(0 : n-1);
angle_rad_lut = atan(tan_table);
%angle_deg_lut = rad2deg( atan(tan_table) );
 
k = 1;
for i = 0 : n-1
    k = k*(1/sqrt(1 + 2^(-2*i)));
end
 

x = k;
y = 0;
theta_tar = theta*pi/180;  % to be rad
z=theta_tar;
 
% preprocess
if (theta_tar > pi/2)  
    theta_tar = theta_tar - pi;
    sign_x = -1;
    sign_y = -1;
elseif (theta_tar < -pi/2)
    theta_tar = theta_tar + pi;
    sign_x = -1;
    sign_y = -1;
else
    sign_x = 1;
    sign_y = 1;
end
     
 for i = 0 : n-1   
         if (z > 0)  
             d =1;
         else
             d=-1;
         end
        x_temp = x;
        y_temp =y;
        z_temp = z;
        x= x_temp - d*y_temp*2^(-i);
        y = y_temp + d*x_temp*2^(-i);
        z = z_temp - d*angle_rad_lut(i+1);
        
end
    

    sin_theta = sign_y*y;
    cos_theta = sign_x*x;

end

3.2 已知坐标求相位

function [theta] = cordic_sincos(x,y,n)
%function:know coordinate,calculate theta
% n:iterations
% theta: -180~180

tan_table = 2.^-(0 : n-1);
angle_rad_lut = atan(tan_table);
 

z = 0;
 
% preprocess
if (x<0 && y>0)  
    x_vec = y;
    y_vec = -x;
    theta_shift =-pi/2;
elseif (x<0 && y<0)
    x_vec = -y;
    y_vec = x;
    theta_shift =pi/2;
else
    x_vec = x;
    y_vec = y;
    theta_shift =0;
end
     
 for i = 0 : n-1   
         if (y_vec > 0)  
             d = -1;
         else
             d=1;
         end
        x_temp = x_vec;
        y_temp =y_vec;
        z_temp = z;
        x_vec= x_temp - d*y_temp*2^(-i);
        y_vec = y_temp + d*x_temp*2^(-i);
        z = z_temp + d*angle_rad_lut(i+1);
        
end
    

    theta = rad2deg(-z - theta_shift);


end

四、cordic算法的verilog实现

待实现