旋转
旋转矩阵是
R ( θ ) = [ cos θ − sin θ sin θ cos θ ] R(\theta) = \begin{bmatrix} \cos \theta & -\sin \theta\\\sin\theta&\cos\theta \end{bmatrix} R(θ)=[cosθsinθ−sinθcosθ]
CORDIC使用迭代算法将v转到某个角度。
进行一次旋转是 v i = R i ∗ v i − 1 v_i = R_i * v_{i-1} vi=Ri∗vi−1的矩阵向量乘法,以线性方程表示为
x
i
=
x
i
−
1
cos
θ
−
y
i
−
1
sin
θ
x_i = x_{i-1}\cos\theta - y_{i-1}\sin\theta
xi=xi−1cosθ−yi−1sinθ
y
i
=
x
i
−
1
sin
θ
+
y
i
−
1
cos
θ
y_i = x_{i-1}\sin\theta + y_{i-1}\cos\theta
yi=xi−1sinθ+yi−1cosθ
各种状况
90°
比如90°的话
有
R ( θ ) = [ 0 − 1 1 0 ] R(\theta) = \begin{bmatrix} 0 & -1\\1&0 \end{bmatrix} R(θ)=[01−10]
实际上就是
[ 0 − 1 1 0 ] [ x y ] = [ − y x ] \begin{bmatrix} 0 & -1\\1&0 \end{bmatrix} \begin{bmatrix} x\\y \end{bmatrix} = \begin{bmatrix} -y\\x \end{bmatrix} [01−10][xy]=[−yx]
实际运算可以简单的交换xy的值,并取一个相反数。
45°
如果是45°
有
R ( 45 ° ) = [ 2 / 2 − 2 / 2 2 / 2 2 / 2 ] R(45°) = \begin{bmatrix} \sqrt{2}/2 & -\sqrt{2}/2\\\sqrt{2}/2&\sqrt{2}/2 \end{bmatrix} R(45°)=[2/22/2−2/22/2]
这样计算效率不高,0,1-1这样的倍数乘法更快。
如果变成
[ 1 − 1 1 1 ] [ x y ] = [ x − y x + y ] \begin{bmatrix} 1 & -1\\1&1 \end{bmatrix} \begin{bmatrix} x\\y \end{bmatrix} = \begin{bmatrix} x-y\\x+y \end{bmatrix} [11−11][xy]=[x−yx+y]
旋转没问题,但是放大了根号2倍,(那么,古尔丹,代价是什么呢)
暂时先只考虑角度问题。
tan θ i = 2 i \tan\theta_i = 2^i tanθi=2i
由于
cos ( θ i ) = 1 1 + tan 2 ( θ i ) \cos (\theta_i) =\frac{1}{1+\tan^2(\theta_i)} cos(θi)=1+tan2(θi)1
sin ( θ i ) = tan ( θ i ) 1 + tan 2 ( θ i ) \sin (\theta_i) =\frac{\tan(\theta_i)}{1+\tan^2(\theta_i)} sin(θi)=1+tan2(θi)tan(θi)
旋转矩阵可以变为
R i = 1 1 + tan 2 ( θ i ) [ 1 − tan θ i t a n θ i 1 ] R_i = \frac{1}{\sqrt{1+\tan^2(\theta_i)}}\begin{bmatrix} 1 & -\tan\theta_i\\tan \theta_i&1 \end{bmatrix} Ri=1+tan2(θi)1[1tanθi−tanθi1]
如果 tan θ i = 2 i \tan\theta_i = 2^i tanθi=2i,那么运算变为
v i = K i [ 1 − 2 − i 2 − i 1 ] [ x i − 1 y j − 1 ] , K = 1 1 + 2 − 2 i v_i=K_i\begin{bmatrix} 1 & -2^{-i}\\ 2^{-i} &1 \end{bmatrix}\begin{bmatrix} x_{i-1}\\y_{j-1} \end{bmatrix},K=\frac{1}{\sqrt{1+2^{-2i}}} vi=Ki[12−i−2−i1][xi−1yj−1],K=1+2−2i1
其中,
2
i
2^{i}
2i相当于数据向左移动i位,
2
−
i
2^{-i}
2−i相当于数据向右移动i位。
这种移位结构十分省资源。
可以添加值为1或-1的因子表示正转或反转。
v i = K i [ 1 − σ 2 − i σ 2 − i 1 ] [ x i − 1 y j − 1 ] , K = 1 1 + 2 − 2 i v_i=K_i\begin{bmatrix} 1 & -\sigma2^{-i}\\ \sigma2^{-i} &1 \end{bmatrix}\begin{bmatrix} x_{i-1}\\y_{j-1} \end{bmatrix},K=\frac{1}{\sqrt{1+2^{-2i}}} vi=Ki[1σ2−i−σ2−i1][xi−1yj−1],K=1+2−2i1
比例因子通常在旋转完后补偿
K ( n ) = Π i = 0 n − 1 K i , lim n → ∞ K ( n ) = 0.6072529350088812561694 K(n) = \Pi_{i=0}^{n-1}K_i,\lim\limits_{n\to\infty} K(n) = 0.6072529350088812561694 K(n)=Πi=0n−1Ki,n→∞limK(n)=0.6072529350088812561694
每次迭代都会旋转
θ
i
=
arctan
2
−
i
\theta_i = \arctan 2^{-i}
θi=arctan2−i,
可以提取计算每一个i对于的θ并储存在片上内存,通过查找表寻找这些值。
正弦和余弦的计算
比如计算60°角的正弦和余弦。
由上表可知,我们可以先正旋转两次,得到71.565°,再逆向旋转一次,得到57.529°,再正向一次,逆向一次,得到61.078°,不断的旋转的话,每次旋转的角度会越来越小,角度趋近于60°,增益也趋于稳定。
如果我们初始的矢量成都为1,那么现在的x和y值就是余弦值和正弦值了。
计算的代码如下
//The file cordic.h holds definitions for the data types and constant valuse
#include "cordic.h"
//The cordic_phase array holds the angle for the current rotation
THETA_TYPE cordic_phase[NUM_ITERATIONS] = {
45, 25.56, 14.036, 7.125
3.576, 1.790, 0.895, ...
};
void
cordic(THETA_TYPE theta, COS_SIN_TYPE &s, COS_SIM_TYPE &c)
{
//Set the initial vector that we will rotate
//current_cos = I;current = Q
COS_SIN_TYPE current_cos =0.60735;
COS_SIN_TYPE current_sin = 0.0;
//Factor is the 2^(-L) value
COS_SIN_TYPE factor = 1.0;
//This loop iteratively rotates the initial vector to find the
//sine and cosine value corresponding to the input theta angle
for(int j = 0; j < NUM_ITERATIONS; j++){
//Determine if we are rotating by a positive or negative angle
int sigma = (theta < 0) ? -1 : 1;
//Save the current_cos,so that it can be used in the sine calculation
COS_SIN_TYPE temp_cos = current_cos;
//Perform the rotation
current_cos = current_cos - current_sin * sigma * factor;
current_sin = temp_cos * sigma * factor + current_sin;
//Determine the new theta
theta = theta - sigma * cordic_phase[j];
//Calculata next 2^(-L) value
factor = factor >> 1;
}
//Set the final sine and cosine values
s = current_sin;c = current_cos;
}
笛卡尔坐标和极坐标转换
两个坐标系的转换关系为
x = r cos θ y = r sin θ r = x 2 + y 2 θ = a t a n 2 ( y , x ) \begin{aligned} x &= r\cos\theta\\y&=r\sin\theta\\r&=\sqrt{x^2+y^2}\\ \theta&=atan2(y,x) \end{aligned} xyrθ=rcosθ=rsinθ=x2+y2=atan2(y,x)
其中
a t a n 2 ( y , x ) = { arctan ( y / x ) i f x > 0 arctan ( y / x ) + π i f x < 0 a n d y ≥ 0 arctan ( y / x ) − π i f x < 0 a n d y < 0 π 2 i f x = 0 a n d y > 0 − π 2 i f x = 0 a n d y < 0 u n d e f i n e d i f x = 0 a n d y = 0 atan2 (y,x) =\begin{cases} \arctan(y/x)&if & x >0\\ \arctan(y/x)+\pi&if & x <0&and&y\ge 0\\ \arctan(y/x)-\pi&if & x <0&and&y<0\\ \frac{\pi}{2}&if & x = 0&and&y>0\\ -\frac{\pi}{2}&if & x = 0&and&y<0\\ undefined&if & x =0&and&y=0\\ \end{cases} atan2(y,x)=⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎧arctan(y/x)arctan(y/x)+πarctan(y/x)−π2π−2πundefinedififififififx>0x<0x<0x=0x=0x=0andandandandandy≥0y<0y>0y<0y=0
这些复杂的三角函数在硬件中不易实现。
所以可以通过与cordic法计算sine和cosine类似的方法来计算,也就是说,把这个矢量旋转到正半轴上√。
女少口阿!
cordic中的数字
#include "inttypes.h"
uint16_t a = 0x4000;
uing16_t b = 0x4000;
//Danger! p depends on sizeof(int)
uint32_t p = a * b;
尽管a和b的值可以用16位二进制数表示,而且他们的积(0x10000000)可以精确的用32位二进制数表示,但是依据C99的转换协议在计算开始时会把a和b的值强制转换为int类型,然后计算出一个整型的结果,之后再将处理结果拓展到32位。
如果
#include "inttypes.h"
//4.0 represented with 12 fractional bits
uint16_t a = 0x4000;
//4.0 represented with 12 fractional bits
uint16_t b = 0x4000;
//Danger! p depends on sizeof(int)
uint32_t p = (a*b) >> 12;
或者
#include "inttypes.h"
uint16_t a = 0x4000;
uint16_t b = 0x4000;
//p is assigned to 0x10000000
uint32_t p = (uint32_t) a*(uint32_t) b;
以及
#include "inttypes.h"
//4.0 represented with 12 fractional bits
uint16_t a = 0x4000;
//4.0 represented with 12 fractional bits
uint16_t b = 0x4000;
//Danger! p depends on sizeof(int)
uint32_t p = ((uint32_t) a*(uint32_t) b) >> 12;
最好使用用c++和Vivado HLS模板类apint<>, ap_uint<>, ap fixed<>, ap_ufixed<>来表示任意精度数据。ap_int<>和ap_uint<>模板类需要整数参数来定义他们的位宽。计算函数通常产生结果的位宽足够宽,可以表示正确结果。
ap_fixed<>和ap_ufixed<>模板类相似,它们都需要两个整数参数来定义位宽(数据所有位数)和数据整数的位宽。
(可以自行查询定点数和浮点数)
浮点数的优化难以自动实现,一般还是手动调。