HLS 开发学习（三）CORDIC算法

旋转

旋转矩阵是

$$R(\theta )=\left[\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right]$$

CORDIC使用迭代算法将v转到某个角度。

进行一次旋转是$v_i=R_i\ast v_{i-1}$的矩阵向量乘法，以线性方程表示为

$$x_i=x_{i-1}\cos \theta -y_{i-1}\sin \theta$$
$$y_i=x_{i-1}\sin \theta +y_{i-1}\cos \theta$$

各种状况

90°

比如90°的话

有

$$R(\theta )=\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right]$$

实际上就是

$$\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}-y\\ x\end{array}\right]$$

实际运算可以简单的交换xy的值，并取一个相反数。

45°

如果是45°
有

$$R(45°)=\left[\begin{array}{cc}\sqrt{2}/2& -\sqrt{2}/2\\ \sqrt{2}/2& \sqrt{2}/2\end{array}\right]$$

这样计算效率不高，0,1-1这样的倍数乘法更快。

如果变成

$$\left[\begin{array}{cc}1& -1\\ 1& 1\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}x-y\\ x+y\end{array}\right]$$

旋转没问题，但是放大了根号2倍，（那么，古尔丹，代价是什么呢）

暂时先只考虑角度问题。

$\tan {\theta }_i=2^i$

由于

$$\cos ({\theta }_i)=\frac{1}{1+{\tan }^2({\theta }_i)}$$

$$\sin ({\theta }_i)=\frac{\tan ({\theta }_i)}{1+{\tan }^2({\theta }_i)}$$

旋转矩阵可以变为

$$R_i=\frac{1}{\sqrt{1+{\tan }^2({\theta }_i)}}\left[\begin{array}{cc}1& -\tan {\theta }_i\\ tan{\theta }_i& 1\end{array}\right]$$

如果$\tan {\theta }_i=2^i$，那么运算变为

$$v_i=K_i\left[\begin{array}{cc}1& -2^{-i}\\ 2^{-i}& 1\end{array}\right]\left[\begin{array}{c}{x}_{i-1}\\ y_{j-1}\end{array}\right],K=\frac{1}{\sqrt{1+2^{-2i}}}$$

其中，$2^i$相当于数据向左移动i位，$2^{-i}$相当于数据向右移动i位。
这种移位结构十分省资源。

可以添加值为1或-1的因子表示正转或反转。

$$v_i=K_i\left[\begin{array}{cc}1& -\sigma 2^{-i}\\ \sigma 2^{-i}& 1\end{array}\right]\left[\begin{array}{c}{x}_{i-1}\\ y_{j-1}\end{array}\right],K=\frac{1}{\sqrt{1+2^{-2i}}}$$

比例因子通常在旋转完后补偿

$$K(n)={\Pi }_{i=0}^{n-1}{K}_i,\underset{n\to \infty }{\lim }K(n)=0.6072529350088812561694$$

每次迭代都会旋转${\theta }_i=\arctan 2^{-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;
}

笛卡尔坐标和极坐标转换

两个坐标系的转换关系为

$$\begin{array}{rl}x& =r\cos \theta \\ y& =r\sin \theta \\ r& =\sqrt{x^2+y^2}\\ \theta & =atan2(y,x)\end{array}$$

其中

$$atan2(y,x)=\{\begin{array}{llccc}\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{array}$$

这些复杂的三角函数在硬件中不易实现。

所以可以通过与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<>模板类相似，它们都需要两个整数参数来定义位宽（数据所有位数）和数据整数的位宽。

(可以自行查询定点数和浮点数)

浮点数的优化难以自动实现，一般还是手动调。