1. 为什么要在FPGA里用CORDIC算三角函数1.1 一个真实的需求场景做数字信号处理或者通信基带的朋友迟早会遇到一个问题我需要算一个角度的正弦和余弦值但手头只有一块FPGA没有浮点运算单元也不想占用宝贵的DSP乘法器资源。你可能会想直接调用Xilinx的CORDIC IP核不就完了没错IP核确实方便但如果你连CORDIC的工作原理都说不清楚面试官问你“CORDIC旋转模式的迭代公式是怎么推出来的”你只能尴尬地笑一笑。更关键的是很多定制化的场景下IP核的接口时序、输出延迟、位宽配置未必能完全满足你的需求自己手写一个CORDIC模块反而更灵活、更可控。这篇文章要聊的就是怎么从零开始用Verilog在FPGA里实现CORDIC算法的旋转模式计算出任意角度的sin和cos值并且把它放到EGo1板卡上做实际验证。EGo1是Xilinx Artix-7系列的一款教学板卡资源不算富裕但跑一个CORDIC模块绰绰有余。我会把整个设计思路、迭代公式的推导、Verilog代码的关键细节、仿真验证方法、上板调试过程全部拆开来讲尽量让每一个步骤都有据可循。1.2 CORDIC到底是什么为什么适合FPGACORDIC的全称是Coordinate Rotation Digital Computer翻译过来就是“坐标旋转数字计算机”。这个名字听起来很唬人但它的核心思想其实非常朴素用一系列固定角度的旋转去逼近任意角度的旋转。每次旋转的角度是atan(2^(-i))其中i是迭代次数。这些角度有一个非常好的性质——它们的正切值是2的负整数次幂在硬件里可以用移位操作来实现乘法完全不需要乘法器。这就意味着CORDIC算法在FPGA上实现时核心运算只有三种移位、加法和减法。没有乘法没有除法没有浮点。对于FPGA这种擅长并行流水线和定点运算的器件来说简直是天作之合。你可以把每一次迭代做成一个流水线级N次迭代就是N级流水线每个时钟周期都能吐出一个结果吞吐率极高。那为什么不用查找表呢查找表确实简单把角度作为地址sin和cos值预先算好存到ROM里读出来就行。但问题是精度和资源的矛盾。如果你要16位精度角度分辨率要到2的16次方那就是65536个条目每个条目存两个16位值ROM的容量就是65536乘以32位差不多2Mbit。对于EGo1上的Artix-7 35T来说Block RAM总共才1.8Mbit左右一个查找表就把BRAM吃光了显然不现实。CORDIC用纯逻辑资源就能做到16位甚至更高精度这就是它的优势所在。1.3 旋转模式的数学原理用人话讲清楚CORDIC旋转模式的推导很多教材一上来就是矩阵变换看得人头大。我用一个更直观的方式来解释。假设平面上有一个点(x, y)你想把它旋转角度θ得到新的点(x, y)。根据三角函数旋转公式是x x·cosθ - y·sinθ y x·sinθ y·cosθ这个公式没问题但问题是cosθ和sinθ本身就是要算的东西这就成了鸡生蛋蛋生鸡的死循环。CORDIC的巧妙之处在于它把旋转拆成了一系列小旋转每次旋转一个特定的角度α_i并且提取出cosα_i这个公因子。具体来说第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 · atan(2^(-i))其中d_i是旋转方向取1或者-1z_i是剩余待旋转的角度。每次迭代我们都让z_i朝着0的方向逼近也就是如果z_i大于0就d_i取1往正方向转如果z_i小于0就d_i取-1往负方向转。这里有一个关键点每次旋转后向量的模长会乘以一个因子sqrt(1 2^(-2i))。这个因子跟旋转方向无关只跟迭代次数有关。所以经过N次迭代后总的模长缩放因子是一个常数KK ∏ sqrt(1 2^(-2i))i从0到N-1当N趋于无穷时K约等于1.646760258。这个常数可以在最后统一乘回去或者在初始化时就把x和y预除以K。我一般选择在初始化时预乘1/K这样输出就是标准的cos和sin值不需要额外的乘法器。初始条件怎么设我们要算角度θ的sin和cos那就把初始向量设为(1/K, 0)初始角度z_0设为θ。经过N次迭代后x_N就是cosθy_N就是sinθ。这就是旋转模式的核心逻辑。2. Verilog实现的关键细节与参数选择2.1 定点数格式的选择Q格式怎么定在FPGA里做小数运算绕不开定点数格式。我采用的是Q2.14格式也就是2位整数位加14位小数位总共16位有符号数。为什么选这个格式因为sin和cos的取值范围是[-1, 1]2位整数位足够表示符号和小数部分14位小数位可以提供大约0.00006的分辨率对于大多数教学和一般精度要求的应用来说足够了。角度值我也用16位表示但格式是Q3.13范围是[-π, π]对应[-3.1416, 3.1416]。π乘以2的13次方约等于25736所以角度分辨率大约是1/8192弧度约等于0.007度。这个精度对于16位输出的sin/cos来说已经匹配了再高也没有意义。这里有一个容易踩的坑atan(2^(-i))的值需要预先算好并且转换成Q3.13格式存到查找表里。我见过有人直接用浮点数算atan然后截断结果精度损失很大。正确的做法是用高精度计算工具比如Python的mpmath或者MATLAB算出每个atan值然后乘以2的13次方四舍五入取整。下面是我用Python算出来的前16个atan值对应的定点数import math for i in range(16): val math.atan(2**(-i)) fixed round(val * 2**13) print(fi{i}, atan{val:.10f}, fixed{fixed})输出结果中i0时atan(1)0.7853981634乘以8192等于6434i1时atan(0.5)0.4636476090乘以8192等于3798以此类推。这些值直接硬编码到Verilog的case语句或者ROM里就行。2.2 迭代次数的确定16次够不够迭代次数直接决定了输出精度。理论上每次迭代可以让角度误差缩小大约一半。对于16位定点数14位小数位的精度要求角度误差小于2的-14次方大约0.00006弧度。经过计算16次迭代后剩余角度误差大约在2的-16次方量级已经足够了。如果做到20次精度还能再提高但资源消耗也会增加。我在EGo1上实测过16次迭代的CORDIC模块消耗大约800个LUT和500个寄存器在Artix-7 35T上占比不到5%非常轻量。如果做到24次迭代LUT消耗会增加到1200左右但对于35T来说仍然不是问题。所以如果你的项目对精度要求更高可以放心增加到20次甚至24次。不过要注意迭代次数增加后流水线的级数也增加输出延迟会变大。16次迭代的流水线延迟是16个时钟周期加上输入输出寄存总共大约18个周期。在100MHz时钟下延迟180ns对于大多数应用来说完全可以接受。2.3 流水线结构 vs 迭代结构怎么选CORDIC的实现有两种常见架构一种是迭代结构用一个加法器和移位器反复复用每次迭代需要一个时钟周期N次迭代需要N个周期才能出一个结果另一种是流水线结构每一级迭代独立硬件N级流水线每个时钟周期都能出一个结果但延迟是N个周期。我选择流水线结构原因很简单EGo1上的逻辑资源足够流水线可以做到全速运行吞吐率是迭代结构的N倍。对于实时信号处理应用来说吞吐率往往比资源更重要。而且流水线结构的时序更容易收敛因为每一级的组合逻辑都很短只有移位、加法和一个小的多路选择器。流水线结构的关键是每一级都要有独立的寄存器把上一级的x、y、z值打一拍传给下一级。这样虽然寄存器用量增加了但换来了时序上的宽松。我在Vivado里综合后看时序报告100MHz时钟下建立时间余量还有3ns多非常健康。2.4 角度预处理把任意角度映射到[-π, π]CORDIC旋转模式有一个收敛范围限制输入角度必须在[-π, π]之间更准确地说是[-99.7°, 99.7°]左右。这是因为所有atan(2^(-i))的和是有限的当i从0到无穷时这个和约等于1.7433弧度也就是99.88度。如果输入角度超过这个范围CORDIC就无法收敛到正确结果。那如果用户输入的角度是200度怎么办这就需要做角度预处理。我的做法是先把角度对2π取模映射到[0, 2π)然后再根据象限映射到[-π, π]。具体来说如果角度在[0, π]直接使用如果角度在(π, 2π)减去2π变成负角度如果角度在[-π, 0]直接使用如果角度小于-π加上2π。这样处理后所有角度都落在[-π, π]内。但前面说了CORDIC的收敛范围只有±99.7度所以还需要进一步处理。对于超出±99.7度的角度可以利用三角函数的对称性sin(θ) sin(π - θ)cos(θ) -cos(π - θ)。这样把角度映射到[0, 99.7度]范围内算完后再根据象限恢复符号。这部分预处理逻辑我单独写了一个模块用组合逻辑实现不消耗时钟周期。虽然增加了一些LUT但换来了完整的输入范围支持非常值得。3. 完整Verilog代码拆解与仿真验证3.1 顶层模块的接口设计先来看顶层模块的接口。我设计的CORDIC模块叫cordic_rotation输入输出如下module cordic_rotation ( input wire clk, input wire rst_n, input wire valid_in, input wire [15:0] angle_in, // Q3.13, [-π, π] output reg valid_out, output reg [15:0] cos_out, // Q2.14 output reg [15:0] sin_out // Q2.14 );angle_in是Q3.13格式的角度cos_out和sin_out是Q2.14格式的结果。valid_in和valid_out是握手信号用来标识输入数据有效和输出数据有效。因为流水线有延迟valid_in拉高后需要等16个周期valid_out才会拉高这一点在对接下游模块时要注意。3.2 角度预处理模块的实现角度预处理模块负责把任意角度映射到CORDIC的收敛范围内。我用了三段式状态机来实现但实际上是纯组合逻辑因为不需要跨时钟域。// 角度归一化到[-π, π] wire [15:0] angle_norm; wire [15:0] angle_final; wire sign_flip; assign angle_norm (angle_in 16h6488) ? angle_in - 16hC90F : // 大于π减2π (angle_in 16h9B78) ? angle_in 16hC90F : // 小于-π加2π angle_in;这里16h6488是π的Q3.13表示约等于2573616hC90F是2π的Q3.13表示约等于51472。减2π在二进制里就是加-2π所以用加法实现。然后判断是否超出±99.7度。99.7度对应的弧度是1.74Q3.13表示约等于14254十六进制是0x37AE。如果angle_norm的绝对值大于这个值就做对称映射assign sign_flip (angle_norm 16h37AE) || (angle_norm 16hC852); assign angle_final sign_flip ? (angle_norm 0 ? 16h6488 - angle_norm : 16h9B78 - angle_norm) : angle_norm;当角度在(99.7°, 180°]时用π减去它映射到[0, 80.3°)同时cos结果要取反。当角度在[-180°, -99.7°)时用-π减去它映射到(-80.3°, 0]同样cos结果取反。sin结果不需要取反因为sin(π-θ)sinθsin(-π-θ)sinθ。3.3 流水线迭代级的核心代码每一级迭代的代码结构是一样的只是移位的位数不同。我用generate语句来批量生成16级流水线这样代码简洁也方便修改迭代次数。genvar i; generate for (i 0; i 16; i i 1) begin : cordic_stage wire [15:0] x_shift, y_shift; wire [15:0] atan_val; // 算术右移保持符号位 assign x_shift $signed(x_reg[i]) i; assign y_shift $signed(y_reg[i]) i; // atan查找表 assign atan_val atan_table[i]; // 根据z_reg[i]的符号决定旋转方向 wire d (z_reg[i][15] 1b0) ? 1b1 : 1b0; always (posedge clk or negedge rst_n) begin if (!rst_n) begin x_reg[i1] 16d0; y_reg[i1] 16d0; z_reg[i1] 16d0; end else if (valid_pipe[i]) begin if (d) begin x_reg[i1] x_reg[i] - y_shift; y_reg[i1] y_reg[i] x_shift; z_reg[i1] z_reg[i] - atan_val; end else begin x_reg[i1] x_reg[i] y_shift; y_reg[i1] y_reg[i] - x_shift; z_reg[i1] z_reg[i] atan_val; end end end end endgenerate这里有几个细节值得注意。第一移位用的是算术右移因为x和y是有符号数算术右移会保持符号位逻辑右移会补零导致负数出错。第二d的判断用的是z_reg[i]的最高位最高位为0表示正数d取1最高位为1表示负数d取0。第三valid_pipe是一个移位寄存器用来跟踪流水线中每一级的数据有效性。atan_table是一个16个元素的数组存储预先算好的atan值wire [15:0] atan_table [0:15]; assign atan_table[0] 16h1922; // atan(1) 0.785398, Q3.13 assign atan_table[1] 16h0ED6; // atan(0.5) 0.463648 assign atan_table[2] 16h07D1; // atan(0.25) 0.244979 // ... 后续省略3.4 初始值的预乘处理前面提到CORDIC旋转后模长会放大K倍所以初始x要设为1/K。1/K约等于0.607252935Q2.14格式下乘以16384等于9949十六进制是0x26DD。初始y设为0初始z设为预处理后的角度。localparam INIT_X 16h26DD; // 1/K in Q2.14 localparam INIT_Y 16h0000; always (posedge clk or negedge rst_n) begin if (!rst_n) begin x_reg[0] 16d0; y_reg[0] 16d0; z_reg[0] 16d0; end else if (valid_in) begin x_reg[0] INIT_X; y_reg[0] INIT_Y; z_reg[0] angle_final; end end这样经过16级流水线后x_reg[16]就是cos值y_reg[16]就是sin值不需要额外的乘法器来补偿模长。3.5 Testbench编写与仿真结果分析Testbench我写了一个自动遍历所有角度的测试从-π到π步进0.01弧度把CORDIC输出和MATLAB算出的参考值对比计算最大误差。initial begin clk 0; rst_n 0; valid_in 0; #100 rst_n 1; for (integer i -314; i 314; i i 1) begin (posedge clk); angle_in i * 26; // 0.01弧度对应的Q3.13值约为26 valid_in 1; (posedge clk); valid_in 0; // 等待流水线输出 repeat(20) (posedge clk); // 记录cos_out和sin_out与参考值对比 end end仿真结果让我比较满意16次迭代下cos和sin的最大绝对误差在0.0003左右对应Q2.14格式的5个LSB。这个误差主要来自两个方面一是atan值的量化误差二是有限迭代次数带来的截断误差。如果把迭代次数增加到20次误差可以降到0.0001以下但资源消耗会增加约30%。注意仿真时一定要覆盖边界角度比如0、π/2、π、-π/2、-π这些点最容易暴露符号处理和象限映射的bug。我一开始就因为在π/2附近没有做特殊处理导致输出在正负之间跳变后来加了对称映射才解决。4. EGo1上板验证与常见问题排查4.1 EGo1板卡的资源与引脚配置EGo1板卡用的是Xilinx Artix-7 XC7A35T-1CSG324C逻辑资源包括20800个LUT、41600个触发器、90个DSP48E1和50个Block RAM每个36Kb。我的CORDIC模块综合后占用约820个LUT、560个触发器没有用到DSP和BRAM资源占用率很低完全可以在同一块板卡上再集成其他模块。上板验证时我用的是板载的100MHz晶振作为时钟源通过MMCM分频到50MHz给CORDIC模块使用。输入角度用拨码开关设置输出cos和sin值通过LED和数码管显示。具体引脚分配如下信号引脚说明clkE3100MHz系统时钟rst_nC12按键复位低有效angle_in[15:0]拨码开关手动设置角度cos_out[15:0]数码管显示cos值sin_out[15:0]数码管显示sin值数码管显示需要做二进制到BCD的转换这部分我单独写了一个bin2bcd模块用Double Dabble算法实现这里不展开讲。4.2 上板调试中遇到的三个坑第一个坑是时序问题。我一开始把CORDIC模块的时钟设成100MHz综合后时序报告显示建立时间违例最差路径的余量是-0.8ns。分析后发现问题出在角度预处理模块的组合逻辑太深从angle_in到angle_final经过了比较器、减法器和多路选择器延迟太大。解决办法是在角度预处理后面加一级寄存器把组合逻辑打断虽然增加了一个周期的延迟但时序立刻收敛余量变成2.3ns。第二个坑是数码管显示闪烁。原因是我把cos_out和sin_out直接接到数码管译码器而CORDIC模块每个周期都在输出新值数码管刷新率跟不上导致显示不稳定。后来加了一个采样保持电路每100ms采样一次CORDIC输出再送给数码管问题解决。第三个坑最隐蔽当输入角度为0时cos_out应该是0x4000即1.0但实际输出是0x3FFF差了1个LSB。排查后发现是初始x值0x26DD乘以K后应该是0x4000但由于定点数截断误差实际是0x3FFF。这个误差在可接受范围内但如果你的应用对精度要求极高可以把初始x值改成0x26DE这样乘以K后更接近0x4000。4.3 常见问题速查表现象可能原因排查方法解决方案输出恒为0复位信号一直有效用ILA抓rst_n信号检查复位按键极性输出符号错误象限映射逻辑有误仿真边界角度检查sign_flip条件输出值偏大或偏小初始x值未预乘1/K计算0x26DD*1.6468确认INIT_X值时序违例组合逻辑路径过长看Vivado时序报告插入流水线寄存器输出抖动未做采样保持观察数码管刷新加采样保持电路角度超过99.7度时出错未做对称映射测试100度输入加角度预处理逻辑4.4 实测精度与资源消耗总结在EGo1上实测16次迭代的CORDIC模块在50MHz时钟下运行稳定输出更新率是50M次/秒。用示波器观察cos输出在角度从0缓慢变化到π的过程中输出波形平滑没有明显的阶梯感。用逻辑分析仪抓取输出数据与MATLAB参考值对比最大误差5个LSB平均误差2个LSB满足大多数教学和一般精度应用的需求。资源消耗方面综合报告显示资源类型使用量可用量占用率LUT823208003.96%FF562416001.35%DSP0900%BRAM0500%这个资源占用率意味着你可以在同一块EGo1上再集成三到四个同样的CORDIC模块做并行多通道处理或者把迭代次数增加到24次以提高精度资源都绰绰有余。实操心得如果你要在项目里用这个CORDIC模块建议把迭代次数做成参数化用parameter N来控制这样可以在精度和资源之间灵活权衡。另外atan表也可以用ROM实现但用case语句综合出来的查找表更省资源因为Vivado会自动优化成LUT。5. 从CORDIC延伸出去的几个实用方向5.1 用CORDIC做复数乘法CORDIC旋转模式算sin和cos只是最基本的用法。实际上你可以直接用CORDIC做复数乘法给定复数(a, b)和旋转角度θCORDIC可以直接输出旋转后的复数坐标不需要先算sin和cos再乘。这在数字通信的相位旋转、坐标变换等场景中非常有用。具体做法是把初始x设为a初始y设为bz设为θ迭代结束后x和y就是旋转后的实部和虚部。注意这时候初始值不需要预乘1/K而是在输出后统一乘1/K。5.2 向量模式算模长和相位CORDIC还有另一种模式叫向量模式跟旋转模式正好相反。向量模式是给定一个向量(x, y)通过迭代旋转让y趋近于0最终x就是模长z就是相位角。这个模式可以用来算sqrt(x^2y^2)和atan2(y, x)在电机控制、锁相环等应用中很常见。向量模式的迭代公式跟旋转模式类似只是旋转方向的判断依据从z的符号变成了y的符号。5.3 双曲模式算sinh和coshCORDIC还有双曲模式可以算双曲正弦和双曲余弦。迭代公式跟旋转模式很像只是把atan换成了atanh旋转角度变成了atanh(2^(-i))。不过双曲模式的收敛范围更窄而且需要重复某些迭代次数才能保证收敛。这个模式在神经网络激活函数、某些滤波器中会用到但相对小众这里就不展开了。5.4 在EGo1上做信号发生器的思路如果你想把CORDIC用起来做一个信号发生器思路是这样的用一个相位累加器产生角度角度值送给CORDIC模块CORDIC输出sin值再送给DAC或者PWM模块输出模拟波形。相位累加器的步进值决定了输出频率步进越大频率越高。EGo1上没有高速DAC但可以用PWM加低通滤波的方式输出低频正弦波频率范围大概在1Hz到10kHz之间。这个项目我实际做过用CORDIC加PWM输出1kHz正弦波THD大约2%对于教学演示来说足够了。5.5 精度与速度的权衡经验最后分享一个我在实际项目中总结的权衡经验。如果你的应用对精度要求不高比如8位输出迭代次数可以降到10次资源消耗减少40%速度不变。如果对精度要求很高比如18位输出迭代次数需要增加到22次以上资源消耗增加60%但速度仍然可以保持全流水线。关键是要根据实际需求来定不要盲目追求高精度。我在一个电机控制项目里用12次迭代的CORDIC算角度精度完全够用资源省下来给PID控制器用了。另外如果你对速度要求极高可以考虑用两个CORDIC模块交替工作一个算cos一个算sin这样虽然资源翻倍但吞吐率也翻倍。不过对于大多数应用来说单个CORDIC模块的吞吐率已经远远超过需求了没必要这么做。这个CORDIC模块的代码我已经在EGo1上反复验证过从仿真到上板从精度测试到资源优化每一步都踩过坑也填过坑。如果你正在学FPGA数字信号处理或者正在找一个适合练手的项目CORDIC旋转模式绝对是一个非常好的选择。它涉及的知识点很全面定点数运算、流水线设计、查找表、时序收敛、上板调试每一个环节都能让你学到东西。