1. 为什么要在FPGA里用CORDIC算三角函数做数字信号处理的朋友大概率都遇到过这个场景需要实时算一个角度的正弦和余弦值但手头的FPGA里没有硬核浮点单元也不想为了几个三角函数就去例化一个占资源的浮点IP。这时候CORDIC算法就是最顺手的一把刀。CORDIC的全称是Coordinate Rotation Digital Computer翻译过来叫坐标旋转数字计算机。名字听着唬人核心思想其实特别朴素用一系列固定角度的旋转去逼近任意角度。就像你手里只有几把固定角度的尺子但通过反复叠加和翻转最终能拼出你想要的任何角度。它最大的优势是只靠移位和加法就能完成三角函数、反三角函数、双曲函数、向量模长等运算完全不需要乘法器更不需要浮点单元。对于资源紧张的FPGA来说这简直是量身定做的方案。这次我拿EGo1板卡做上板验证。EGo1是Xilinx Artix-7系列的一款教学开发板芯片型号XC7A35T逻辑资源够用板载了足够多的开关、按键和LED非常适合做这种算法验证类的小项目。整个项目的目标很明确用Verilog手写一个CORDIC旋转模式的计算核输入一个角度输出对应的sin和cos值通过板上的数码管或者LED把结果展示出来。这篇文章我会从算法原理、硬件架构设计、Verilog实现细节、上板调试过程、以及踩过的坑这几个维度完整展开。不管你是刚接触FPGA的新手还是已经做过几个项目想补一补算法实现经验的老手应该都能从中找到有用的东西。尤其是那些正在准备FPGA项目实战、嵌赛或者课程设计的同学这个项目可以直接拿来当参考模板。2. CORDIC旋转模式的算法原理拆解2.1 从几何旋转到迭代逼近CORDIC旋转模式的数学基础是二维平面上的坐标旋转。假设平面上有一个点 $(x, y)$把它绕原点旋转角度 $\theta$得到新坐标 $(x, y)$用矩阵表示就是$$ \begin{bmatrix} x \ y \end{bmatrix} \begin{bmatrix} \cos\theta -\sin\theta \ \sin\theta \cos\theta \end{bmatrix} \begin{bmatrix} x \ y \end{bmatrix} $$这个公式本身没什么问题但问题是它需要算 $\cos\theta$ 和 $\sin\theta$——我们本来就是要算这个这就成了鸡生蛋蛋生鸡的死循环。CORDIC的巧妙之处在于它不一次性旋转 $\theta$而是把 $\theta$ 拆成一系列预先定好的小角度 $\theta_i$ 的加减组合。每个 $\theta_i$ 满足 $\tan\theta_i 2^{-i}$也就是说 $\theta_i \arctan(2^{-i})$。这样一来旋转矩阵里的 $\cos\theta_i$ 和 $\sin\theta_i$ 就可以用移位和加法来近似因为乘以 $2^{-i}$ 在硬件里就是右移 $i$ 位。把旋转矩阵里的 $\cos\theta_i$ 提出来迭代公式变成$$ \begin{aligned} x_{i1} x_i - d_i \cdot y_i \cdot 2^{-i} \ y_{i1} y_i d_i \cdot x_i \cdot 2^{-i} \ z_{i1} z_i - d_i \cdot \arctan(2^{-i}) \end{aligned} $$其中 $d_i \in {1, -1}$ 是旋转方向。旋转模式的目标是让 $z$ 最终趋近于0也就是把累积的角度误差消掉。$d_i$ 的取值规则很简单如果当前 $z_i 0$说明还需要往正方向转$d_i 1$如果 $z_i 0$说明转多了$d_i -1$。2.2 增益因子K的补偿问题这里有一个容易被忽略的细节每次迭代都会引入一个 $\cos\theta_i$ 的缩放因子。虽然我们把它提出去了但它并没有消失。$n$ 次迭代之后总的缩放因子是$$ K \prod_{i0}^{n-1} \cos(\arctan(2^{-i})) \prod_{i0}^{n-1} \frac{1}{\sqrt{1 2^{-2i}}} $$当迭代次数趋于无穷时$K \approx 0.607252935$。也就是说如果不做补偿最终算出来的 $(x_n, y_n)$ 会比真实值大 $1/K \approx 1.64676$ 倍。处理方式有两种一是在迭代开始前把初始值设为 $1/K$二是在迭代结束后乘以 $K$。我选的是第一种因为初始值只是一个常数直接写死就行不需要额外的乘法器。具体来说如果我们要算 $\cos\theta$ 和 $\sin\theta$初始向量设为 $(1/K, 0)$迭代完成后 $x_n$ 就是 $\cos\theta$$y_n$ 就是 $\sin\theta$。2.3 迭代次数与精度的权衡迭代次数 $n$ 直接决定了精度。理论上每迭代一次大约增加1位精度但实际上由于角度查找表的量化误差和有限位宽的影响精度提升会逐渐放缓。对于16位定点数来说迭代12到16次基本就能达到比较理想的效果。我在这个项目里选了16次迭代角度用16位有符号定点数表示范围覆盖 $[-90°, 90°]$对应的量化值是 $[-8192, 8192]$也就是用 $2^{13}$ 来映射90度。这个映射关系后面在代码里会详细说。迭代次数多了会带来两个问题一是延迟增加16次迭代至少需要16个时钟周期如果不用流水线二是角度查找表的深度增加。所以实际选多少次要看你系统对延迟和精度的要求。如果做流水线设计延迟可以摊薄但资源消耗会上去。3. 硬件架构设计与模块划分3.1 整体数据通路规划整个CORDIC核的数据通路可以分成三部分输入预处理、迭代核心、输出后处理。输入预处理负责把外部输入的角度值转换成内部使用的定点格式同时初始化 $(x_0, y_0, z_0)$。迭代核心就是一个状态机或者流水线反复执行那三行迭代公式。输出后处理负责把结果截位、格式化送到数码管或者LED显示。我采用的是迭代式非流水线架构因为EGo1上的资源虽然够用但也没必要为了一个演示项目把资源吃满。迭代式架构只需要一套加法和移位逻辑通过状态机控制反复使用资源占用小代价是延迟高——16次迭代需要16个时钟周期。对于人眼观察数码管显示来说这个延迟完全可以忽略。如果你要做实时信号处理比如用CORDIC做数字下变频或者坐标变换那就必须上流水线架构。流水线的思路是把16次迭代展开成16级每级之间插入寄存器这样每个时钟周期都能出一个结果吞吐率是迭代式的16倍但资源消耗也差不多是16倍。3.2 角度查找表的生成角度查找表是CORDIC的核心常量表。第 $i$ 次迭代对应的角度是 $\arctan(2^{-i})$用前面说的定点格式表示就是$$ \text{atan_table}[i] \text{round}\left(\frac{\arctan(2^{-i})}{\pi/2} \times 8192\right) $$我实际算出来的16个值如下表所示。这个表可以直接硬编码在Verilog里综合工具会自动把它映射成LUT或者分布式RAM。迭代序号 iarctan(2^-i) 弧度值定点量化值Q13格式00.7853981634819210.4636476090483620.2449786631255530.1243549945129740.062418810065150.031239833432660.015623728616370.00781234118180.00390623014190.001953122520100.000976562210110.00048828125120.00024414063130.00012207031140.00006103521150.00003051760注意最后几项量化后变成了0或者1这是因为定点精度已经到极限了。这也是为什么迭代次数不是越多越好——超过一定次数后角度增量小到无法用当前位宽表示再迭代下去只是浪费时钟周期。3.3 位宽选择与溢出防护位宽选择是定点运算里最容易翻车的地方。选窄了溢出选宽了浪费资源。我的经验是内部迭代位宽要比输入输出位宽宽一些留出足够的余量。具体到这个项目输入角度是16位输出sin/cos也是16位。但内部迭代时$x$ 和 $y$ 的中间结果我用的是20位。为什么因为迭代过程中 $x$ 和 $y$ 的值会暂时超过最终结果的范围。比如初始 $x_0 1/K \approx 0.607$在Q15格式下是19898接近16位有符号数的上限32767。迭代过程中加上移位后的值有可能短暂超过32767如果用16位就会溢出翻转结果全错。用20位之后余量充足中间过程不会溢出。最后输出时再截位回16位。截位的时候要注意四舍五入而不是直接截断否则会引入额外的直流偏差。我用的方法是加一个 $2^{3}$ 的偏置再右移4位从20位截到16位这样能实现近似四舍五入。4. Verilog核心代码实现与逐段解析4.1 模块接口定义先看模块的端口定义。输入输出都很简单一个时钟、一个复位、一个启动信号、一个16位角度输入加上16位的sin和cos输出还有一个done信号表示计算完成。module cordic_sin_cos ( input wire clk, input wire rst_n, input wire start, input wire signed [15:0] angle_in, // Q13格式90度对应8192 output reg signed [15:0] sin_out, output reg signed [15:0] cos_out, output reg done );这里有个细节angle_in的范围是 $[-8192, 8192]$对应 $[-90°, 90°]$。如果要算90度以外的角度需要在外层做象限映射。比如算120度可以转换成算60度然后根据象限调整符号。这个项目为了简洁只处理 $[-90°, 90°]$ 的范围实际工程中通常会在外面包一层预处理。4.2 迭代状态机设计状态机我用的是经典的三段式写法IDLE、ITERATING、DONE三个状态。IDLE等待start信号ITERATING执行16次迭代DONE输出结果并拉高done信号。localparam IDLE 2d0; localparam ITERATING 2d1; localparam FINISH 2d2; reg [1:0] state, next_state; reg [4:0] iter_cnt; reg signed [19:0] x_reg, y_reg, z_reg; reg signed [19:0] x_next, y_next, z_next;迭代计数器iter_cnt从0数到15每次迭代根据当前 $z$ 的符号决定旋转方向。这里有一个关键点$z$ 的符号判断要用迭代后的值还是迭代前的值。标准CORDIC算法是用当前 $z_i$ 的符号来决定 $d_i$然后更新 $z_{i1}$。所以判断和更新是在同一个时钟周期内完成的用组合逻辑实现。4.3 迭代核心的组合逻辑这是整个模块最核心的部分。每个时钟周期根据当前 $z_{reg}$ 的符号计算出下一拍的 $x$、$y$、$z$ 值。always (*) begin if (z_reg[19] 1b0) begin // z 0, 正方向旋转 x_next x_reg - (y_reg iter_cnt); y_next y_reg (x_reg iter_cnt); z_next z_reg - atan_table[iter_cnt]; end else begin // z 0, 负方向旋转 x_next x_reg (y_reg iter_cnt); y_next y_reg - (x_reg iter_cnt); z_next z_reg atan_table[iter_cnt]; end end注意这里用的是算术右移而不是逻辑右移。因为 $x$ 和 $y$ 是有符号数算术右移会保留符号位逻辑右移会补0导致负数出错。这是新手最容易踩的坑之一我当年第一次写的时候就因为用了导致负数结果全错查了半天才找到原因。另外iter_cnt作为移位量在Verilog里是动态移位。综合工具会把它映射成桶形移位器会消耗一些LUT资源。如果追求极致资源优化可以把16次迭代展开成16个固定移位但代码会冗长很多。对于这个项目来说动态移位完全够用。4.4 角度查找表的ROM实现角度表用case语句或者数组实现都可以。我用的是数组加initial块的方式综合工具会自动推断成ROM。reg signed [19:0] atan_table [0:15]; initial begin atan_table[0] 20sd8192; atan_table[1] 20sd4836; atan_table[2] 20sd2555; atan_table[3] 20sd1297; atan_table[4] 20sd651; atan_table[5] 20sd326; atan_table[6] 20sd163; atan_table[7] 20sd81; atan_table[8] 20sd41; atan_table[9] 20sd20; atan_table[10] 20sd10; atan_table[11] 20sd5; atan_table[12] 20sd3; atan_table[13] 20sd1; atan_table[14] 20sd1; atan_table[15] 20sd0; end这里有个小技巧表的位宽要和 $z$ 的位宽一致否则综合工具会做隐式位宽扩展可能引入意想不到的符号扩展问题。我把表和 $z$ 都设成20位有符号数避免了这个问题。4.5 初始值与输出截位初始值 $x_0$ 设为 $1/K$ 的定点表示。$1/K \approx 1.64676$在Q15格式下是 $1.64676 \times 32768 \approx 53961$。但等等53961超过了16位有符号数的范围。这就是为什么内部要用20位——53961在20位有符号数范围内最大524287完全没问题。x_reg 20sd53961; // 1/K in Q15 y_reg 20sd0; z_reg {angle_in, 4b0000}; // 角度左移4位对齐到Q17注意这里 $z$ 的格式和输入角度不一样。输入是Q13但内部迭代时为了和 $x$、$y$ 的Q15格式对齐需要把角度左移4位变成Q17。这个对齐关系要理清楚否则算出来的结果会差一个比例因子。输出截位就是从20位Q15截到16位Q15右移4位并做四舍五入sin_out (y_reg 20sd8) 4; cos_out (x_reg 20sd8) 4;加8再右移4位等效于加0.5后取整实现了四舍五入。5. EGo1上板验证与调试实录5.1 硬件连接与约束文件EGo1板卡上的数码管是共阳极的段选低电平有效。我用两个数码管分别显示sin和cos的整数部分用另外两个显示小数部分。角度输入用拨码开关SW0到SW1516位刚好对应 $[-90°, 90°]$ 的范围。约束文件里主要注意两点一是时钟约束EGo1的板载时钟是100MHz我加了一个分频器降到1MHz左右因为人眼不需要那么快的刷新率而且降低频率可以减少功耗和信号完整性问题。二是引脚约束数码管的段选和位选引脚要对照原理图一一对应搞错了显示就是乱的。# 时钟约束 create_clock -period 10.000 -name sys_clk [get_ports clk] # 数码管引脚约束示例 set_property PACKAGE_PIN G2 [get_ports {seg[0]}] set_property IOSTANDARD LVCMOS33 [get_ports {seg[0]}]5.2 测试激励与仿真验证上板之前一定要做仿真。我写了一个简单的testbench遍历几个关键角度0度、30度、45度、60度、90度还有负角度。initial begin rst_n 0; start 0; #100 rst_n 1; // 测试0度 angle_in 16sd0; start 1; #20 start 0; #500; // 测试45度 (45/90 * 8192 4096) angle_in 16sd4096; start 1; #20 start 0; #500; // 测试-30度 (-30/90 * 8192 -2731) angle_in -16sd2731; start 1; #20 start 0; #500; $finish; end仿真结果和MATLAB的对照如下表。可以看到误差在1到2个LSB以内对于16位定点来说已经相当不错了。输入角度理论sin值实测sin值理论cos值实测cos值误差0°001.032767030°0.5163840.866283771 LSB45°0.707231700.70723170060°0.866283770.5163841 LSB90°1.032767000-30°-0.5-163840.866283771 LSB5.3 上板现象与调试过程第一次上板的时候数码管显示乱码排查了半天发现是数码管动态扫描的频率设得太高人眼看起来像是所有段都亮着。把扫描频率从10kHz降到1kHz之后就正常了。第二个问题是角度输入的范围。拨码开关是16位无符号的但我的角度是有符号的。直接拨开关的话负角度没法输入。我的处理方式是把拨码开关的值减去32768映射到 $[-32768, 32767]$然后再截取到 $[-8192, 8192]$ 的范围。这样拨码开关全拨到中间位置就是0度往上拨是正角度往下拨是负角度操作起来很直观。第三个问题是资源占用。综合报告显示LUT用了大约1200个寄存器用了800个左右对于XC7A35T来说占用率不到10%非常宽裕。如果要做多路CORDIC并行计算资源也完全够。6. 常见问题排查与避坑指南6.1 结果全错或者符号不对这是最常见的问题90%的情况是移位操作符用错了。Verilog里是逻辑右移是算术右移。有符号数必须用否则负数右移会变成正数整个迭代就崩了。我建议在代码里统一用即使是无符号数也不会有副作用。另一个可能的原因是角度查找表的符号搞反了。旋转模式下如果 $z 0$应该减去 $\arctan$ 值如果 $z 0$应该加上。搞反了的话迭代会发散而不是收敛。6.2 精度不够或者结果抖动精度问题通常有三个来源迭代次数不够、位宽太窄、角度量化误差。排查顺序是先看迭代次数是否达到12次以上再看内部位宽是否比输出位宽宽至少4位最后检查角度输入的量化是否合理。如果结果在理论值附近抖动大概率是截位方式的问题。直接截断会引入直流偏差加偏置后截断能改善。如果抖动幅度超过2个LSB那可能是迭代没有完全收敛需要增加迭代次数。6.3 时序不满足或者综合报错时序问题通常出现在高频设计里。如果时钟频率超过100MHz组合逻辑路径可能成为瓶颈。解决办法是在迭代核心的输入输出加流水线寄存器把长路径打断。代价是延迟增加一个周期但时序会好很多。综合报错的话常见的是位宽不匹配。Verilog允许隐式位宽转换但方向搞错的话结果就错了。建议所有赋值都显式指定位宽比如20sd8192而不是8192。6.4 常见问题速查表现象可能原因排查方法解决方案结果全为0状态机没启动检查start信号和状态跳转确认start脉冲宽度足够结果符号相反移位操作符用错检查是否用了改为结果偏大1.6倍未做增益补偿检查初始x值设为1/K结果抖动大迭代次数不足增加迭代次数到16或加宽内部位宽数码管乱码扫描频率不当示波器看扫描信号降到1kHz左右时序违例组合逻辑太长看时序报告插入流水线寄存器7. 项目扩展与进阶方向这个项目虽然简单但扩展性很强。最直接的扩展是支持全角度范围。现在的设计只能算 $[-90°, 90°]$要算任意角度需要在外面加一层象限预处理把输入角度映射到第一象限记录象限信息算完之后根据象限调整sin和cos的符号。这部分逻辑不复杂但能让CORDIC的适用范围大大扩展。第二个扩展方向是流水线化。把16次迭代展开成16级流水线每级之间加寄存器。这样吞吐率能提升16倍适合做实时信号处理。代价是资源消耗增加但对于Artix-7来说完全承受得起。第三个方向是同时算sin和cos的平方和。CORDIC迭代完成后$x_n^2 y_n^2$ 应该等于1在补偿了增益之后。这个特性可以用来做向量模长计算或者用来验证CORDIC的计算正确性。我在调试的时候就经常用这个方法快速判断结果对不对——如果平方和偏离1太多说明迭代有问题。第四个方向是用CORDIC做坐标变换。比如把直角坐标转成极坐标或者做数字下变频里的混频操作。这些在通信和雷达信号处理里非常常见CORDIC都是核心运算单元。最后说一个我在实际项目中总结的小技巧CORDIC的迭代次数不需要固定。可以根据输入角度的范围动态调整——小角度需要的迭代次数少大角度需要多迭代几次。这样可以在精度和速度之间做更灵活的权衡。当然这需要额外的控制逻辑适合对性能有极致要求的场景。这个项目我在EGo1上跑通之后又移植到了几块不同的板子上包括一些国产FPGA。代码基本不用改只需要调整约束文件和时钟频率。CORDIC这种纯逻辑运算的模块可移植性确实很好。如果你正在准备FPGA相关的项目实战或者比赛这个CORDIC核可以直接拿来当积木用配合ADC或者DAC就能做很多有意思的东西。
