cordic算法详解Word文档格式.docx
《cordic算法详解Word文档格式.docx》由会员分享,可在线阅读,更多相关《cordic算法详解Word文档格式.docx(22页珍藏版)》请在冰点文库上搜索。
1、CORDIC的几何原理介绍
假设在xy坐标系中有一个点P1(x1,y1),将P1点绕原点旋转θ角后得到点P2(x2,y2)。
于是可以得到P1和P2的关系。
x2=x1cosθ–y1sinθ=cosθ(x1–y1tanθ)
y2=y1cosθ+x1sinθ=cosθ(y1+x1tanθ)
以上就是CORDIC的几何原理部分,而我们该如何深入理解这个几何原理的真正含义呢?
从原理中,我们可以知道,当已知一个点P1的坐标,并已知该点P1旋转的角度θ,则可以根据上述公式求得目标点P2的坐标。
然后,麻烦来了,我们需要用FPGA去执行上述运算操作,而FPGA的Verilog语言根本不支持三角函数运算。
因此,我们需要对上述式子进行简化操作,将复杂的运算操作转换为一种单一的“迭代位移”算法。
那么,接下来我们介绍优化算法部分。
2、CORDIC的优化算法介绍
我们先介绍算法的优化原理:
将旋转角θ细化成若干分固定大小的角度θi,并且规定θi满足tanθi=2-i,因此∑θi的值在[-99.7°
,99.7°
]范围内,如果旋转角θ超出此范围,则运用简单的三角运算操作即可(加减π)。
然后我们需要修改几何原理部分的假设,假设在xy坐标系中有一个点P0(x0,y0),将P0点绕原点旋转θ角后得到点Pn(xn,yn)。
于是可以得到P0和Pn的关系。
xn=x0cosθ–y0sinθ=cosθ(x0–y0tanθ)
yn=y0cosθ+x0sinθ=cosθ(y0+x0tanθ)
然后,我们将旋转角θ细化成θi,由于每次的旋转角度θi是固定不变的(因为满足tanθi=2-i),如果一直朝着一个方向旋转则∑θi一定会超过θ(如果θ在[-99.7°
]范围内)。
因此我们需要对θi设定一个方向值di。
如果旋转角已经大于θ,则di为-1,表示下次旋转为顺时针,即向θ靠近;
如果旋转角已经小于θ,则di为1,表示下次旋转为逆时针,即也向θ靠近。
然后我们可以得到每次旋转的角度值diθi,设角度剩余值为zi+1,则有zi+1=zi-diθi,其中z0为θ。
如此随着i的增大,角度剩余值zi+1将会趋近于0,此时运算结束。
(注:
可以发现,di与zi的符号位相同)
第一次旋转θ0,d0为旋转方向:
x1=cosθ0(x0–d0y0tanθ0)
y1=cosθ0(y0+d0x0tanθ0)
第二次旋转θ1,d1为旋转方向:
x2=cosθ1(x1–d1y1tanθ1)=cosθ1cosθ0(x0–d0y0tanθ0–d1y0tanθ1–d1d0x0tanθ1tanθ0)
y2=cosθ1(y1+d1x1tanθ1)=cosθ1cosθ0(y0+d0x0tanθ0+d1x0tanθ1–d1d0y0tanθ1tanθ0)
大家可能已经发现了,在我们旋转的过程中,式子里一直会有tanθi和cosθi,而每次都可以提取出cosθi。
虽然我们的FPGA无法计算tanθi,但我们知道tanθi=2-i,因此可以执行和tanθi效果相同的移位操作2-i来取代tanθi。
而对于cosθi,我们可以事先全部提取出来,然后等待迭代结束之后(角度剩余值zi+1趋近于0,一般当系统设置最大迭代次数为16时zi+1已经很小了),然后计算出∏cosθi的值即可。
总结一下:
迭代公式有三:
xi+1=xi–diyi2-i,提取了cosθi,2-i等效替换了tanθi之后
yi+1=yi+dixi2-i,提取了cosθi,2-i等效替换了tanθi之后
zi+1=zi-diθi
其中i从0开始迭代,假设当i=n-1时,zn趋近于0,迭代结束。
然后对结果乘上∏cosθi(i从0至n-1),于是得到点Pn(xn∏cosθi,yn∏cosθi),此时的点Pn就近似等于之前假设中的点Pn(xn,yn)了,所以此时的点Pn同样满足之前假设得到的公式:
xn∏cosθi=x0cosθ–y0sinθ
yn∏cosθi=y0cosθ+x0sinθ
由于i从0至n-1,所以上式可以转化成下式:
xn=1/∏cosθi(x0cosθ–y0sinθ),(其中i从0至n-1)
yn=1/∏cosθi(y0cosθ+x0sinθ),(其中i从0至n-1)
注意:
上式中的xn,yn是经过迭代后的结果,而不是之前假设中的点Pn(xn,yn)。
了解这点是十分重要的。
根据高中学的三角函数关系,可以知道cosθi=1/[(1+tan2θi)^0.5]=1/[(1+2-2i)^0.5],而1/[(1+2-2i)^0.5]的极值为1,因此我们可以得出一个结论:
当i的次数很大时,∏cosθi的值趋于一个常数。
关于如何计算∏cosθi的代码如下所示:
closeall;
clear;
clc;
%初始化
die=16;
%迭代次数
jiao=zeros(die,1);
%每次旋转的角度
cos_value=zeros(die,1);
%每次旋转的角度的余弦值
K=zeros(die,1);
%余弦值的N元乘积
K_1=zeros(die,1);
%余弦值的N元乘积的倒数
fori=1:
die
a=2^(-(i-1));
jiao(i)=atan(a);
cos_value(i)=cos(jiao(i));
if(i==1)
K(i)=cos_value(i);
K_1(i)=1/K(i);
else
K(i)=K(i-1)*cos_value(i);
end
end
jiao=vpa(rad2deg(jiao)*256,10)
cos_value=vpa(cos_value,10)
K=vpa(K,10)
K_1=vpa(K_1,10)12345678910111213141516171819202122232425
从上表也可以看出,当迭代次数为16,i=15时,cosθi的值已经非常趋近于1了,∏cosθi的值则约等于0.607253,1/∏cosθi为1.64676。
所以当迭代次数等于16时,通过迭代得到的点Pn坐标已经非常接近之前假设中的点Pn。
所以,当迭代次数等于16时,这个式子是成立的。
此时,已知条件有三个x0、y0和θ。
通过16次迭代,我们可以得到xn和yn。
而式中的∏cosθi是个随i变化的值,我们可以预先将其值存入系统中。
然后,我们人为设置x0=∏cosθi,y0=0,则根据等式,xn=cosθ,yn=sinθ。
其中1/∏cosθi的值我们也同样预先存入系统中。
如此,我们就实现了正弦和余弦操作了。
二、CORDIC的具体操作流程介绍
1、CORDIC的旋转模式
由于算法较复杂,一休哥再总结一些具体的操作流程。
1、设置迭代次数为16,则x0=0.607253,y0=0,并输入待计算的角度θ,θ在[-99.7°
]范围内。
2、根据三个迭代公式进行迭代,i从0至15:
xi+1=xi–diyi2-i
yi+1=yi+dixi2-i
注:
z0=θ,di与zi同符号。
3、经过16次迭代计算后,得到的x16和y16分别为cosθ和sinθ。
至此,关于CORDIC的三角函数cosθ和sinθ的计算原理讲解结束。
关于CORDIC算法计算三角函数cosθ和sinθ的MATLAB代码如下所示:
x=zeros(die+1,1);
y=zeros(die+1,1);
z=zeros(die+1,1);
x
(1)=0.607253;
%初始设置
z
(1)=pi/4;
%待求角度θ
%迭代操作
fori=1:
die
ifz(i)>
=0
d=1;
d=-1;
x(i+1)=x(i)-d*y(i)*(2^(-(i-1)));
y(i+1)=y(i)+d*x(i)*(2^(-(i-1)));
z(i+1)=z(i)-d*atan(2^(-(i-1)));
cosa=vpa(x(17),10)
sina=vpa(y(17),10)
c=vpa(z(17),10)123456789101112131415161718192021222324
2、CORDIC的向量模式
讲完了旋转模式后,我们接着讲讲向量模式下的圆坐标系。
在这里,我们需从头来过了,假设在xy坐标系中有一个点P0(x0,y0),将P0点绕原点旋转θ角后得到点Pn(xn,0),θ在[-99.7°
于是可以得到P0和Pn的关系:
yn=y0cosθ+x0sinθ=cosθ(y0+x0tanθ)=0
如何得到Pn(xn,yn)一直是我们的目标。
而此时,我们还是列出那几个等式:
根据三个迭代公式进行迭代,i从0至15:
不过此时我们尝试改变初始条件:
设置迭代次数为16,则x0=x,y0=y,z0=0,di与yi的符号相反。
表示,经过n次旋转,使Pn靠近x轴。
因此,当迭代结束之后,Pn将近似接近x轴,此时yn=0,可知旋转了θ,即zn=θ=arctan(y/x)。
而
因此,可得ycosθ+xsinθ=0,
xn=1/∏cosθi(xcosθ–ysinθ)=1/∏cosθi{[(xcosθ–ysinθ)^2]^(1/2)}
=1/∏cosθi{[x2cos2θ+y2sin2θ–2xysinθcosθ]^(1/2)}
=1/∏cosθi{[x2cos2θ+y2sin2θ+y2cos2θ+x2sin2θ]^(1/2)}
=1/∏cosθi{[x2+y2]^(1/2)}
由上可以知道,我们通过迭代,就算出了反正切函数zn=θ=arctan(y/x),以及向量OP0(x,y)的长度d=xn*∏cosθi。
关于反正切函数,一休哥要多啰嗦几句了,由于θ在[-99.7°
]范围内,所以我们输入向量OP0(x,y)时,需要保证其在第一、四象限。
关于CORDIC算法计算反三角函数arctanθ的MATLAB代码如下所示:
x
(1)=100;
y
(1)=200;
k=0.607253;
ify(i)>
d=vpa(x(17)*k,10)
a=vpa(y(17),10)
c=vpa(rad2deg(z(17)),10)12345678910111213141516171819202122232425
三、CORDIC的旋转模式——Verilog仿真
一休哥在编写CORDIC算法时,采用了16级流水线,仿真效果十分明显。
以下是顶层文件的代码。
为了避免浮点运算,为了满足精度要求,一休哥对每个变量都放大了2^16倍,并且引入了有符号型reg和算术右移。
关于Verilog代码的编写,一休哥已经不想多说了,因为代码是完全符合我之前所讲的CORDIC的原理与MATLAB仿真代码。
相信大家在看完本文的前两个部分之后,对Verilog的理解应该不是难事儿。
moduleCordic_Test
(
CLK_50M,RST_N,
Phase,
Sin,Cos,Error
);
inputCLK_50M;
inputRST_N;
input[31:
0]Phase;
output[31:
0]Sin;
0]Cos;
0]Error;
`definerot032'
d2949120//45度*2^16
`definerot132'
d1740992//26.5651度*2^16
`definerot232'
d919872//14.0362度*2^16
`definerot332'
d466944//7.1250度*2^16
`definerot432'
d234368//3.5763度*2^16
`definerot532'
d117312//1.7899度*2^16
`definerot632'
d58688//0.8952度*2^16
`definerot732'
d29312//0.4476度*2^16
`definerot832'
d14656//0.2238度*2^16
`definerot932'
d7360//0.1119度*2^16
`definerot1032'
d3648//0.0560度*2^16
`definerot1132'
d1856//0.0280度*2^16
`definerot1232'
d896//0.0140度*2^16
`definerot1332'
d448//0.0070度*2^16
`definerot1432'
d256//0.0035度*2^16
`definerot1532'
d128//0.0018度*2^16
parameterPipeline=16;
parameterK=32'
h09b74;
//K=0.607253*2^16,32'
h09b74,
regsigned[31:
0]x0=0,y0=0,z0=0;
0]x1=0,y1=0,z1=0;
0]x2=0,y2=0,z2=0;
0]x3=0,y3=0,z3=0;
0]x4=0,y4=0,z4=0;
0]x5=0,y5=0,z5=0;
0]x6=0,y6=0,z6=0;
0]x7=0,y7=0,z7=0;
0]x8=0,y8=0,z8=0;
0]x9=0,y9=0,z9=0;
0]x10=0,y10=0,z10=0;
0]x11=0,y11=0,z11=0;
0]x12=0,y12=0,z12=0;
0]x13=0,y13=0,z13=0;
0]x14=0,y14=0,z14=0;
0]x15=0,y15=0,z15=0;
0]x16=0,y16=0,z16=0;
reg[1:
0]Quadrant[Pipeline:
0];
always@(posedgeCLK_50MornegedgeRST_N)
begin
if(!
RST_N)
begin
x0<
=1'
b0;
y0<
z0<
=K;
=32'
d0;
=Phase[15:
0]<
<
16;
x1<
y1<
z1<
elseif(z0[31])
=x0+y0;
=y0-x0;
=z0+`rot0;
=x0-y0;
=y0+x0;
=z0-`rot0;
x21);
z21);
z2<
=z1-`rot1;
x32);
z32);
z3<
=z2-`rot2;
x43);
z43);
z4<
=z3-`rot3;
x54);
z54);
z5<
=z4-`rot4;
x65);
z65);
z6<
=z5-`rot5;
x76);
z76);
z7<
=z6-`rot6;
x87);
z87);
z8<
=z7-`rot7;
x98);
z98);
z9<
=z8-`rot8;
always@(posedgeCLK_50Mor