FPGA算法學習(1) -- Cordic（圓周系統之向量模式）

旋轉模式用來解決三角函數，實現極坐標到直角坐標的轉換，基礎理論請參考Cordic算法——圓周系統之旋轉模式。那么，向量模式則用來解決反三角函數的問題，體現的應用主要是直角坐標向極坐標轉換，即已知一點的直角坐標(x,y)，求其極坐標(α,γ)，實際上是求arctan(y/x)。

旋轉模式下，每次迭代使z趨近於α（α-z趨近於0），而向量模式下，則使y趨近於0，這一點很好理解，即從坐標位置，旋轉到x正半軸，一共旋轉了多少角度，則該角度即為α，從而知道了極角。

如圖所示，在單位圓上，向量OP與X軸的正半軸夾角為α，故P點的坐標可表示為

根據開頭描述，我們需要轉動向量OP，先順時針旋轉θ角至向量OQ，Q點的坐標可表示為

這里定義θ為目標旋轉角度。根據三角函數公式可將上式展開為

現在已經有點 Cordic 算法的樣子了，但是我們看到每次旋轉都要計算 4 次浮點數的乘法運算，運算量還是太大了。還需要進一步的改進，改進的切入點當然還是坐標變換的過程。

將式(1.1)代入到式(1.3)中可得

用矩陣形式表示為：

旋轉了i次以后，可以得到：

最終需將y_Q_i+1轉為0，先按45°的二分法查找來解釋過程，用C語言描述過程為：

#include <stdio.h>
#include <stdlib.h>

double cordic_v(double x, double y);

int main(viod)
{
	double alfa = cordic_v(120.0,200.0);		//直角坐標（x,y）
	printf("\n 極角為 = %f \n",alfa);
	return 0;
}
double cordic_v(double x, double y)
{
	const double sine[] = {0.7071067811865,0.3826834323651,0.1950903220161,
						0.09801714032956,0.04906767432742,0.02454122852291,0.01227153828572,
						0.006135884649154,0.003067956762966,0.001533980186285,
						7.669903187427045e-4,3.834951875713956e-4,1.917475973107033e-4,
						9.587379909597735e-5,4.793689960306688e-5,2.396844980841822e-5
						};

	const double cosine[] = {0.7071067811865,0.9238795325113,0.9807852804032,0.9951847266722,
							0.9987954562052,0.9996988186962,0.9999247018391,0.9999811752826,0.9999952938096,
							0.9999988234517,0.9999997058629,0.9999999264657,0.9999999816164,0.9999999954041,
							0.999999998851,0.9999999997128
						};
	int i = 0;
	double x_new, y_new;
	double angleSum = 0.0;
	double angle = 45.0;		//第一次旋轉角度為45°
	for( i=0; i<15;i++)
	{
		if(y > 0)
		{
			x_new = x * cosine[i] + y * sine[i];
			y_new = y * cosine[i] - x * sine[i];
			x = x_new;
			y = y_new;
			angleSum += angle;
		}
		
		else
		{
			x_new = x * cosine[i] - y * sine[i];
			y_new = y * cosine[i] + x * sine[i];
			x = x_new;
			y = y_new;
			angleSum -= angle;
		}
		printf("旋轉次數: i = %2d 旋轉角度 = %10f, 累計旋轉角度 = %10f\n", i+1, angle,angleSum );
		angle /= 2;
	}
	return angleSum;
}

經過旋轉模式的推導，向量模式的偽旋轉公式，可表示為

C語言描述過程，如下：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double cordic_v(double x, double y);
double r = 0.0;	//定義一個模長全局變量		

int main(viod)
{
	double alfa = cordic_v(120.0,200.0);		//直角坐標（x,y）
	printf("\n極角 = %5f, 模長 = %5f\n",alfa,r);
	return 0;
}
double cordic_v(double x, double y)
{
	const double theta[] = { 45.0, 26.56505118, 14.03624347, 7.125016349, 
						3.576334375, 1.789910608, 0.8951737102, 0.4476141709, 
						0.2238105004, 0.1119056771, 0.05595289189, 0.02797645262, 
						0.01398822714, 0.006994113675, 0.003497056851, 0.001748528427
						}; //旋轉角度
	int i = 0;
	double x_new, y_new;
	double angleSum = 0.0;
	r = sqrt(x*x+y*y);
	for( i=0; i<16;i++)
	{
		if(y > 0)
		{
			x_new = x + y/(1<<i);
			y_new = y - x/(1<<i);
			x = x_new;
			y = y_new;
			angleSum += theta[i];
		}
		
		else
		{
			x_new = x - y/(1<<i);
			y_new = y + x/(1<<i);
			x = x_new;
			y = y_new;
			angleSum -= theta[i];
		}
		printf("旋轉次數: i = %2d 旋轉角度 = %10f, 累計旋轉角度 = %10f, y = %5f\n", i+1,theta[i],angleSum,y );
	}
	return angleSum;
}

同樣，向量模式的cordic算法適用於第一、四象限的坐標變換，在第二、三象限的坐標需要進行預處理。

參考

• 三角函數計算，Cordic 算法入門--Ivan

• 《基於FPGA的數字信號處理（第2版）》——高亞軍著