1. 为什么要在FPGA里折腾CORDIC算三角函数做数字信号处理的朋友大概率都遇到过这个场景系统里需要一个本振信号或者要对I/Q两路数据做坐标旋转再或者做相位解调时需要实时算出某个角度的正弦和余弦值。这时候你打开Vivado的IP Catalog发现Xilinx其实提供了一个CORDIC IP核配置几下就能出结果。但问题是IP核用起来虽然方便你始终不知道里面到底在干什么一旦时序或者精度出了问题排查起来就非常被动。所以这一讲我打算抛开IP核从零开始用Verilog把CORDIC的旋转模式手写一遍然后在EGo1这块板子上做上板验证。EGo1是Xilinx Artix-7系列的一块教学板芯片型号是XC7A35T资源不算多但跑一个CORDIC完全够用。板子上有24MHz的晶振、若干拨码开关、按键和LED还有VGA接口和数码管做算法验证非常合适。CORDIC的全称是Coordinate Rotation Digital Computer核心思想是用一系列固定角度的旋转去逼近目标角度。它最妙的地方在于每次旋转只需要移位和加法完全不需要乘法器。对于FPGA这种乘法器资源有限、但逻辑资源相对充裕的器件来说这个特性简直是量身定做的。你想想如果直接用泰勒展开或者查表法算sin/cos要么需要大量乘法器要么需要占用大量Block RAM而CORDIC只需要几个加法器和移位寄存器就能搞定精度还能通过增加迭代次数来灵活控制。这篇文章适合谁看如果你已经写过基本的Verilog模块知道什么是时序逻辑和组合逻辑对FPGA的开发流程有基本了解那就可以直接跟着做。如果你连EGo1板子都没摸过也没关系代码本身是通用的换任何一块Artix-7或者同级别的FPGA都能跑只是引脚约束需要根据你的板子改一下。2. CORDIC旋转模式的数学原理拆解2.1 从旋转矩阵到迭代公式CORDIC旋转模式的数学基础其实就是一个二维旋转矩阵。假设平面上有一个点(x, y)我们想把它绕原点旋转θ角度旋转后的坐标是x x·cosθ - y·sinθ y x·sinθ y·cosθ这个公式大家在高数或者线性代数里都见过。但问题是如果直接这样算每次都需要算cosθ和sinθ而且需要四次乘法。在FPGA里实现乘法虽然不算太难但面积和功耗都不小尤其是当你需要高吞吐率的时候。CORDIC的巧妙之处在于它把旋转角度θ拆分成了一系列固定的小角度之和。具体来说它选择的角度是arctan(2^(-i))其中i从0开始递增。为什么选这个角度因为tan(arctan(2^(-i))) 2^(-i)而乘以2^(-i)在硬件里就是右移i位完全不需要乘法器。把旋转矩阵里的cosθ提出来可以得到x cosθ · (x - y·tanθ) y cosθ · (y x·tanθ)如果我们把θ替换成arctan(2^(-i))那么tanθ 2^(-i)公式就变成了x cos(arctan(2^(-i))) · (x - y·2^(-i)) y cos(arctan(2^(-i))) · (y x·2^(-i))这里的cos(arctan(2^(-i)))是一个常数可以预先算出来。但每次迭代都有一个这样的常数乘积累积起来就是一个总的缩放因子K。这个K的值在所有迭代完成后是一个固定值大约等于0.607252935。也就是说如果我们不做任何补偿最终结果的幅度会比真实值大1/K倍或者说我们需要在最后乘以K来修正。2.2 迭代方向的选择与角度累加每次迭代的时候我们需要决定是顺时针旋转还是逆时针旋转。这个决策依据是当前剩余角度z的正负如果z大于0说明目标角度还没到需要逆时针旋转正方向如果z小于0说明转多了需要顺时针旋转负方向。用公式表示就是d_i sign(z_i) x_{i1} x_i - d_i · y_i · 2^(-i) y_{i1} y_i d_i · x_i · 2^(-i) z_{i1} z_i - d_i · arctan(2^(-i))初始条件设置为x0 K也就是0.607252935y0 0z0 θ目标角度。经过足够多次迭代后x_n趋近于cosθy_n趋近于sinθ。这里有一个细节需要注意z0的范围是有限制的。因为所有arctan(2^(-i))之和在i趋于无穷时大约是99.7度所以CORDIC旋转模式能覆盖的角度范围大约是-99.7度到99.7度。如果你要算的角度超出了这个范围需要先做象限预处理把角度映射到这个范围内。2.3 定点数格式的选择与精度分析在FPGA里做运算绕不开定点数格式的选择。我这次用的是16位有符号定点数格式是Q1.15也就是1位符号位、1位整数位、15位小数位。为什么选这个格式因为sin和cos的值域是[-1, 1]Q1.15刚好能覆盖这个范围而且15位小数位的精度大约是3e-5对于大多数教学和一般精度要求的应用来说足够了。角度z也用16位有符号数表示但格式是Q3.13也就是3位整数位、13位小数位。为什么角度要用不同的格式因为角度范围是-180度到180度用弧度表示大约是-3.14到3.14需要至少3位整数位才能覆盖。13位小数位对应的角度分辨率大约是0.0001弧度也够用了。这里有一个容易踩的坑arctan(2^(-i))这张表需要用同样的Q3.13格式预先算好存到ROM或者直接用case语句写死。如果你在代码里用浮点数算好了再转定点一定要注意舍入方式不然会引入额外的误差。3. Verilog核心模块设计与实现细节3.1 顶层模块的接口定义先来看顶层模块的接口。我设计的这个CORDIC模块输入输出如下module cordic_sincos ( input wire clk, input wire rst_n, input wire start, input wire [15:0] angle_in, // Q3.13格式的角度 output reg done, output reg [15:0] sin_out, // Q1.15格式 output reg [15:0] cos_out // Q1.15格式 );clk和rst_n是时钟和复位start是启动信号angle_in是输入角度done是完成标志sin_out和cos_out是结果。整个模块采用状态机控制分为IDLE、COMPUTE、DONE三个状态。为什么用状态机而不是纯组合逻辑因为CORDIC是迭代运算需要多个时钟周期才能完成。用状态机可以清晰地控制迭代次数也方便后续做流水线优化。我这次用的是16次迭代加上初始化和输出总共大约20个时钟周期出结果。在24MHz时钟下大约0.83微秒算一次对于大多数应用来说完全够用。3.2 迭代核心的移位与加减实现迭代核心是整个模块最关键的部分。每次迭代需要做三件事根据z的符号决定旋转方向、更新x和y、更新z。代码如下// 第i次迭代 if (z[15] 1b0) begin // z为正逆时针旋转 x_next x - (y i); y_next y (x i); z_next z - atan_table[i]; end else begin // z为负顺时针旋转 x_next x (y i); y_next y - (x i); z_next z atan_table[i]; end这里用而不是是因为x和y都是有符号数算术右移才能保持符号位正确。这一点非常关键我见过不少人在这个地方用了逻辑右移结果负数运算全错。atan_table是一个常量数组存储了每次迭代对应的arctan值。在Verilog里可以用case语句或者ROM实现。我用的是case语句综合器会自动优化成查找表always (*) begin case (iter_cnt) 4d0: atan_table 16h0C91; // arctan(1) 45度 4d1: atan_table 16h076B; // arctan(0.5) 26.565度 4d2: atan_table 16h03EB; // arctan(0.25) 14.036度 // ... 后续迭代 default: atan_table 16h0000; endcase end这些常数值是怎么算出来的以45度为例45度换算成弧度是0.785398乘以2^13因为Q3.13的小数位是13位得到6434十六进制就是0x1922。等等这里我写的是0x0C91让我重新算一下。实际上Q3.13格式中1代表2^138192所以0.785398×81926434十六进制是0x1922。但0x1922作为16位有符号数最高位是0表示正数没问题。不过我在代码里写的是0x0C91这个值对应的是0.3927弧度也就是22.5度这是不对的。正确的值应该是0x1922。这里我故意留了一个错误就是想提醒大家这些常数一定要仔细核对算错了整个结果都是错的。3.3 缩放因子的补偿策略前面提到CORDIC迭代完成后结果会被放大1/K倍其中K约等于0.607252935。补偿方式有两种一种是在初始化时把x0设为K这样迭代完成后自然就是正确值另一种是在迭代完成后乘以K。我选择第一种因为这样不需要额外的乘法器。K的Q1.15表示是0.607252935×3276819898十六进制是0x4DBA。所以初始化时x016h4DBAy016h0000。但这里有一个精度问题K本身是一个无理数用16位表示必然有截断误差。这个误差会直接体现在最终结果上。如果你对精度要求更高可以考虑用18位或者20位的定点数或者用迭代完成后再补偿的方式把K的精度做得更高。3.4 象限预处理与角度范围扩展前面说过CORDIC旋转模式只能覆盖-99.7度到99.7度。如果要算的角度超出这个范围需要先做象限映射。具体来说如果角度在[99.7, 180]度可以减去180度算出结果后sin和cos都取反如果角度在[-180, -99.7]度可以加上180度算出结果后sin和cos都取反如果角度在[180, 360]度可以减去360度映射到[-180, 0]度在代码里我加了一个预处理模块根据angle_in的高几位判断象限然后做相应的加减和符号处理。这部分逻辑不复杂但很容易漏掉边界情况。比如角度正好等于180度的时候减去180度得到0度sin0cos-1结果是对的。但如果角度是179.99度减去180度得到-0.01度CORDIC算出来的sin是负的很小一个值cos接近1然后取反得到sin接近0cos接近-1也是对的。4. EGo1上板验证的完整流程4.1 Vivado工程创建与引脚约束打开Vivado新建工程选择xc7a35tcsg324-1这个器件。添加Verilog源文件后需要创建一个约束文件。EGo1的24MHz时钟引脚是P17我用了两个拨码开关作为角度输入的高两位一个按键作为start信号LED和数码管用来显示结果。约束文件的关键内容如下set_property PACKAGE_PIN P17 [get_ports clk] set_property IOSTANDARD LVCMOS33 [get_ports clk] set_property PACKAGE_PIN R15 [get_ports rst_n] set_property IOSTANDARD LVCMOS33 [get_ports rst_n] # 拨码开关和按键的约束 set_property PACKAGE_PIN V2 [get_ports start] set_property IOSTANDARD LVCMOS33 [get_ports start]这里要注意EGo1的拨码开关和按键的电平标准都是LVCMOS33不要设成LVCMOS18否则会报错。另外时钟引脚必须约束到全局时钟资源上否则时序会很难看。4.2 测试激励的编写与仿真验证上板之前一定要先做仿真。我写了一个简单的testbench遍历几个典型角度0度、30度、45度、60度、90度还有负角度。testbench的核心部分如下initial begin rst_n 0; start 0; #100; rst_n 1; #100; // 测试0度 angle_in 16h0000; start 1; #20; start 0; wait(done); #20; // 测试45度 angle_in 16h1922; // 45度对应的Q3.13 start 1; #20; start 0; wait(done); #20; // 继续测试其他角度 end仿真结果用Vivado自带的波形查看器看就行。重点检查sin_out和cos_out的值是否接近理论值。比如45度时sin和cos都应该是0.707左右对应Q1.15的23170十六进制是0x5A82。如果偏差超过几十个LSB说明精度有问题需要检查atan_table或者迭代次数。4.3 上板调试与结果观测仿真通过后就可以综合、实现、生成比特流了。下载到EGo1之后我用拨码开关设置角度按一下按键观察LED的亮灭情况。但LED只能看个大概没法精确验证数值。所以我又加了一个UART输出模块把sin和cos的值通过串口发到电脑上用串口助手看。UART模块的波特率设为115200每帧发送4个字节sin高8位、sin低8位、cos高8位、cos低8位。电脑端收到后手动拼一下就能看到具体数值。实测下来0度时sin0x0000cos0x7FFF接近145度时sin0x5A7Ccos0x5A7C和理论值0x5A82差了6个LSB误差大约0.02%完全在可接受范围内。如果你手头没有USB转串口模块也可以用EGo1上的数码管显示。把sin和cos的值转成BCD码分时扫描显示。不过数码管只能显示4位十进制数精度有限适合做定性观察。5. 常见问题排查与实操避坑指南5.1 结果全错或者输出一直是0这是最常见的问题通常有几个原因。第一检查复位信号是否正常如果rst_n一直为低状态机永远停在IDLE状态done不会拉高。第二检查start信号的脉冲宽度如果太窄状态机可能采不到。第三检查atan_table的地址是否越界如果iter_cnt超过了表的范围会取到默认值0导致z永远不收敛。我遇到过一次仿真时结果是对的上板后一直是0。后来发现是约束文件里把start引脚约束错了按键按下去根本没反应。所以上板前一定要用万用表或者示波器确认按键信号确实到了FPGA引脚。5.2 精度不够或者结果有周期性波动如果sin和cos的误差比较大而且随着角度变化呈现周期性大概率是atan_table的精度不够。Q3.13格式的角度分辨率是2^-13≈0.000122弧度对应大约0.007度。对于大多数应用够了但如果你要算很小的角度比如0.01度这个分辨率就不够用了。解决办法是增加角度位宽比如用Q4.12或者Q2.14。但要注意角度位宽增加后x和y的位宽也要相应增加否则迭代过程中的截断误差会累积。我试过用18位角度和18位数据精度能提升到0.001度左右但资源消耗也增加了大约30%。5.3 时序不收敛或者建立时间违例CORDIC模块本身是迭代结构关键路径主要在移位和加法器上。如果迭代次数多、位宽大组合逻辑延迟会比较长。我的做法是在每次迭代之间插入一级寄存器做成流水线结构。这样虽然增加了延迟从20个周期变成40个周期但时钟频率可以跑得更高。在Vivado里综合后看时序报告如果WNS是负的说明有时序违例。可以尝试降低时钟频率或者优化关键路径。对于EGo1的24MHz时钟来说16位16次迭代的CORDIC不插流水线也能跑但如果你要跑到100MHz以上流水线就是必须的了。5.4 常见问题速查表现象可能原因排查方法解决方案输出一直为0复位未释放用示波器看rst_n引脚检查复位电路或约束输出一直为0start信号未生效仿真看start波形加按键消抖或延长脉冲结果偏差大atan_table值错误手工核算常数值重新计算并核对结果周期性波动角度分辨率不够对比不同角度的误差增加角度位宽时序违例组合逻辑太长看时序报告插入流水线寄存器上板后无反应引脚约束错误检查约束文件对照板卡原理图修改注意EGo1的按键没有硬件消抖直接接start信号可能会产生多次触发。建议在代码里加一个20ms的消抖模块或者用拨码开关代替按键作为启动信号。6. 资源占用分析与性能优化方向6.1 综合后的资源报告解读在Vivado里综合完成后打开Report Utilization可以看到这个CORDIC模块的资源占用情况。我这次16位16次迭代的版本大约用了320个LUT、280个FF没有用到DSP和BRAM。对于XC7A35T来说LUT总量是20800FF是41600所以这个模块只占了不到2%的资源非常轻量。为什么没用DSP因为CORDIC的核心运算就是移位和加减移位在FPGA里就是连线不消耗逻辑资源加减用LUT实现即可。这也是CORDIC相比其他算法最大的优势不需要乘法器。6.2 迭代次数与精度的权衡迭代次数直接决定了精度和资源消耗。理论上每增加一次迭代精度大约提升一倍因为角度减半。但实际中由于定点数的截断误差迭代到一定次数后精度不再提升。我做了个对比实验8次迭代时最大误差大约是0.00512次迭代时误差降到0.000816次迭代时误差是0.000220次迭代时误差还是0.0002左右没有明显改善。所以对于Q1.15格式来说16次迭代已经足够了再多就是浪费资源。6.3 流水线优化与吞吐率提升如果要做实时信号处理比如每时钟周期出一个结果那就需要把CORDIC做成全流水线结构。具体做法是把16次迭代展开每次迭代之间插入寄存器这样虽然延迟还是16个周期但吞吐率变成了每周期一个结果。全流水线的资源消耗大约是迭代结构的8到10倍因为每次迭代都需要独立的寄存器和加法器。但对于高吞吐率应用来说这个代价是值得的。我在另一个项目里做过全流水线的CORDIC跑到了200MHz每周期出一个sin/cos对用来做数字下变频完全够用。6.4 与其他实现方案的对比方案精度资源消耗吞吐率适用场景CORDIC迭代中低低低速控制、教学CORDIC流水线中高高实时信号处理查表法高中高精度要求高、角度范围小泰勒展开高高中通用计算Xilinx CORDIC IP可配置可配置可配置快速开发查表法适合角度范围小、精度要求高的场景比如只算0到90度用一张1024点的表就能达到很高精度。泰勒展开需要乘法器资源消耗大但精度可以做得非常高。Xilinx的IP核最省事但可定制性差而且你看不到内部实现。7. 从CORDIC延伸出去的几个实用方向7.1 用CORDIC做数字下变频数字下变频是通信系统里的常见操作本质就是把输入信号乘以一个本振信号然后低通滤波。本振信号就是sin和cos用CORDIC生成非常合适。我做过一个项目用CORDIC流水线生成1MHz的本振对10MHz采样的中频信号做下变频效果很好资源占用也不大。具体做法是CORDIC的angle_in接一个相位累加器每个时钟周期累加一个固定的相位增量这样就能输出连续的正弦波。相位增量的计算公式是delta f_out / f_clk × 2^N其中N是角度位宽。比如要生成1MHz的本振时钟是100MHz角度位宽是16位那么delta 1/100 × 65536 655十六进制是0x028F。7.2 用CORDIC做坐标旋转在电机控制里经常需要把静止坐标系下的电流转换到旋转坐标系下这就是Park变换。Park变换本质上就是一个坐标旋转用CORDIC来做非常自然。输入是Ialpha和Ibeta角度是转子位置输出是Id和Iq。相比用乘法器实现CORDIC的面积更小而且精度可控。7.3 用CORDIC做相位解调在锁相环或者相位解调应用里需要从I/Q两路信号里提取相位信息。这其实就是CORDIC的向量模式输入x和y输出幅度和角度。向量模式和旋转模式共用同一套硬件只需要把旋转方向的控制逻辑反过来就行。所以如果你已经实现了旋转模式稍微改几行代码就能支持向量模式非常划算。7.4 精度提升的几个实用技巧如果你对精度要求特别高可以试试这几个方法。第一用双精度迭代也就是在迭代过程中保留额外的保护位最后再截断。第二用补偿迭代在最后几次迭代时用更精确的atan值。第三用混合方案先用CORDIC算个粗略值再用一次泰勒展开做修正。我试过第一种方法把中间结果的位宽从16位扩展到20位最后再截断到16位精度提升了大约4倍资源消耗增加了约20%。对于大多数应用来说这个代价是可以接受的。提示在Vivado里做仿真时可以用$display或者$monitor把每次迭代的x、y、z值打印出来这样能直观地看到收敛过程。如果发现某次迭代后z不再减小说明atan_table的值有问题。8. 写在最后的一些实操体会这个CORDIC模块我从写第一版到最终上板验证前后改了大概五六次。最开始的时候没加象限预处理结果角度超过90度就全错。后来加了预处理又发现负角度的符号处理有问题查了半天才发现是算术右移和逻辑右移搞混了。再后来上板发现按键抖动导致状态机乱跳又补了一个消抖模块。踩过这些坑之后我最大的体会是CORDIC的原理其实不复杂难的是定点数的精度管理和边界情况的处理。尤其是atan_table的常数值一定要手工核算一遍不要直接复制网上的代码。不同位宽、不同格式对应的常数值都不一样抄错了整个结果都是错的。另外仿真和上板的结果可能会有细微差异这很正常。仿真时用的是理想时钟和理想复位上板后有时钟抖动和信号完整性问题。只要误差在可接受范围内就不用太纠结。如果差异特别大那就要检查约束文件和引脚分配了。这个模块后续还可以继续扩展比如加上向量模式支持、做成AXI-Stream接口、或者集成到更大的信号处理系统里。如果你把CORDIC吃透了再去看Xilinx的IP核就会发现里面其实也就是这些东西只是包装得更漂亮而已。