CORDIC 原理

CORDIC（Coordinate Rotation Digital Computer）算法即坐标旋转数字计算方法，是J.D.Volder1于1959年首次提出，主要用于三角函数、双曲线、指数、对数的计算。该算法通过基本的加和移位运算代替乘法运算，使得矢量的旋转和定向的计算不再需要三角函数、乘法、开方、反三角、指数等函数。

CORDIC 理论编辑

1.1、坐标旋转数字计算机CORDIC

坐标旋转数字计算机CORDIC(COordinate Rotation DIgital Computer)算法，通过移位和加减运算，能递归计算常用函数值，如Sin，Cos，Sinh，Cosh等函数，由J. Volder于1959年提出，首先用于导航系统，使得矢量的旋转和定向运算不需要做查三角函数表、乘法、开方及反三角函数等复杂运算。J. Walther在1974年用它研究了一种能计算出多种超越函数的统一算法。

1.2、CORDIC原理

如图所示，初始向量（X(0)，Y(0)）旋转θ角度之后得到向量（X1，Y1），此向量有如下关系：

CORDIC算法

X1=X0*cos(θ)-Y0*sin(θ)=cos(θ)(X0-Y0*tan(θ))

Y1=Y0*cos(θ)+X0*sin(θ)=cos(θ)(Y0+X0*tan(θ))

注：θ为待求角

假设初始向量经过N次旋转之后得到新向量，且每次旋转角度δ正切值都为2的倍数，则第i次旋转角度为δ=arctan(2^(-i))，即cosδ=（1/（1+2^(-2i)））^0.5。

容易得到角度θ≈∑S(i)●δ(i)，其中S(i)=1或－1，表示旋转角度的方向，

第i步旋转可以表示为：

X(i+1)=(（1/（1+2^(-2i)））^0.5)●(X(i)-S(i)Y(i)2^(-i))

Y(i+1)=(（1/（1+2^(-2i)））^0.5)●(Y(i)+S(i)X(i)2^(-i))

其中（1/（1+2^(-2i)））^0.5)称为校模因子，当旋转次数一定时，趋于一个常数，Π（1/（1+2^(-2i)））^0.5)≈0.6073

这样，算法每一步就可以简化为：

X(i+1)=0.6073●(X(i)-S(i)Y(i)2^(-i))

Y(i+1)=0.6073●(Y(i)+S(i)X(i)2^(-i))

从而可以看出，对于移动的角度θ，现在只需要硬件加减法器和移位器就可以算出结果。引入Z，表示i次旋转后相位累加的部分和，则：

Z(i+1)=Z(i)-S(i)arctan(2^(-i))

经过n次旋转之后，Z→0，即与目标角重合，即：

X(n)=X1=X0*cos(θ)-Y0*sin(θ)

Y(n)=Y1=Y0*cos(θ)+X0*sin(θ)

1.3、三角函数的计算

以sin/cos计算为例，可利用正/余弦的和角公式递归进行：

cos(a+b) = cos(a)cos(b) – sin(a)sin(b) = cos(a) [cos(b) – tan(a)sin(b)]

sin(a+b) = sin(a)cos(b) + cos(a)sin(b) = cos(a) [tan(a)cos(b) +sin(b)]

取a=arctan(2^-k), 即tan(a)=2^-k, 则cos(b) – tan(a)sin(b) 可通过移位和减法来实现。

如果角度z可以表示为z = s0 arctan(2^0) + s1 arctan(2^-1) + ... + sn arctan(2^-n), 其中s0, s1, ..., sn取+1或-1(+1可以理解为逆时针转角，即加上一个角度； -1则相反) ，那么角度z的sin/cos计算可以通过一系列的移位和加减运算来实现。注意到cos(sk arctan(2^-k))=cos(arctan(2^-k)) 与转角方向无关。此外，z应取第一项限角度(收敛域)，对於其他项限角度，可由其第一项限对应角度变换得到。

相类似地，sinh/cosh的计算利用以下公式：

cosh(a+b) = cosh(a)cosh(b) + sinh(a)sinh(b) = cosh(a) [cosh(b) + tanh(a)sinh(b)]

sinh(a+b) = sinh(a)cosh(b) + cosh(a)sinh(b) = cosh(a) [tanh(a)cosh(b) + sinh(b)]

取a=arctanh(2^-k), 即tanh(a)=2^-k, 则cosh(b) + tanh(a)sinh(b) 可通过移位和减法来实现。如果参数z可以表示为z = s1 arctanh(2^-1) + s2 arctanh(2^-2) + ... + sn arctanh(2^-n), 其中s1, s2, ..., sn取+1或-1 ，那么z的sinh/cosh计算可以通过一系列的移位和加减运算来实现。

z应取[-ln2, ln2]范围内的值，否则应先预处理 z = z’– pln2, 求得cosh(z’)/sinh(z’)的值，则

cosh(z) = cosh(z’)cosh(pln2) + sinh(z’)sinh(pln2) = ½[cosh(z’) + sinh(z’)]2^p + ½[cosh(z’) – sinh(z’)]2^-p

sinh(z) = sinh(z’)cosh(pln2) + cosh(z’)sinh(pln2) = ½[cosh(z’) + sinh(z’)]2^p + ½[sinh(z’) – cosh(z’)]2^-p 。

sin/cos和sinh/cosh的计算是CORDIC算法的两个特例，CORDIC算法可描述如下：

给定初始值x(0), y(0), z(0),

x(k+1) = x(k) – ms(k)y(k)2^-q(m,k), y(k+1) = y(k) + s(k)x(k) 2^-q(m,k), z(k+1) = z(k) – s(k)d(k),

其中m表示模式，q(m,k) 为移位序列，s(k) 取+1或-1表示旋转方向，d(k) 为递进角度。

1.4、CORDIC的MATLAB 运算

以cos(a)/sin(a)计算为例，m = 1, x(0) = 1, y(0) = 0, z(0) = a, s(k) = sign(z(k)),移位序列q(1,k): 0, 1, 2, ..., 递进角度为d(k)=arctan(2^-q(1,k)) 。

下面是实现的Matlab程序（保存为m文件）：

function [sin,cos] =cordic(angle);

% 初始化

x = 1;

y = 0;

z = angle;

a = 0;

d = 1;

k = 0.6073; %K 增益

x = k*x;

while a<100 %此处不能判断d的符号控制循环,会死循环.应用a次数控制.

if z >=0

d=1;

else

d=-1;

end

%迭代

xNew=x;

x=xNew-(y*d*(1/2^a));

y=y+(xNew*d*(1/2^a));

z=z-(d*(atan(1/2^a)));

a=a+1;

end

cos= x

sin= y

k值为cos(arctan(1)), cos(arctan(2^-1)), ..., cos(arctan(2^-K) 的连积值，收敛为0.6072529350。

以cosh(a)/sinh(a)计算为例，m = -1, x(0) = 1, y(0) = 0, z(0) = a, s(k) = sign(z(k)),移位序列q(-1,k): 1, 2, 3, 4, 4, 5, ... (3n+1重复两次以保证收敛， 4, 13, 40, ...), 递进角度为d(k)=arctanh(2^-q(-1,k)) 。

通过对初始值和旋转方向s(k) 的选择，模式m=0可以计算乘法和除法； 模式m=1可以计算sin/cos/arcsin/arccos/arctan; 模式m=-1可直接计算sinh/cosh/exp/arctanh/ln/sqrt, 间接计算arcsinh/arccosh，参见^[1]。

MATLAB 代码实现2：

m=1;

u=1;

K=1.6468;

a=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ];

y0=0;

x0=1/K;

angle = 30;

z0=angle*pi/180;

for i=1:length(a)

x=x0-m*u*y0*2^(-a(i));

y=y0+u*x0*2^(-a(i));

z=z0-u*atan(2^(-a(i)));

x0=x;

y0=y;

z0=z;

u=sign(z0);

end

cosine = x;

sine = y;

MATLAB 代码实现3：

function [cordic_x,cordic_y,cordic_z]=cordic(x,y,z,m,rotation,iterate)

%% [cordic_x,cordic_y,cordic_z]=cordic(x,y,z,m,rotation,iterate)

%% iterate: iterate times

%% rotation: set to 1 with rotation mode; otherwise set to 0 with vector mode.

%% cordic Algorithm for rotation:Z->0

%% m = 1 :(x,y,z)=(1,0,z); x = cos(z) ;y =sin(z)

%% m = 0 :(x,y,z)=(x,0,z);y = x*z

%% m = -1 :(x,y,z)=(1,0,z); x = cosh(z) ; y = sinh(z) ;

%% exp(z) = cosh(z) + sinh(z)

%%-------------------------------------------------------------------------

%% cordic Algorithm for Vector:y->0

%% m = 1 :(x,y,z)=(x,y,0); z = atan(y/x) ; x = (x^2+y^2)^(1/2)/scale

%% m = 0 :(x,y,z)=(x,y,0); z = y/z

%% m = -1 : (y<x),(x,y,z)=(x,y,0); z =atanh(y/x) ; x = (x^2-y^2)^(1/2)/scale

if(m==-1) %% hyperbolic

k_initial = 1; %% in hyperbolic mode , the k's initial value should be "1"

scale = 1.2075 ;

elseif(m==1) %% circyle

k_initial = 0;

scale = 0.6073;

else %% line

k_initial = 0;

scale = 1;

end

x=x*scale;

for k = k_initial : iterate

if(m==-1)&&(mod(k,3)==1)&&( k~=1) %% in hyperbolic mode ,should iterate twice while k == 3N+1

for i = 1:2

if(m == 1)

u=atan(2^(-k))*180/pi;

elseif(m == 0)

u=2^(-k);

elseif(m == -1)

u=atanh(2^(-k));

end

xpre=x;

ypre=y;

zpre=z;

upre=u;

if(rotation) d = sign(zpre);

else d = -sign(ypre);

end

x=xpre-m*d*ypre*2^(-k);

y=ypre+d*xpre*2^(-k);

z=zpre-d*upre;

end

else

if(m == 1)

u=atan(2^(-k))*180/pi;

elseif(m == 0)

u=2^(-k);

elseif(m == -1)

u=atanh(2^(-k));

end

xpre=x;

ypre=y;

zpre=z;

upre=u;

if(rotation) d = sign(zpre);

else d = -sign(ypre);

end

x=xpre-m*d*ypre*2^(-k);

y=ypre+d*xpre*2^(-k);

z=zpre-d*upre;

end

end

cordic_x=x;

cordic_y=y;

cordic_z=z;

　　

现实意义编辑

由于具有频率精度高、转换时间短、频谱纯度高以及频率相位易编程等特点，数控振荡器(NCO)被广泛应用于软件无线电数字上、下变频以及各种频率和相位数字调制解调系统中。NCO传统的实现方法主要有查表法、多项式展开法或近似法，但这些方法在速度、精度、资源方面难以兼顾。而采用CORDIC算法来实现超函数时，则无需使用乘法器，它只需要一个最小的查找表(LUT)，利用简单的移位和相加运算，即可产生高精度的正余弦波形，尤其适合于FPGA的实现。

数字控制振荡器（NCO，numerical controlled oscillator）是软件无线电、直接数据频率合成器（DDS，Direct digital synthesizer）、快速傅立叶变换（FFT，Fast Fourier Transform）等的重要组成部分，同时也是决定其性能的主要因素之一，随着芯片集成度的提高、在信号处理、数字通信领域、调制解调、变频调速、制导控制、电力电子等方面得到越来越广泛的应用。