CORDIC算法：用移位与加减实现硬件高效三角函数计算

1. 从“算不了”到“算得快”：CORDIC算法的工程魅力

最近因为一个FPGA项目，需要实现高精度的三角函数和反三角函数运算。在嵌入式领域，尤其是FPGA和低端MCU上，直接调用数学库进行浮点运算往往是奢侈的，资源消耗大、速度慢。于是，我不得不重新捡起那个久闻大名却一直敬而远之的算法——CORDIC。和很多朋友一样，我最初查资料时也是一头雾水：满篇的矩阵、公式推导，看完了也不知道怎么写到代码里，更别提用硬件描述语言（如Verilog）去实现了。那种感觉，就像给你一张复杂的地图，却没告诉你现在站在哪儿。

经过几天的折腾，我终于把这块硬骨头啃了下来。我发现，CORDIC（坐标旋转数字计算机）的精髓，恰恰在于它用一套极其简单、规则的迭代操作，解决了复杂的非线性函数计算问题。它把复杂的乘除、开方、三角函数运算，巧妙地转化为了一系列的移位和加减运算。这对于硬件实现来说，简直是天作之合。今天，我就抛开那些让人望而生畏的纯数学推导，从一个工程师的视角，带你“分分钟”看懂CORDIC的核心思想，并手把手带你走一遍计算流程，让你明白它到底是怎么“转”出正弦余弦值的。

这篇文章适合所有被数学公式吓退的嵌入式工程师、FPGA开发者，以及对高效算法实现感兴趣的朋友。我们不讲高深理论，只关注“怎么用”和“为什么这么用”。

2. CORDIC算法的核心思想：像圆规一样“转”出答案

2.1 问题本质：坐标旋转与三角函数

让我们先回到最根本的问题：我们想求sin(30°)和cos(30°)。在几何上，这等价于什么呢？想象一个单位圆（半径为1），圆上有一个点，初始位置在(1, 0)，也就是0度角的位置。如果我们把这个点逆时针旋转30度，那么它新的坐标(x’, y’)是多少？根据三角函数的定义，这个新点的横坐标x’就是cos(30°)，纵坐标y’就是sin(30°)。

所以，计算三角函数的问题，转化为了旋转一个向量并求其新坐标的问题。CORDIC算法正是抓住了这个几何本质。它的全称“坐标旋转数字计算”也由此而来。但关键来了：计算机不会“连续地”旋转，它只能进行离散的、分步的操作。CORDIC的天才之处在于，它设计了一套特殊的旋转规则。

2.2 核心策略：化整为零的“微旋转”

CORDIC算法的核心策略是：将一个大角度的旋转，分解为许多次预先定义好的、固定小角度的旋转。这些固定的小角度有一个巧妙的选择：它们使得每次旋转的计算变得异常简单。

具体来说，算法预先计算好一系列角度值：θ_i = arctan(2^{-i})。这里的i从0开始递增。我们来看一下前几个值：

• i=0: θ₀ = arctan(2⁰) = arctan(1) = 45°
• i=1: θ₁ = arctan(2⁻¹) = arctan(0.5) ≈ 26.565°
• i=2: θ₂ = arctan(2⁻²) = arctan(0.25) ≈ 14.036°
• i=3: θ₃ = arctan(2⁻³) = arctan(0.125) ≈ 7.125°
• ...

你会发现，这些角度值在迅速减小。算法的目标就是：用这一系列角度θ_i的正向或反向旋转，去逼近我们想要的目标角度（比如30度）。

2.3 旋转的简化：从复杂乘法到简单移位

在平面上旋转一个向量 (x, y) 一个角度 θ，标准的旋转公式是：

x' = x * cosθ - y * sinθ y' = y * cosθ + x * sinθ

这里涉及四次乘法和两次加法，对于硬件并不友好。

CORDIC的第二个神来之笔是：它强制每次旋转时，让旋转角度的正切值等于2的负整数次幂，即tan(θ_i) = 2^{-i}。这意味着sinθ_i ≈ 2^{-i}，cosθ_i ≈ 1（当i稍大时）。如果我们把每次旋转的cosθ_i因子暂时忽略（最后统一补偿），那么旋转公式就简化成了：

x_{i+1} = x_i - d_i * y_i * 2^{-i} y_{i+1} = y_i + d_i * x_i * 2^{-i} z_{i+1} = z_i - d_i * θ_i

其中：

• (x_i, y_i)是当前点的坐标。
• z_i是当前剩余需要旋转的角度（初始值为目标角度）。
• d_i是本次旋转的方向，取+1（逆时针）或-1（顺时针）。它的选择原则很简单：让剩余角度z_i趋近于0。如果z_i > 0，说明还需要逆时针转，则d_i = +1；如果z_i < 0，说明转过头了需要顺时针回调，则d_i = -1。
• 2^{-i}这个乘法操作，在二进制计算机或数字电路中，等价于将y_i或x_i向右移位 i 位！这是整个算法硬件友好的关键。

注意：这里我们忽略了cosθ_i因子，多次旋转累积会引入一个固定的伸缩因子K = Π cosθ_i。当迭代次数足够多时（例如16次），K约等于0.607252935。因此，在最终得到结果(x_n, y_n)后，我们需要将其除以K（或乘以1/K≈1.64676），才能得到正确的坐标值。在“向量模式”（用于计算模长和角度）下，这个因子同样重要。

3. 手算演示：CORDIC如何“转”出sin30°和cos30°

理论说得再多，不如亲手算一遍。我们就以计算sin(30°)和cos(30°)为例，目标角度Z0 = 30°。初始向量设在(1, 0)。我们进行前几次迭代，看看算法是如何逼近的。

初始化:

• x0 = 1
• y0 = 0
• z0 = 30°(剩余角度)
• 预定义角度表:θ0=45°, θ1=26.565°, θ2=14.036°, θ3=7.125°, θ4=3.576°...

迭代过程:

第1次迭代 (i=0, θ0=45°):

• 判断方向：z0 = 30° > 0，需要逆时针转，所以d0 = +1。
• 更新坐标：x1 = x0 - d0 * y0 * 2⁰ = 1 - 1*0*1 = 1y1 = y0 + d0 * x0 * 2⁰ = 0 + 1*1*1 = 1
• 更新剩余角度：z1 = z0 - d0 * θ0 = 30° - 45° = -15°
• 解读：我们试图用45°的大步去逼近30°，结果转过了头，变成了-15°。当前坐标变成了(1,1)，对应角度45°。

第2次迭代 (i=1, θ1=26.565°):

• 判断方向：z1 = -15° < 0，剩余角度为负，说明需要顺时针回调，所以d1 = -1。
• 更新坐标：x2 = x1 - d1 * y1 * 2⁻¹ = 1 - (-1)*1*0.5 = 1 + 0.5 = 1.5y2 = y1 + d1 * x1 * 2⁻¹ = 1 + (-1)*1*0.5 = 1 - 0.5 = 0.5
• 更新剩余角度：z2 = z1 - d1 * θ1 = -15° - (-1)*26.565° = -15° + 26.565° = 11.565°
• 解读：这次我们顺时针回调了26.565°，剩余角度变成了正的11.565°。坐标变成了(1.5, 0.5)。

第3次迭代 (i=2, θ2=14.036°):

• 判断方向：z2 = 11.565° > 0，需要逆时针转，d2 = +1。
• 更新坐标：x3 = x2 - d2 * y2 * 2⁻² = 1.5 - 1*0.5*0.25 = 1.5 - 0.125 = 1.375y3 = y2 + d2 * x2 * 2⁻² = 0.5 + 1*1.5*0.25 = 0.5 + 0.375 = 0.875
• 更新剩余角度：z3 = z2 - d2 * θ2 = 11.565° - 14.036° = -2.471°
• 解读：逆时针转14.036°，稍微又转过了头一点。坐标更新为(1.375, 0.875)。

第4次迭代 (i=3, θ3=7.125°):

• 判断方向：z3 = -2.471° < 0，需要顺时针回调，d3 = -1。
• 更新坐标：x4 = x3 - d3 * y3 * 2⁻³ = 1.375 - (-1)*0.875*0.125 = 1.375 + 0.109375 = 1.484375y4 = y3 + d3 * x3 * 2⁻³ = 0.875 + (-1)*1.375*0.125 = 0.875 - 0.171875 = 0.703125
• 更新剩余角度：z4 = z3 - d3 * θ3 = -2.471° - (-1)*7.125° = -2.471° + 7.125° = 4.654°
• 解读：顺时针回调后，剩余角度又变为正。坐标变为(1.484375, 0.703125)。

我们可以继续迭代下去。每次迭代，旋转的步长θ_i都在减半（近似），剩余角度z_i的绝对值也越来越小，就像用一把刻度越来越精细的尺子去测量一个长度，最终会无限逼近真实值。

最终结果与补偿： 假设我们迭代了N次（比如16次），得到了最终坐标(x_N, y_N)和剩余角度z_N ≈ 0。别忘了我们之前忽略的伸缩因子K。我们最终的单位圆坐标应该是：

cosθ = x_N / K sinθ = y_N / K

对于大量迭代，K趋近于一个常数。如果我们初始向量不是(1,0)，而是(K, 0)，那么经过旋转后，最终的(x_N, y_N)就直接是(cosθ, sinθ)了！这是工程上常用的技巧，称为初始值补偿。

实操心得：在手算或理解阶段，可以忽略K因子，专注于角度逼近的过程。但在实际硬件实现时，必须考虑K因子。通常有两种处理方式：1) 将初始x值设为1/K，这样迭代后输出即为正确值；2) 在迭代结束后，将结果乘以K的倒数。在FPGA中，这个乘法通常可以用一个固定的乘法器或更复杂的移位加法来实现。

4. 从MATLAB仿真到硬件实现思路

理解了原理，我们就能看懂文章开头那段MATLAB代码了。它正是实现了上述的迭代过程。我们重新审视并解读一下关键部分：

N = 30; % 迭代30次 x_i = [1, zeros(1,N-1)]; % 初始化x坐标数组，从1开始 y_i = [0, zeros(1,N-1)]; % 初始化y坐标数组，从0开始 z_i = [pi/6, zeros(1,N-1)]; % 初始化剩余角度数组，目标30度（π/6弧度） for n = 0:N-1 % 注意循环索引从0开始更符合公式中的i if(z_i(n+1) > 0) d = 1; else d = -1; end % 核心迭代公式 x_i(n+2) = x_i(n+1) - ( y_i(n+1) * d * 2^(-n) ); y_i(n+2) = y_i(n+1) + ( x_i(n+1) * d * 2^(-n) ); z_i(n+2) = z_i(n+1) - ( atan(2^(-n)) * d ); % 使用预存的arctan(2^-i)值 end % 计算正弦余弦（这里假设伸缩因子K已通过初始值处理，代码中直接归一化） sin_30 = y_i ./ sqrt(x_i.^2 + y_i.^2); % 实际上对于单位圆，sqrt(x^2+y^2)应约等于K cos_30 = x_i ./ sqrt(x_i.^2 + y_i.^2);

这段代码清晰地展示了算法的流程。运行后，sin_30和cos_30向量的最后一个值会无限逼近0.5和0.8660。

4.1 FPGA硬件实现架构

在FPGA中实现CORDIC，我们追求的是高速、流水线和资源效率。其核心是一个可重复使用的迭代单元（Processing Element, PE）。下面是一个典型的流水线型CORDIC结构思路：

1. 预计算角度表：将arctan(2^{-i})的值预先计算好，量化成固定位宽（如16位、24位）后，存储在FPGA的ROM或作为常量。这是唯一的“表”，非常小。

2. 迭代单元设计：每个PE完成一次迭代操作。

   • 输入：x_i,y_i,z_i,arctan_table[i]。
   • 方向决策：根据z_i的符号位（MSB）决定d_i。d_i = ~sign(z_i)（符号位取反，或直接判断大于/小于0）。
   • 移位器：根据当前迭代次数i，将y_i和x_i进行右移i位操作。在硬件中，这通过固定布线实现，无需桶形移位器。
   • 加减法单元：执行x_i ± (y_i>>i)和y_i ± (x_i>>i)以及z_i ± arctan_table[i]。
   • 输出：x_{i+1},y_{i+1},z_{i+1}。

3. 流水线结构：将N个PE串联起来，每个PE对应一次迭代。数据像流水一样依次通过每个PE。在每个时钟周期，流水线的每一级都在处理不同数据的一次迭代。这样，虽然单个结果需要N个时钟周期才能输出（延迟），但每个时钟周期都能输出一个新的结果（吞吐量高）。

4. 初始值与缩放因子：

   • 旋转模式（求三角函数）：初始化x0 = K(或1/K，取决于设计)，y0 = 0,z0 = 目标角度。输出x_N ≈ cosθ,y_N ≈ sinθ。
   • 向量模式（求模长和角度）：初始化x0 = 输入x,y0 = 输入y,z0 = 0。通过旋转使y趋向于0，最终z_N ≈ arctan(y/x)，x_N ≈ K * sqrt(x^2+y^2)。

注意事项：硬件实现时必须考虑数据位宽和量化误差。随着迭代进行，数据需要保持足够的精度。通常采用符号扩展和定点数表示。迭代次数N决定了精度和资源消耗的平衡，一般16次迭代能达到16位精度，足够多数应用。

4.2 与查找表法的对比

在FPGA中，计算三角函数还有一种常见方法是查找表（LUT）。我们来简单对比一下：

选择建议：如果需要同时计算多种函数（如同时需要幅值和相位），或者FPGA的BRAM资源非常紧张，CORDIC是更优选择。如果只需求某一个固定精度的三角函数，且对单周期延迟有严格要求，查找表可能更简单直接。

5. 常见问题与工程实践陷阱

在实际将CORDIC算法投入使用的过程中，我踩过不少坑。这里总结几个典型问题，希望能帮你绕开这些弯路。

5.1 精度不够：迭代次数怎么选？

这是最常遇到的问题。迭代次数N直接决定了最终精度。理论上，第i次迭代引入的角度误差最大约为θ_i ≈ 2^{-i}(弧度)。经过N次迭代后，最大的角度误差约为2^{-N}弧度。同时，由于忽略cos因子引起的幅度误差也会随着N增加而减小。

经验公式：对于输出位宽为B bit的定点数系统，通常选择N >= B。例如，要获得16位有效精度的输出，迭代16次是常见的起点。你可以通过MATLAB或Python建立定点数模型进行仿真，直接观察不同N值下的输出误差，从而确定满足你系统要求的最小N。

避坑技巧：不要盲目追求高迭代次数。在FPGA中，每增加一级迭代，就增加一级流水线延迟和相应的逻辑资源。在满足精度要求的前提下，选择最小的N。对于通信中的载波同步等应用，12-14次迭代可能就够了。

5.2 结果不对：伸缩因子K忘了吗？

这是我第一次实现时犯的错误。仿真结果总是比预期值大一个固定的倍数（约1.646）。这就是忘了处理伸缩因子K。

解决方案：

1. 初始化补偿（推荐）：在旋转模式下，将初始向量设置为(x0, y0) = (A, 0)，其中A = 1/K ≈ 0.607(如果希望输出范围为±1)，或者A = K ≈ 1.647(如果希望输出最大幅值为±K)。这样迭代结束后，(x_N, y_N)就是正确的正弦余弦值。
2. 后处理乘法：迭代结束后，将结果(x_N, y_N)乘以1/K。这需要一个常数乘法器。

在硬件中，K和1/K都是常数，它们的乘法可以通过移位加法的常系数乘法来高效实现，不一定需要通用的乘法器。

5.3 收敛范围：能算任意角度吗？

基本CORDIC旋转模式的收敛范围是[-99.7°, 99.7°]（约±π/2弧度）。这是因为arctan(2^0)=45°，而所有预定义角度θ_i的和收敛于约99.7°。如果输入角度超出这个范围，算法将无法收敛。

工程处理：

• 角度预处理：利用三角函数的周期性，将所有输入角度θ归算到[-π, π]或[0, 2π]内。
• 象限处理：利用三角函数的对称性，将第二象限(π/2, π]的角度用π - θ转换到第一象限计算，然后调整符号。通常的做法是：
  1. 检查角度是否在[π/2, π)，若是，则θ' = π - θ，最终结果sinθ = sinθ',cosθ = -cosθ'。
  2. 检查角度是否在[π, 3π/2)，若是，则θ' = θ - π，最终结果sinθ = -sinθ',cosθ = -cosθ'。
  3. 检查角度是否在[3π/2, 2π)，若是，则θ' = 2π - θ，最终结果sinθ = -sinθ',cosθ = cosθ'。 这样，所有角度都被映射到了第一象限[0, π/2)进行计算，完美契合CORDIC的收敛范围。

5.4 定点数仿真：软件验证的关键一步

在写RTL代码之前，必须用定点数模型进行充分的软件仿真。MATLAB或Python的fixedpoint库（或手动量化）是你的好朋友。

仿真要点：

1. 确定量化格式：例如，Q2.14格式表示1个符号位，2个整数位，14个小数位。
2. 模拟移位操作：定点数的右移即向下取整或舍入，这会引入截断误差。
3. 监控溢出：迭代过程中数据可能增长，需要足够的整数位宽防止溢出。通常需要比输出精度多几位（保护位）。
4. 建立黄金参考模型：使用浮点运算的双精度结果作为参考，计算定点数模型的误差（如绝对误差、信噪比）。

只有软件定点模型的结果正确且满足误差要求，才能保证RTL实现的功能正确性。

5.5 资源与速度的权衡：迭代型 vs 流水线型

CORDIC在硬件中有两种主要实现方式：

• 迭代型：只使用一套计算单元，通过状态机控制进行N个时钟周期的迭代，完成一次计算。资源省，速度慢（吞吐率低）。
• 流水线型：展开为N级流水线，每级一个PE。资源消耗大，速度快（每个时钟周期输出一个结果）。

选择依据：

• 迭代型：适用于对吞吐率要求不高、面积敏感的低速应用，如传感器数据的慢速处理。
• 流水线型：适用于高速数据流处理，如软件无线电、实时图像处理、通信解调。这是FPGA上最常用的形式，能充分发挥FPGA的并行优势。

在最终决定采用哪种方式之前，一定要明确系统的数据吞吐率要求。