作者：桂。
時間：2017-08-14 19:22:26
鏈接：http://www.cnblogs.com/xingshansi/p/7359940.html
前言
CORDIC算法常用來求解信號的幅度與相位，它的優勢在於借助：移位寄存器+加法器/減法器便可以實現求解，而無需乘法器。大大簡化了運算。本文圍繞CORDIC整理用到的知識，先做個引子，不定期更新。
一、CORDIC算法
CORDIC(Coordinate Rotation Digital Computer) 算法由Volder於1959年提出，該算法利用移位和加減運算，實現常用三角函數計算。
A-二分查找
現在給定一坐標點（X = 100,Y = 200）可知theta = 63.435°。
這里只討論第一象限，當四個象限同時討論時，只需要根據x,y的符號位給定一個標志位，對theta進行補償即可。CORDIC的核心思想：將（X,Y）旋轉一定的角度，當縱坐標為0時，旋轉的角度就是theta。順時針旋轉theta后的坐標為（x',y'）：
如何旋轉呢，可以借鑒二分查找法的思想。我們知道θ的范圍是0到90度。那么就先旋轉45度試試。這里為什么選擇45°呢？個人的理解：
tan(x/y) = 90° - tan(y/x)
可見45°是一個分界線。
旋轉之后縱坐標為70.71，還是大於0，說明旋轉的度數不夠，接着再旋轉22.5度（45度二分）
這時總共旋轉了45+22.5=67.5度。結果縱坐標變為了負數，說明θ<67.5度，這時就要往回轉，還是二分查找法的思想，這次轉-11.25度（回轉為負角度）。
依次類推，直到迭代的結果滿足要求為止。什么時候滿足要求呢？這個需要根據項目要求的精度確定，例如誤差＜0.5°，即45°/2(n-1)＜0.5°，二分次數n最小取8。可見只需要存取8個特定角度的正弦+8個特定角度的余弦，便可以實現角度的估計，這少量的數值可以提前寫入存儲器，運行時直接調用。
B-角度運算的初步改進
從旋轉矩陣可以看出，每二分查找一次，就需要4次乘法運算。為了減小運算量，對旋轉矩陣進一步處理：
求解相位時，cos(theta)只影響x,y的比例關系，而比例並不影響theta的取值，因此可以忽略：
這樣便將每次二分的4次乘法運算縮減為2次。
C-角度運算的進一步改進
注意到第一次循環時，tan(45)=1，所以第一次循環實際上是不需要乘法運算的。第二次運算tan(22.5)=0.4142135623731,是個很不整的小數，如果參與乘法的小數比較有規律/比較整，則乘法運算會簡化一些，因此這里轉化思路：不再嚴格的二分區間查找，而是非均勻分割角度區間，使得tan（theta）盡可能整→這樣一來，乘法不再需要乘法器，而是移位便可以實現。
tan(45)= 1 tan(26.565051177078)= 1/2 tan(14.0362434679265)= 1/4 tan(7.1250163489018)= 1/8 tan(3.57633437499735)= 1/16 tan(1.78991060824607)= 1/32 tan(0.8951737102111)= 1/64 tan(0.4476141708606)= 1/128 tan(0.2238105003685)= 1/256 ...
MATLAB仿真（迭代次數由項目指標給定的精度確定）：
function theta = cordic(x,y)
%僅以第一象限為例
angle = atand(2.^([0:-1:-15]));
len = length(angle);
theta = 0;
for i = 0:len-1
if (y>0)
x_new = x+y*2^-i;%此處移位寄存器實現
y_new = y-x*2^-i;%此處移位寄存器實現
x = x_new;
y = y_new;
theta = theta+angle(i+1);
else
x_new = x-y*2^-i;%此處移位寄存器實現
y_new = y+x*2^-i;%此處移位寄存器實現
x = x_new;
y = y_new;
theta = theta-angle(i+1);
end
end
至此便完成了CORDIC的整體思路。其他函數可仿照該思路類推。
CORDIC求解幅值的思路：
可近似處理：
具體參考：《Design of Jacobi EVD processor based on CORDIC for DOA estimation with MUSIC algorithm》
二、求解三角思路延伸
以arctan為例，迭代次數隨着theta精度的增加而增加，迭代次數過多會導致時間延遲過大，因此可以考慮函數逼近的思路。常用的有泰勒逼近和切比雪夫逼近。
因為tan(45°) = 1，且tan(x/y) = 90°-tan(y/x)，因此arctan(x)只需要討論x屬於[0,1]的情況，這一點剛好滿足切比雪夫近似。
A-泰勒逼近
即高等數學里的泰勒展開：
對應的指令：
syms x y = atan(x); p50=taylor(y,'Order',3,'ExpansionPoint',0.5);
Order：指定階數，expansionpoint：指定展開的位置。
求解：
double(subs(p50,x,i))；%i為指定位置的數值
B-切比雪夫逼近
切比雪夫多項式可表示為：
設
，可以求得以下遞推關系：
函數f(x)可以近似表達為傅里葉級數：
其中：
這就完成了切比雪夫近似的思想。
給出代碼實現：
function f = Chebyshev(y,k,x0)
%用切比雪夫多項式逼近已知函數
%已知函數：y（類型：表達式）
%逼近已知函數所需項數：k
%逼近點的x坐標：x0
%求得的切比雪夫逼近多項式或在x0處的逼近值：f
syms t;
T(1:k+1) = t;
T(1) = 1;
T(2) = t;
c(1:k+1) = 0.0;
c(1)=quad(matlabFunction(y(t)*T(1)/sqrt(1-t^2)),-1,1)/pi;
c(2)=2*quad(matlabFunction(y(t)*T(2)/sqrt(1-t^2)),-1,1)/pi;
f = c(1)+c(2)*t;
for i=3:k+1
T(i) = 2*t*T(i-1)-T(i-2);
c(i) = 2*quad(matlabFunction(y(t)*T(i)/sqrt(1-t^2)),-1,1)/pi;
f = f + c(i)*T(i);
f = vpa(f,6);
if(i==k+1)
if(nargin == 3)
f = subs(f,'t',x0);
else
f = vpa(expand(f),6);
end
end
end
函數中vpa用來控制精度，例如：
digits(3);%控制精度 vpa(1/3);
結果便是0.333，而不是0.33333333...
調用形式：
Chebyshev(y,3,i);%i為對應橫坐標
這里對於taylor、chebv展開，都按3階處理（其他階數類似），得出逼近誤差：
可以看出chebv的逼近更加穩健。
三、CORDIC仿真
分兩個思路整理，首先是調用IP核，其次是不使用IP核。
A-調用IP核
直接看使用手冊即可。
B-不使用IP核
//CORDIC implementation for sine and cosine for Final Project
//Claire Barnes
module CORDIC(clock, cosine, sine, x_start, y_start, angle);
parameter width = 16;
// Inputs
input clock;
input signed [width-1:0] x_start,y_start;
input signed [31:0] angle;
// Outputs
output signed [width-1:0] sine, cosine;
// Generate table of atan values
wire signed [31:0] atan_table [0:30];
assign atan_table[00] = 'b00100000000000000000000000000000; // 45.000 degrees -> atan(2^0)
assign atan_table[01] = 'b00010010111001000000010100011101; // 26.565 degrees -> atan(2^-1)
assign atan_table[02] = 'b00001001111110110011100001011011; // 14.036 degrees -> atan(2^-2)
assign atan_table[03] = 'b00000101000100010001000111010100; // atan(2^-3)
assign atan_table[04] = 'b00000010100010110000110101000011;
assign atan_table[05] = 'b00000001010001011101011111100001;
assign atan_table[06] = 'b00000000101000101111011000011110;
assign atan_table[07] = 'b00000000010100010111110001010101;
assign atan_table[08] = 'b00000000001010001011111001010011;
assign atan_table[09] = 'b00000000000101000101111100101110;
assign atan_table[10] = 'b00000000000010100010111110011000;
assign atan_table[11] = 'b00000000000001010001011111001100;
assign atan_table[12] = 'b00000000000000101000101111100110;
assign atan_table[13] = 'b00000000000000010100010111110011;
assign atan_table[14] = 'b00000000000000001010001011111001;
assign atan_table[15] = 'b00000000000000000101000101111100;
assign atan_table[16] = 'b00000000000000000010100010111110;
assign atan_table[17] = 'b00000000000000000001010001011111;
assign atan_table[18] = 'b00000000000000000000101000101111;
assign atan_table[19] = 'b00000000000000000000010100010111;
assign atan_table[20] = 'b00000000000000000000001010001011;
assign atan_table[21] = 'b00000000000000000000000101000101;
assign atan_table[22] = 'b00000000000000000000000010100010;
assign atan_table[23] = 'b00000000000000000000000001010001;
assign atan_table[24] = 'b00000000000000000000000000101000;
assign atan_table[25] = 'b00000000000000000000000000010100;
assign atan_table[26] = 'b00000000000000000000000000001010;
assign atan_table[27] = 'b00000000000000000000000000000101;
assign atan_table[28] = 'b00000000000000000000000000000010;
assign atan_table[29] = 'b00000000000000000000000000000001;
assign atan_table[30] = 'b00000000000000000000000000000000;
reg signed [width:0] x [0:width-1];
reg signed [width:0] y [0:width-1];
reg signed [31:0] z [0:width-1];
// make sure rotation angle is in -pi/2 to pi/2 range
wire [1:0] quadrant;
assign quadrant = angle[31:30];
always @(posedge clock)
begin // make sure the rotation angle is in the -pi/2 to pi/2 range
case(quadrant)
2'b00,
2'b11: // no changes needed for these quadrants
begin
x[0] <= x_start;
y[0] <= y_start;
z[0] <= angle;
end
2'b01:
begin
x[0] <= -y_start;
y[0] <= x_start;
z[0] <= {2'b00,angle[29:0]}; // subtract pi/2 for angle in this quadrant
end
2'b10:
begin
x[0] <= y_start;
y[0] <= -x_start;
z[0] <= {2'b11,angle[29:0]}; // add pi/2 to angles in this quadrant
end
endcase
end
// run through iterations
genvar i;
generate
for (i=0; i < (width-1); i=i+1)
begin: xyz
wire z_sign;
wire signed [width:0] x_shr, y_shr;
assign x_shr = x[i] >>> i; // signed shift right
assign y_shr = y[i] >>> i;
//the sign of the current rotation angle
assign z_sign = z[i][31];
always @(posedge clock)
begin
// add/subtract shifted data
x[i+1] <= z_sign ? x[i] + y_shr : x[i] - y_shr;
y[i+1] <= z_sign ? y[i] - x_shr : y[i] + x_shr;
z[i+1] <= z_sign ? z[i] + atan_table[i] : z[i] - atan_table[i];
end
end
endgenerate
// assign output
assign cosine = x[width-1];
assign sine = y[width-1];
endmodule
參考：
- http://blog.csdn.net/liyuanbhu/article/details/8458769