任何適合產(chǎn)品實現(xiàn)的算法，都是將簡易實現(xiàn)作為第一目標(biāo)。CORDIC算法是建立在適應(yīng)性濾波器、FFT、解調(diào)器等眾多應(yīng)用基礎(chǔ)上計算超越函數(shù)的方法。其核心思想是二分逐次逼近。

CORDIC（Coordinate Rotation Digital Computer）算法即坐標(biāo)旋轉(zhuǎn)數(shù)字計算方法，是J.D.Volder1于1959年首次提出，主要用于三角函數(shù)、雙曲線、指數(shù)、對數(shù)的計算。該算法通過基本的加和移位運算代替乘法運算，使得矢量的旋轉(zhuǎn)和定向的計算不再需要三角函數(shù)、乘法、開方、反三角、指數(shù)等函數(shù)。

本文是基于使用Verilog HDL設(shè)計實現(xiàn)Cordic算法，實現(xiàn)正弦、余弦、反正切函數(shù)的實現(xiàn)。將復(fù)雜的運算轉(zhuǎn)化成RTL擅長的加減法和乘法，而乘法運算可以用移位運算代替。Cordic算法有兩種模式，旋轉(zhuǎn)模式和向量模式。可以在圓坐標(biāo)系、線性坐標(biāo)系、雙曲線坐標(biāo)系使用。本文初步實現(xiàn)在圓坐標(biāo)系下的兩種模式的算法實現(xiàn)。

旋轉(zhuǎn)模式，迭代位移算法。假設(shè)有一點P0（x0,y0），經(jīng)過逆時針旋轉(zhuǎn)角度θ，到達點Pm（xm，ym），我們根據(jù)數(shù)學(xué)運算可以得到公式如下:

xm = x0cosθ - y0sinθ

= cosθ(x0 – y0tanθ)

ym = x0sinθ +y0cosθ

=cosθ(x0tanθ+y0)

如果不考慮旋轉(zhuǎn)后的向量模值，只考慮旋轉(zhuǎn)角度，即去掉cosθ，得到如下方程式。這里旋轉(zhuǎn)的角度的正確的，但xm和ym的值增加。cosθ值是小于等于1的，它的倒數(shù)值大于等于1，所以xm和ym模值增大了。去掉cosθ項可以方便我們后面的坐標(biāo)平面旋轉(zhuǎn)的計算。這里稱為偽旋轉(zhuǎn)。

xm = x0 – y0tanθ ym =x0tanθ+y0

Cordic的方法核心就是偽旋轉(zhuǎn)，將旋轉(zhuǎn)角θ細(xì)化成若干個大小固定的角度θi，規(guī)定θi滿足tanθi = 2^-i，通過一系列的迭代旋轉(zhuǎn)，每次旋轉(zhuǎn)θi，i為迭代次數(shù)，規(guī)定∑θi的范圍即旋轉(zhuǎn)角度θ的范圍為[-99.7, 99.7]。如果θ的大于這個范圍則可通過三角運算操作轉(zhuǎn)化到該范圍的角度。

我們通過事先將所有每次旋轉(zhuǎn)的角度計算出來，由于每次旋轉(zhuǎn)的角度是固定的，所以經(jīng)過i次旋轉(zhuǎn)的∑θi可能會超過θ，所以就必須設(shè)置一個方向值di，如果旋轉(zhuǎn)角度之和已經(jīng)小于θ，則di為1，下次旋轉(zhuǎn)繼續(xù)為逆時針旋轉(zhuǎn)，如果旋轉(zhuǎn)角度之和大于θ，則di為-1，下次旋轉(zhuǎn)為順時針。設(shè)置zi+1為旋轉(zhuǎn)剩余角度，zi+1 = z0 – di *zi，z0 = θ，隨著i值的增大，zi+1會趨向于0時，即旋轉(zhuǎn)結(jié)束。di與zi的符號位相同。

采用偽旋轉(zhuǎn)的方法，每次提出一個cosθi，旋轉(zhuǎn)結(jié)束后會產(chǎn)生一個∏cosθi的累乘，一旦我們確定了迭代次數(shù)，∏cosθi就是一個常數(shù)，迭代公式可寫為。這是將cosθi提出、tanθi 替換成 2^-i后的結(jié)果。di與zi的符號位相同。

xi+1 = xi - di * yi * 2^-i yi+1 = yi + di * xi * 2^-i zi+1 = z0 - di * θi

設(shè)迭代i = n - 1，那么旋轉(zhuǎn)n次后得到Pm的坐標(biāo)應(yīng)該為（xn * ∏cosθi, yn * ∏cosθi）。因為每次迭代都會提出一個cosθi，旋轉(zhuǎn)n次后的xn和yn就會少乘一個∏cosθi，所以實際上最終的Pm坐標(biāo)角度近似于(xn * ∏cosθi, yn * ∏cosθi)。

xn * ∏cosθi = x0cosθ - y0sinθ yn * ∏cosθi = y0cosθ + x0sinθ xn = 1/∏cosθi (x0cosθ – y0sinθ) yn = 1/∏cosθi (y0cosθ +x0sinθ)

伸縮因子，KN = 1 / ∏cosθi，已知迭代次數(shù)，我們可以預(yù)先計算KN的值。如下這是使用MATLAB計算出的迭代結(jié)果數(shù)值。

xn =KN (x0cosθ – y0sinθ) yn = KN (y0cosθ +x0sinθ)

從上表可以得出，我們預(yù)先計算出KN的值，然后令x0 = ∏cosθi，y0 = 0，則上述公式可化簡為

xn = cosθ yn = sinθ 即可實現(xiàn)正弦、余弦操作了。旋轉(zhuǎn)模式

Cordic算法旋轉(zhuǎn)模式使用Verilog HDL的實現(xiàn)流程

確定迭代次數(shù)，將每次迭代的角度計算出來，預(yù)先定義為參數(shù)，為了避免浮點運算，將角度值向左移位16位，取整數(shù)部分。

根據(jù)迭代公式進行迭代計算，本設(shè)計取16次迭代，從上表可以看出，當(dāng)?shù)螖?shù)越大時，1/∏cosθi會趨向于一個確定值。如果對結(jié)果精度要求更高，可以設(shè)置更高的迭代次數(shù)，根據(jù)迭代次數(shù)，可以將伸縮因子KN = 1/∏cosθi計算出來。同樣將其左移16位。

xi+1 = xi - di * yi * 2^-i yi+1 = yi + di * xi * 2^-i zi+1 = z0 - di * θi

設(shè)置x0 = ∏cosθi，y0 = 0，則求出x16 = cosθ，y16 = sinθ。

這里需要注意的是，我們在進行迭代運算的時候，將2^-i變成移位運算，對于正余弦來說是有正負(fù)的，所以在一開始定義的時候，就應(yīng)該定義成有符號數(shù)，Verilog中也可以定義有符號數(shù)，最高位表示符號位，定義如下

迭代寄存器定義為有符號數(shù)，那么我們移位運算就不能用>>邏輯右移<<邏輯左移或來移位了，而是用>>>算術(shù)右移和<<<算術(shù)左移。邏輯左移也就相當(dāng)于算數(shù)左移，右邊統(tǒng)一添0?，邏輯右移，左邊統(tǒng)一添0?，算數(shù)右移，左邊添加的數(shù)和符號有關(guān)。

例如1010_1010， []是添加的位 邏輯左移一位：0101_010[0] 算數(shù)左移一位：0101_010[0] 邏輯右移一位：[0]101_0101 算數(shù)右移一位：[1]101_0101

迭代運算采用16級流水線進行運算，最終需要判斷輸出的正余弦值在哪個象限，前面講旋轉(zhuǎn)角度θ的范圍為[-99.7,99.7]，不在這個范圍我們要進行三角運算使其滿足這個范圍，當(dāng)輸入的角度小于90度即可進行計算，當(dāng)輸入角度大于90度小于180度，將輸入角度減去90度并設(shè)定當(dāng)前角度處于第二象限，然后進行計算，當(dāng)輸入角度大于180度小于270度，將輸入的角度減去180度設(shè)置當(dāng)前角度處于第三象限，進行計算，當(dāng)輸入的角度大于270度，減去270設(shè)置當(dāng)前角度處于第四象限，進行計算。象限的設(shè)定通過quarant寄存器實現(xiàn)。

如果角度在第一象限，sin(x) = sin(a),cos(x) = sin(a)最后的結(jié)果x16 = cosθ, y16 = sinθ，這里我想起了那句口訣，一全正，二正弦，三正切，四余弦

如果角度在第二象限，

sin(x) = sin(a+90) = cos(a),

cos(x) = cos(a+90) = -sin(a)

如果角度在第三象限，

sin(x) = sin(a+180) = -sin(a),

cos(x) = cos(a+180) = -cos(a)

如果角度在第四象限，

sin(x) = sin(a+270) = cos(a),

cos(x) = cos(a+270) = -sin(a)

對于正數(shù)，我們直接賦值輸出，負(fù)數(shù)，這里使用有符號數(shù)表示，將其取反加1即可。最終使用modelsim對算法進行仿真，從波形圖上看已經(jīng)初步實現(xiàn)了sin，cos函數(shù)。

向量模式

Cordic算法在向量模式下的計算方法和旋轉(zhuǎn)模式基本上是類似的，設(shè)有一點P0(x0, y0)，經(jīng)過旋轉(zhuǎn)一定角度到與x軸重合，得到點Pm(xm, ym)，即ym = 0。

xm = x0cosθ - y0sinθ

= cosθ(x0 – y0tanθ)

ym = y0cosθ + x0sinθ

= cosθ(y0 +x0tanθ) = 0

我們設(shè)置x0 = x, y0 = y, z0 = 0,迭代次數(shù)為16，經(jīng)過16次迭代后得到zn = θ = arctan(y/x)和坐標(biāo)所代表的向量的模值d = xm = xn * ∏cosθi，di與yi方向相反，即當(dāng)時結(jié)束運算。實現(xiàn)方法為判斷yi的符號位，符號位為1，di為1，符號位為0，di為-1。

xi+1 = xi - di * yi * 2^-i

yi+1 = yi + di * xi * 2^-i

zi+1 = z0 - di * θi

關(guān)于反正切函數(shù)，由于在[-99.7°，99.7°]范圍內(nèi)，所以我們輸入向量P0(x0, y0)時，需要保證其在第一、四象限。

下面是使用MATLAB計算出來的數(shù)據(jù)和FPGA計算出來的數(shù)據(jù)進行比較。

從FPGA計算出的結(jié)果與MATLAB來比較，和實際結(jié)果之間的誤差還是挺小的，畢竟是硬件計算出來的數(shù)據(jù)，向量的誤差就比較大了，如果對于精度比較高的計算，我們可以通過提高迭代次數(shù)來提高精度。

使用ISE進行綜合并下載工程到開發(fā)板上實驗結(jié)果比較。

旋轉(zhuǎn)模式

Sin的理論數(shù)據(jù)

Chipscope抓取到的數(shù)據(jù)

Cos的理論數(shù)據(jù)

Chipscope抓取到的數(shù)據(jù)

向量模式

Arctan(y/x)理論數(shù)值

Chipscope抓取到的數(shù)據(jù)

坐標(biāo)的模值的理論數(shù)據(jù)

責(zé)任編輯：haq