1. 为什么要在FPGA里用CORDIC算三角函数1.1 一个真实的需求场景做数字信号处理或者通信基带的朋友大概率都遇到过这样的问题系统里需要实时计算sin和cos比如做数字下变频、正交解调、坐标旋转、相位检测甚至是电机控制里的Park变换。这时候你有几个选择。第一种直接调用厂商提供的IP核。Xilinx的CORDIC IP核确实好用配置几个参数就能出结果但问题是它是个黑盒你只知道输入输出中间怎么算的不清楚而且换到其他平台比如安路、高云、紫光同创就没有这个IP了移植起来很麻烦。第二种用查找表。把0到90度的sin值提前算好存到ROM里用的时候查表加插值。这个方法简单粗暴但精度和资源是一对矛盾——想要精度高ROM就得大Block RAM消耗惊人想要省资源精度又不够看。第三种用泰勒展开或者多项式逼近。理论上可行但阶数高了乘法器不够用阶数低了精度又不行而且定点数运算里截断误差很难控制。第四种就是CORDIC算法。它最大的好处是只用移位和加法就能算出三角函数不需要乘法器不需要大容量ROM精度可以通过迭代次数灵活控制。对于FPGA这种逻辑资源丰富但乘法器DSP Slice有限的器件来说CORDIC简直是量身定做的方案。我这次用的板卡是EGo1上面是Xilinx Artix-7系列的芯片。选它的原因很简单学校实验室里最常见价格便宜资源够用而且配套的Vivado工程模板很成熟。这篇文章就把我从算法推导到上板验证的完整过程拆开来讲包括代码、参数计算、仿真结果和踩过的坑。1.2 CORDIC旋转模式的本质是什么CORDIC的全称是Coordinate Rotation Digital Computer坐标旋转数字计算机。名字听着唬人核心思想其实特别朴素。想象你在一个二维平面上有一个点(x, y)你想把它旋转一个角度θ。按照正常的旋转矩阵x x·cosθ - y·sinθ y x·sinθ y·cosθ这个公式里出现了cosθ和sinθ问题来了——我要算sin和cos结果公式里又需要sin和cos这不是死循环了吗CORDIC的聪明之处在于它把旋转角度θ拆解成一系列预先定义好的小角度的累加。这些小角度满足一个特殊条件tan(θi) 2^(-i)。也就是说第i次旋转的角度θi arctan(2^(-i))。这样一来旋转矩阵里的cosθi和sinθi就可以用2的负幂次来表示cosθi 1 / sqrt(1 2^(-2i)) sinθi 2^(-i) / sqrt(1 2^(-2i))代入旋转公式提取出公共因子1/sqrt(1 2^(-2i))剩下的部分就变成了x x - y·di·2^(-i) y y x·di·2^(-i)其中di是方向因子取1或者-1代表这次旋转是逆时针还是顺时针。乘2^(-i)在硬件里就是右移i位一个移位操作就搞定了连乘法器都不需要。那个被提取出来的公共因子1/sqrt(1 2^(-2i))每次迭代都会产生一个n次迭代之后总的缩放因子K是K ∏ 1/sqrt(1 2^(-2i)) ≈ 1.646760258这个值是固定的跟旋转角度无关只跟迭代次数有关。所以实际使用的时候要么在最后把结果除以K要么在初始化的时候把输入预先乘以1/K约等于0.607252935。1.3 旋转模式的工作流程CORDIC有两种工作模式旋转模式和向量模式。旋转模式解决的是已知角度求sin和cos的问题向量模式解决的是已知坐标求幅值和相位的问题。这里我们只聊旋转模式。旋转模式的输入是初始坐标(x0, y0)和目标角度z0。通常我们把x0设成1/K也就是0.607252935y0设成0z0设成你想要计算的角度。然后开始迭代每一步我们看当前剩余角度z还剩多少。如果z 0说明还没转够di取1继续逆时针转如果z 0说明转多了di取-1往回顺时针转。每次迭代的角度是固定的arctan(2^(-i))这个值可以提前算好存成一张小表。迭代n次之后z趋近于0此时的x和y就是cos(z0)和sin(z0)。这里有个关键点CORDIC能覆盖的角度范围是有限的。因为所有arctan(2^(-i))加起来当i从0到无穷大时总和约为99.88度。也就是说CORDIC旋转模式只能处理-99.88度到99.88度之间的角度。超出这个范围怎么办利用三角函数的周期性把角度预先折叠到[-π/4, π/4]或者[-90°, 90°]范围内算完之后再根据象限恢复符号。这一步叫角度预处理是实际工程中必不可少的一环。2. 算法参数计算与定点数设计2.1 迭代次数怎么选迭代次数直接决定了精度和资源消耗。理论上每迭代一次精度大约提高1个比特。如果你需要16位精度那至少得迭代16次。但实际工程中由于定点数的量化误差和截断误差通常需要比目标精度多迭代2到3次。我这次的目标是16位输出精度所以选了16次迭代。实测下来16次迭代的误差在±2个LSB以内对于大多数信号处理应用已经足够了。如果你做的是高精度测量或者需要24位输出那就得迭代20次以上资源消耗会明显增加。这里给一个经验公式迭代次数n ≈ 目标精度位数 log2(n)。后面那个log2(n)是补偿缩放因子K的累积误差。比如目标16位n取16时log2(16)416420但实际测试16次就够了因为K的误差在16次迭代后已经小于1个LSB了。2.2 角度表的定点化CORDIC的核心是一张角度表存储每次迭代对应的arctan(2^(-i))。这张表用定点数表示格式很关键。我采用的是Q2.14格式也就是2位整数位加14位小数位总共16位。为什么选这个格式因为最大的角度arctan(1) 45度用弧度表示是0.785398需要至少1位整数位。但考虑到角度预处理后角度范围在[-90°, 90°]用弧度表示最大约1.57需要2位整数位Q2格式能表示-2到1.999...的范围。14位小数位对应精度为2^(-14) ≈ 6.1e-5弧度对于16位输出来说足够了。角度表的具体数值以弧度为单位Q2.14格式的十进制值迭代序号iarctan(2^(-i)) 弧度Q2.14定点值十进制Q2.14定点值十六进制00.785398163128680x324410.46364760975960x1DAC20.24497866340140x0FAE30.12435499520380x07F640.06241881010230x03FF50.0312398335120x020060.0156237292560x010070.0078123411280x008080.003906230640x004090.001953123320x0020100.000976562160x0010110.00048828180x0008120.00024414140x0004130.00012207020x0002140.00006103510x0001150.00003051810x0001注意最后两行的值都是1这是因为Q2.14格式的精度限制。从i14开始arctan(2^(-14))已经小于1个LSB了所以直接取1或者0都行。实际工程中如果迭代次数超过14次后面的角度值对结果影响微乎其微可以省略。2.3 缩放因子K的处理前面提到CORDIC迭代过程中会引入一个固定的缩放因子K ≈ 1.646760258。处理方式有两种方案一预先缩放输入。把初始x0设为1/K ≈ 0.607252935y0设为0。这样迭代完直接得到cos和sin不需要额外处理。这个方案的优点是输出就是最终结果缺点是初始值是个无理数定点化会引入误差。方案二后缩放输出。初始x0设为1y0设为0迭代完得到的是K·cos和K·sin最后再乘以1/K。这个方案的优点是初始值简单缺点是最后需要一次乘法。我选的是方案一因为EGo1上的DSP资源有限能省一个乘法器就省一个。1/K的Q2.14定点值是0.607252935 × 16384 ≈ 9949十六进制是0x26DD。2.4 数据位宽的选择位宽的选择是个权衡。位宽太窄精度不够位宽太宽资源浪费。我的设计是输入角度16位Q2.14内部迭代数据20位Q4.16输出16位Q1.15。为什么内部用20位因为迭代过程中会有累积误差中间多留4位作为保护位Guard Bits可以有效减少截断误差对最终结果的影响。输出的时候截取高16位低4位直接丢弃。为什么输出用Q1.15因为sin和cos的值域是[-1, 1]用1位符号位加15位小数位正好。Q1.15格式能表示的最小值是2^(-15) ≈ 3.05e-5对于16位精度来说刚好。3. Verilog代码实现与关键细节3.1 顶层模块设计顶层模块负责角度预处理、CORDIC迭代和输出后处理。接口很简单输入时钟、复位、角度值输出sin和cos。module cordic_sin_cos ( input wire clk, input wire rst_n, input wire start, input wire [15:0] angle_in, // Q2.14格式范围[-pi, pi] output reg done, output reg [15:0] sin_out, // Q1.15格式 output reg [15:0] cos_out // Q1.15格式 );角度预处理模块负责把输入角度折叠到[-π/4, π/4]范围内并记录象限信息。这一步很关键因为CORDIC只能处理[-99.88°, 99.88°]的角度而输入可能是任意角度。// 角度预处理将角度折叠到[-pi/4, pi/4] // 输入范围[-pi, pi]Q2.14格式 // pi/4的Q2.14值是0.785398 * 16384 12868 // pi/2的Q2.14值是1.570796 * 16384 25736 // pi的Q2.14值是3.141593 * 16384 51472 reg [1:0] quadrant; reg [15:0] angle_folded; reg sign_sin, sign_cos; always (*) begin if (angle_in[15]) begin // 负角度 // 取绝对值 if (angle_in 16h8000) begin // 角度在[-pi, -pi/2] quadrant 2b11; angle_folded 16h0000 - angle_in - 16h3FFF; // 加pi/2 sign_sin 1b1; sign_cos 1b0; end else begin // 角度在[-pi/2, 0] quadrant 2b10; angle_folded 16h0000 - angle_in; sign_sin 1b1; sign_cos 1b1; end end else begin if (angle_in 16h3FFF) begin // 角度在[pi/2, pi] quadrant 2b01; angle_folded angle_in - 16h3FFF; // 减pi/2 sign_sin 1b0; sign_cos 1b1; end else begin // 角度在[0, pi/2] quadrant 2b00; angle_folded angle_in; sign_sin 1b0; sign_cos 1b0; end end end这段预处理逻辑看起来有点绕但核心思想很简单把任意角度映射到第一象限的[0, π/4]范围内然后根据原始象限恢复符号。具体来说第一象限[0, π/2]直接算sin和cos都为正第二象限[π/2, π]角度减π/2算出来的sin就是coscos就是-sin第三象限[-π, -π/2]角度加π/2取反符号都取反第四象限[-π/2, 0]角度取反sin取反cos不变3.2 CORDIC迭代核心迭代核心是一个状态机每个时钟周期完成一次迭代。16次迭代需要16个时钟周期加上预处理和后处理总共约20个周期出结果。// CORDIC迭代核心 reg [4:0] iter_cnt; reg signed [19:0] x, y, z; reg signed [19:0] x_next, y_next, z_next; // 角度表Q2.14格式扩展到20位有符号数 wire signed [19:0] atan_table [0:15]; assign atan_table[0] 20sd12868; assign atan_table[1] 20sd7596; assign atan_table[2] 20sd4014; assign atan_table[3] 20sd2038; assign atan_table[4] 20sd1023; assign atan_table[5] 20sd512; assign atan_table[6] 20sd256; assign atan_table[7] 20sd128; assign atan_table[8] 20sd64; assign atan_table[9] 20sd32; assign atan_table[10] 20sd16; assign atan_table[11] 20sd8; assign atan_table[12] 20sd4; assign atan_table[13] 20sd2; assign atan_table[14] 20sd1; assign atan_table[15] 20sd1; // 迭代方向判断 wire di ~z[19]; // z为正时di1z为负时di0 always (posedge clk or negedge rst_n) begin if (!rst_n) begin iter_cnt 5d0; x 20sd9949; // 1/K的Q4.16定点值 y 20sd0; z 20sd0; done 1b0; end else if (start) begin if (iter_cnt 5d0) begin // 初始化 x 20sd9949; y 20sd0; z {angle_folded[15], angle_folded, 4b0000}; // 扩展到20位 iter_cnt 5d1; done 1b0; end else if (iter_cnt 5d16) begin // 迭代计算 if (di) begin x_next x - (y (iter_cnt - 1)); y_next y (x (iter_cnt - 1)); z_next z - atan_table[iter_cnt - 1]; end else begin x_next x (y (iter_cnt - 1)); y_next y - (x (iter_cnt - 1)); z_next z atan_table[iter_cnt - 1]; end x x_next; y y_next; z z_next; iter_cnt iter_cnt 1b1; end else begin // 迭代完成输出结果 sin_out {y[19], y[18:4]}; // 截取高16位 cos_out {x[19], x[18:4]}; done 1b1; iter_cnt 5d0; end end end这里有几个细节值得展开说。移位操作的处理。y (iter_cnt - 1)是算术右移对于有符号数来说右移会保留符号位。但要注意当iter_cnt1时移位量是0也就是不移位这是正确的因为第一次迭代的旋转角度是arctan(1)45度对应的因子是2^01。符号位的扩展。角度表的值是Q2.14格式的16位数但内部数据是20位所以需要符号扩展到20位。在Verilog里20sd12868这种写法会自动处理符号扩展。迭代次数的边界。当iter_cnt16时执行最后一次迭代然后iter_cnt变成17进入else分支输出结果。这里要注意iter_cnt是5位宽最大能表示31所以不会溢出。3.3 输出后处理与象限恢复迭代完成后根据预处理阶段记录的象限信息恢复sin和cos的符号。// 象限恢复 always (posedge clk or negedge rst_n) begin if (!rst_n) begin sin_final 16sd0; cos_final 16sd0; end else if (done) begin case (quadrant) 2b00: begin // 第一象限 sin_final sin_out; cos_final cos_out; end 2b01: begin // 第二象限 sin_final cos_out; cos_final -sin_out; end 2b10: begin // 第三象限 sin_final -sin_out; cos_final -cos_out; end 2b11: begin // 第四象限 sin_final -cos_out; cos_final sin_out; end endcase end end这里有个容易搞混的地方第二象限的sin是正的cos是负的。但我们的预处理是把角度减了π/2所以算出来的sin实际上是cos(θ-π/2) sin(θ)cos实际上是cos(θ-π/2) -sin(θ)。所以恢复的时候要对应交换。4. EGo1上板验证与实测结果4.1 硬件连接与约束文件EGo1板卡上的时钟是100MHz我用的是板载的50MHz晶振经过PLL倍频得到的。输入角度通过拨码开关或者串口输入输出通过LED或者数码管显示。为了验证方便我直接用Vivado的ILA集成逻辑分析仪抓取内部信号。约束文件的关键部分# 时钟约束 create_clock -period 20.000 -name sys_clk [get_ports clk] # 输入输出延迟约束 set_input_delay -clock sys_clk -max 2.0 [get_ports angle_in*] set_output_delay -clock sys_clk -max 2.0 [get_ports sin_out*] set_output_delay -clock sys_clk -max 2.0 [get_ports cos_out*]EGo1的引脚分配要注意拨码开关和LED的引脚在板卡手册里都有这里不展开。重点说一下时序约束CORDIC的迭代逻辑是组合逻辑加寄存器关键路径在移位和加法器上。100MHz时钟下16位加法器的延迟大约2-3ns加上布线延迟时序应该能过。如果时序不收敛可以把迭代拆成两级流水线代价是延迟增加一个周期。4.2 仿真验证仿真用的是Vivado自带的Simulator。Testbench里遍历了0到2π的所有角度步进π/180也就是1度对比CORDIC输出和MATLAB计算的参考值。// Testbench关键部分 initial begin rst_n 0; start 0; #100 rst_n 1; for (integer i 0; i 360; i i 1) begin angle_in i * 91; // 1度对应的Q2.14值约为91 start 1; #20 start 0; #500; // 等待计算完成 $display(Angle%d, sin%d, cos%d, i, sin_out, cos_out); end $finish; end仿真结果和MATLAB对比最大误差出现在45度附近约为±2个LSB。这个误差主要来自两个方面一是角度表的定点化误差二是迭代过程中的截断误差。对于16位输出来说2个LSB的误差相当于0.006%的相对误差完全可接受。4.3 资源消耗分析在Artix-7 XC7A35T上综合实现后的资源报告资源类型使用量可用量利用率LUT312208001.5%FF256416000.6%DSP0900%Block RAM0500%DSP使用量为0这是CORDIC最大的优势。整个设计只用了312个LUT和256个触发器对于Artix-7来说简直是九牛一毛。这意味着你可以在同一个芯片里实例化几十个CORDIC核并行处理多路信号。对比一下如果用Xilinx的CORDIC IP核默认配置下会消耗1个DSP48和若干LUT。虽然IP核的精度和速度可能更好但资源消耗也更高。对于成本敏感或者需要大量并行通道的应用手写CORDIC的优势很明显。4.4 实测波形与精度分析用ILA抓取的波形显示从start信号拉高到done信号拉高总共用了18个时钟周期100MHz下是180ns。这个延迟对于大多数实时信号处理应用来说完全可以接受。精度方面我做了两组测试测试一单点精度。输入角度45度Q2.14值12868理论sin和cos都是0.70710678。CORDIC输出sin0x5A82十进制23170cos0x5A82换算成浮点数是23170/327680.707092误差约为1.5e-5。测试二全角度扫描。0到360度每度一个点最大误差出现在接近0度和90度的位置约为±2个LSB。这个误差分布符合CORDIC算法的理论预期因为角度表的精度在两端最差。5. 常见问题与排查技巧实录5.1 输出一直是0或者不变这是新手最常遇到的问题。排查思路按顺序来第一步检查start信号。CORDIC核心是边沿触发的如果start一直为高状态机可能卡在某个状态。用ILA抓一下start和iter_cnt确认状态机在正常跳转。第二步检查复位信号。rst_n如果是低电平所有寄存器都被复位输出自然是0。确认复位释放的时机是否正确。第三步检查角度输入。如果angle_in一直是0那算出来的sin是0cos是1/K×K1这是正确的结果。但如果angle_in有变化而输出不变那可能是预处理模块的逻辑有问题。第四步检查时钟。EGo1的板载晶振是50MHz还是100MHz如果约束文件写错了时钟根本就没跑起来。5.2 精度不达标精度问题通常有三个来源角度表精度不够。如果你用的是Q1.15格式存角度表那精度只有2^(-15)≈3e-5弧度对于16位输出来说勉强够用。建议用Q2.14或者Q3.13多留一位整数位。迭代次数不够。前面说过16次迭代对应约16位精度。如果你需要20位精度至少得迭代20次。但注意迭代次数超过16次后角度表的值会变得非常小Q2.14格式已经无法表示需要改用Q3.13或者更高精度的格式。截断误差累积。每次迭代的移位操作都会丢弃低位16次迭代下来误差会累积。解决办法是在内部数据里多留几位保护位比如用24位内部数据输出时截取高16位。5.3 时序不收敛100MHz时钟下16位加法器的延迟加上布线延迟关键路径可能在5-6ns左右理论上能过。但如果你的设计里还有其他逻辑或者布线拥塞时序就可能不收敛。解决办法有两个一是降低时钟频率比如降到50MHz时序压力立刻减半二是插入流水线寄存器把一次迭代拆成两级第一级算移位第二级算加法。代价是延迟增加但吞吐量不变。5.4 象限判断错误这个问题很隐蔽因为大部分角度下输出看起来是对的只有在特定象限才会出错。排查方法是输入一个第二象限的角度比如135度理论sin是0.707cos是-0.707。如果输出sin-0.707cos0.707那就是象限恢复的逻辑写反了。我的经验是预处理和恢复的逻辑要成对设计。预处理时怎么折叠的恢复时就怎么展开最好画个表格对照着写。5.5 常见问题速查表现象可能原因排查方法解决方案输出恒为0复位未释放查rst_n波形确认复位释放时机输出恒为0start信号无脉冲查start波形确保start有上升沿输出不变角度输入未更新查angle_in确认输入数据变化精度差迭代次数不够增加迭代次数至少16次精度差角度表精度低检查Q格式改用Q2.14或更高时序不收敛关键路径太长查时序报告降频或插流水线象限错误恢复逻辑写反测第二象限角度对照表格修正资源超限迭代次数太多查资源报告减少迭代或复用5.6 几个实操心得心得一先用MATLAB验证算法。在写Verilog之前我习惯先用MATLAB或者Python把CORDIC算法跑一遍确认角度表、缩放因子、迭代逻辑都正确。这样上板的时候心里有底出了问题也知道是算法问题还是硬件问题。心得二ILA是调试神器。Vivado的ILA可以实时抓取内部信号比仿真更接近真实情况。我一般会在iter_cnt、x、y、z这几个关键信号上挂ILA一眼就能看出状态机卡在哪里。心得三角度预处理用查找表更省事。如果角度范围固定可以直接用查找表做预处理比一堆if-else清晰得多。比如把0到2π分成1024个点每个点对应一个象限和折叠后的角度查表就行。心得四输出加一级寄存器。CORDIC的输出直接连到外部引脚可能会有毛刺加一级寄存器打一拍输出会稳定很多。这一拍延迟对于大多数应用来说无所谓。心得五注意有符号数的移位。Verilog里的是算术右移会保留符号位这是正确的。但如果你不小心用了那就是逻辑右移负数会变成正数结果全错。这个坑我踩过调试了一下午才发现。6. 扩展应用与性能优化方向6.1 多通道并行处理CORDIC核的资源消耗很低312个LUT和256个FF在Artix-7上可以轻松实例化几十个。如果你需要同时处理多路信号比如8通道的正交解调直接复制8个CORDIC核就行资源完全够用。但要注意多个CORDIC核共享角度表可以进一步节省资源。角度表是常量综合工具会自动把它优化成查找表或者常量逻辑不会重复消耗Block RAM。6.2 流水线设计提升吞吐量当前的迭代式设计是每个时钟周期完成一次迭代16次迭代需要16个周期。如果你需要更高的吞吐量可以把迭代展开成流水线每一级流水线完成一次迭代。这样每个时钟周期都能出一个结果代价是延迟从16个周期变成16级流水线的延迟还是16个周期但吞吐量变成1。流水线设计的资源消耗会明显增加因为每一级都需要独立的寄存器和加法器。但对于高速应用来说这个代价是值得的。6.3 与其他算法的对比算法精度资源消耗速度适用场景CORDIC中高低无DSP中资源受限、多通道查找表中高Block RAM高精度要求不高泰勒展开高中DSP高单通道、高精度厂商IP核高中高快速开发CORDIC的核心竞争力在于零DSP消耗和灵活的精度控制。如果你的FPGA里DSP资源紧张或者需要大量并行通道CORDIC是最优解。6.4 进一步优化的思路角度表压缩。16个角度值里从i8开始后面的值都是2的幂次可以用移位生成不需要存储。这样角度表只需要存前8个值节省一半的查找表资源。自适应迭代次数。对于小角度输入前几次迭代的角度值很大可以跳过一些迭代。比如输入角度小于10度时第一次迭代的45度旋转可以直接跳过从第二次开始。这样可以减少迭代次数提高速度。混合精度。前几次迭代用高精度后几次迭代用低精度因为后面的角度值很小对结果影响不大。这样可以在保证精度的同时减少位宽节省资源。我在实际项目中用过角度表压缩的技巧把16个条目压缩到8个LUT消耗从312降到280左右效果不算特别明显但对于资源极度紧张的场景还是有用的。自适应迭代次数更适合角度范围已知且分布不均匀的场景比如通信系统里相位通常集中在某些区间。最后分享一个小技巧如果你用的是Vivado可以在综合设置里把CORDIC模块的-flatten_hierarchy设为none这样综合工具不会把层次打平方便你在网表里定位问题。这个设置对于调试复杂设计特别有用。