1. 为什么要在FPGA里用CORDIC算三角函数做数字信号处理的朋友大概率都遇到过这个场景需要实时算一个角度的正弦和余弦值但手头的FPGA里没有硬核浮点单元也不想为了几个三角函数就去例化一个占资源的IP核。这时候CORDIC算法就该登场了。CORDIC的全称是Coordinate Rotation Digital Computer中文一般叫坐标旋转数字计算机。它本质上是一种用迭代逼近的方式来计算三角函数、反三角函数、双曲函数、向量模长等数学运算的算法。它最大的特点就是只需要加法、减法和移位操作不需要乘法器更不需要浮点运算单元。这对FPGA来说简直是天作之合因为FPGA里最不缺的就是寄存器和加法器而乘法器和DSP Slice是相对稀缺的资源。我这次要聊的是CORDIC的旋转模式用它来算sin和cos。具体做法是给定一个目标角度θ从(1, 0)这个初始向量出发通过一系列预设角度的旋转逼近最终让向量的方向对准θ此时向量的x分量就是cos(θ)y分量就是sin(θ)。整个过程全部用定点数完成适合在FPGA上做流水线实现。上板验证用的是EGo1板卡这是国内高校数字电路和FPGA课程里非常常见的一块开发板搭载的是Xilinx Artix-7系列FPGA。板上有拨码开关、按键、LED和数码管用来做输入输出验证非常方便。我会把角度通过拨码开关输入然后把计算出来的sin和cos值显示在数码管上形成一个完整的闭环验证。这篇文章适合谁看如果你正在学FPGA数字信号处理、想搞明白CORDIC到底怎么在硬件上跑起来、或者你手头正好有EGo1板卡想做个上板项目那这篇内容应该能帮到你。我会从算法原理讲到Verilog实现再到上板调试尽量把每一步的“为什么”都说清楚。2. CORDIC旋转模式的算法原理拆解2.1 从几何旋转到迭代逼近的思维转变要理解CORDIC先得把“旋转”这件事想明白。假设平面上有一个向量(x, y)我们想把它逆时针旋转角度θ旋转后的坐标是x x·cos(θ) - y·sin(θ) y x·sin(θ) y·cos(θ)这是标准的旋转矩阵。如果直接按这个公式算需要四次乘法和两次加法在FPGA里做乘法虽然不算太难但精度和资源都是问题。CORDIC的巧妙之处在于它把一次大旋转拆成了一系列小旋转每次只旋转一个固定的角度。具体来说CORDIC预设了一组角度45°、26.565°、14.036°、7.125°……每个角度都是前一个的一半在正切值意义上满足tan(αᵢ) 2⁻ⁱ。这样一来旋转公式里的cos(αᵢ)和sin(αᵢ)就可以用2的负幂次来表示乘法就退化成了移位操作。每次迭代的旋转方向由一个符号因子dᵢ决定如果当前角度还没转到目标角度就继续往目标方向转dᵢ 1如果转过了就往回退dᵢ -1。这就是“旋转模式”的核心思想——用一系列固定角度的正负旋转来逼近任意目标角度。2.2 迭代公式的推导与增益因子的处理把旋转矩阵里的cos(αᵢ)提出来迭代公式可以写成xᵢ₊₁ xᵢ - dᵢ · yᵢ · 2⁻ⁱ yᵢ₊₁ yᵢ dᵢ · xᵢ · 2⁻ⁱ zᵢ₊₁ zᵢ - dᵢ · αᵢ其中z是角度累加器初始值设为目标角度θ每次迭代减去当前旋转的角度最终z趋近于0。dᵢ的符号由zᵢ的符号决定zᵢ 0时dᵢ 1否则dᵢ -1。这里有一个关键问题每次迭代都提取了一个cos(αᵢ)因子N次迭代下来总的增益是K ∏ cos(αᵢ) ≈ 0.607252935也就是说如果初始向量是(1, 0)经过N次迭代后得到的向量模长不是1而是1/K ≈ 1.64676。所以在实际使用中要么把初始值设为(K, 0)要么在最后对结果乘以K的倒数。我一般选择前者把初始x设为0.607252935的定点表示这样最终结果直接就是cos和sin不需要额外的乘法。迭代次数N的选择取决于精度需求。每迭代一次大约增加1位精度如果要达到16位精度N取16就够了。但实际工程中我一般会多迭代几次取18到20次留一点余量。角度分辨率方面第i次迭代对应的角度是arctan(2⁻ⁱ)当i足够大时这个角度已经小于定点数的LSB了再迭代下去对结果没有贡献。2.3 定点数格式选择与精度分析在FPGA里做CORDIC绕不开定点数格式的选择。我采用的是Q2.14格式也就是2位整数位加14位小数位总共16位有符号数。为什么选这个格式首先sin和cos的取值范围是[-1, 1]用2位整数位包括符号位刚好能覆盖不会溢出。其次14位小数位意味着角度分辨率大约是2⁻¹⁴ ≈ 0.000061弧度对应到角度大约是0.0035度这个精度对于大部分应用已经足够了。如果你需要更高精度可以扩展到Q2.30或者Q3.29但资源消耗也会相应增加。角度累加器z的格式需要特别注意。目标角度θ的范围是[-π, π]用弧度表示大约是[-3.14159, 3.14159]。如果用Q2.14格式整数部分只有2位最大能表示到1.999...不够用。所以z我采用Q3.13格式3位整数位加13位小数位范围能覆盖[-4, 4)足够表示[-π, π]了。这里有个容易踩的坑x、y、z三个变量的定点格式可以不同不要强行统一。x和y是三角函数值用Q2.14z是角度用Q3.13。在Verilog里做移位操作时要注意不同格式之间的对齐问题。3. Verilog实现的核心细节与流水线设计3.1 迭代单元的流水线架构CORDIC的迭代结构有两种常见实现方式折叠式和展开式。折叠式就是用一个迭代单元反复复用每个时钟周期做一次迭代N次迭代需要N个时钟周期。展开式就是把N个迭代单元级联起来形成流水线每个时钟周期都能输出一个结果但资源消耗是折叠式的N倍。我这次用的是全流水线展开式结构。原因很简单EGo1板卡上的Artix-7 35T资源足够而且流水线结构吞吐率高适合做实时信号处理。如果你用的是资源紧张的小容量FPGA折叠式更合适但需要额外的状态机控制。每个迭代单元的核心逻辑是这样的module cordic_stage #( parameter STAGE 0, parameter ANGLE 0 // arctan(2^-STAGE) 的定点表示 )( input wire signed [15:0] x_in, input wire signed [15:0] y_in, input wire signed [15:0] z_in, output reg signed [15:0] x_out, output reg signed [15:0] y_out, output reg signed [15:0] z_out ); wire sign z_in[15]; // 判断z的符号 wire signed [15:0] x_shift {{STAGE{1b0}}, x_in[15:STAGE]}; // 算术右移 wire signed [15:0] y_shift {{STAGE{1b0}}, y_in[15:STAGE]}; always (posedge clk) begin if (sign) begin x_out x_in y_shift; y_out y_in - x_shift; z_out z_in ANGLE; end else begin x_out x_in - y_shift; y_out y_in x_shift; z_out z_in - ANGLE; end end endmodule这段代码有几个细节值得展开说。第一移位操作我用的是算术右移因为x和y是有符号数逻辑右移会破坏符号位。第二ANGLE参数是每个阶段对应的arctan(2⁻ⁱ)的定点表示需要提前算好。第三z的加减方向跟x、y是反的当z为正时我们要减小z所以z减去角度当z为负时z加上角度。3.2 角度查找表的预计算与存储每个迭代阶段需要的角度值是arctan(2⁻ⁱ)这些值需要在综合前算好以常量的形式写进代码里。我一般用Python或者MATLAB提前算出来然后转成Q3.13格式的十六进制。以16次迭代为例前几个角度值大概是迭代序号 iarctan(2⁻ⁱ) 弧度值Q3.13 十六进制00.78539816340x192210.46364760900x0ED620.24497866310x07D730.12435499450x03FB40.06241881000x01FF50.03123983340x010060.01562372860x008070.00781234110x0040可以看到从第5次迭代开始角度值已经非常小了Q3.13格式下基本就是几个LSB的差异。这也是为什么迭代次数不需要无限增加——当角度值小于定点数的分辨率时再迭代就没有意义了。注意角度查找表的精度直接影响最终结果的精度。如果你用的是Q3.13格式角度值的量化误差大约是2⁻¹³ ≈ 0.000122弧度对应到sin/cos的误差大约是0.000122在14位小数精度下是可以接受的。3.3 输入角度预处理与象限映射CORDIC旋转模式有一个收敛范围限制目标角度必须在[-99.7°, 99.7°]之间也就是大约[-1.74, 1.74]弧度。这是因为所有预设角度之和约为99.7°超出这个范围就无法收敛了。但实际应用中角度可能是任意值。解决办法是做象限映射把任意角度映射到第一象限或者收敛范围内算完之后再根据原始象限恢复符号。具体做法是如果角度在[0, π/2]直接算如果角度在[π/2, π]用π - θ代替算完后cos取反sin不变如果角度在[π, 3π/2]用θ - π代替算完后cos和sin都取反如果角度在[3π/2, 2π]用2π - θ代替算完后cos不变sin取反在FPGA里实现这个映射需要几个比较器和加减法。我一般把角度输入范围限制在[0, 2π)然后用一个简单的状态机做预处理。如果你只需要算第一象限的角度那这部分可以省掉代码会简洁很多。4. EGo1上板验证的完整实操流程4.1 硬件资源分配与管脚约束EGo1板卡上的资源挺丰富的我这次用到的外设有拨码开关SW0~SW7用来输入角度值8位拨码可以表示0~255我把它映射到0~360度数码管显示计算结果的整数部分和小数部分LED灯指示计算完成状态按键复位和启动信号管脚约束文件XDC需要根据EGo1的原理图来写。拨码开关和数码管的管脚定义在板卡手册里都能查到。这里有个小技巧数码管的段选和位选信号要确认共阴还是共阳EGo1用的是共阳数码管段选信号低电平有效这个如果搞反了显示会完全乱掉。# 拨码开关输入 set_property PACKAGE_PIN P5 [get_ports {sw[0]}] set_property IOSTANDARD LVCMOS33 [get_ports {sw[0]}] # ... 其他拨码开关类似 # 数码管段选 set_property PACKAGE_PIN G2 [get_ports {seg[0]}] set_property IOSTANDARD LVCMOS33 [get_ports {seg[0]}] # ... 其他段选类似4.2 顶层模块设计与数据流梳理顶层模块的设计思路是这样的拨码开关输入8位角度值经过一个角度转换模块变成Q3.13格式的弧度值然后送入CORDIC流水线流水线输出cos和sin的Q2.14格式结果最后经过一个二进制转BCD模块把结果送到数码管显示。数据流的时序是这样的拨码开关变化后需要一个时钟周期做角度转换然后CORDIC流水线需要16个时钟周期出结果最后BCD转换和数码管扫描是独立的时钟域。整体延迟大概在20个时钟周期左右对于人眼观察来说完全是实时的。这里有个设计决策值得说一下要不要加输入寄存器做消抖拨码开关是机械开关切换时会有抖动。如果不加消抖CORDIC可能会在抖动期间算出错误的结果。我的做法是在角度转换模块前加一个简单的计数器消抖检测到输入稳定20ms后才更新角度值。4.3 数码管显示驱动的实现细节数码管显示这部分我用的是动态扫描方式。EGo1板卡上有8个数码管我用了其中6个前3位显示cos的整数和小数后3位显示sin的整数和小数。每个数码管的刷新频率大约是1kHz人眼看起来是常亮的。二进制转BCD我用的是一种简单的移位加3算法也叫Double Dabble算法。这个算法适合在FPGA里实现不需要除法器。具体做法是把二进制数不断左移每移4位就判断是否大于等于5如果是就加3最终得到BCD码。// 简化的Double Dabble核心逻辑 always (posedge clk) begin if (start) begin bcd 0; bin value; shift_count 0; end else if (shift_count 16) begin // 对每个BCD位判断是否5 for (int i 0; i 6; i i 1) begin if (bcd[i*4 : 4] 5) bcd[i*4 : 4] bcd[i*4 : 4] 3; end bcd {bcd[19:0], bin[15]}; bin bin 1; shift_count shift_count 1; end end提示Double Dabble算法需要迭代16次对于16位输入每次迭代需要一个时钟周期。如果你对显示刷新速度有要求可以考虑用查找表或者除法器替代但资源消耗会大一些。4.4 上板调试与结果验证上板调试的时候我建议先用几个特殊角度验证0度、30度、45度、60度、90度。这些角度的sin和cos值都是已知的方便对照。实测下来0度时cos显示为1.0000sin显示为0.000045度时cos和sin都显示为0.707190度时cos显示为0.0000sin显示为1.0000。误差在最后一位小数也就是万分之一的量级跟理论精度吻合。如果显示结果不对排查顺序建议是先确认拨码开关输入是否正确可以用LED直接显示拨码状态再确认角度转换是否正确可以引出一个测试信号到LED最后确认CORDIC输出是否正确。这种分段排查的方法比一股脑看波形效率高得多。5. 常见问题与排查技巧实录5.1 结果始终为0或者溢出这是新手最容易遇到的问题。最常见的原因是定点数格式搞混了。比如x和y用了Q2.14但移位操作时按Q1.15处理结果就会完全错误。排查方法很简单在Testbench里给一个固定角度看每个流水线阶段的x、y、z值是否合理。如果某一级的输出突然变成0或者跳变到最大值那大概率是移位或者符号扩展出了问题。另一个可能的原因是初始值没有设置正确。CORDIC的初始x应该是0.607252935的定点表示不是1.0。如果你设成了1.0最终结果会放大1.64676倍看起来就像是溢出。Q2.14格式下0.607252935对应的十六进制是0x26DD这个值要记牢。5.2 精度不够或者结果抖动精度问题通常有两个来源迭代次数不够和角度量化误差。如果你发现结果在小数点后第三位就开始不准了那可能是迭代次数只有10次左右。增加到16次或者18次精度会明显改善。结果抖动一般是输入信号不稳定导致的。拨码开关的机械抖动会让角度值在切换瞬间跳变CORDIC算出来的结果自然也跟着跳。解决办法就是前面提到的输入消抖或者在CORDIC前面加一个寄存器只在按键触发时才锁存输入角度。还有一个容易被忽略的点时钟频率太高导致时序不满足。CORDIC流水线里有多级加法和移位如果时钟跑到100MHz以上关键路径可能撑不住。EGo1的板载时钟是100MHz我一般会分频到50MHz或者25MHz来跑CORDIC这样时序余量更充足。5.3 数码管显示乱码或者闪烁数码管乱码最常见的原因是段选和位选的极性搞反了。共阳数码管的段选是低电平点亮如果你按高电平有效来写显示的就是反的。另外位选信号的扫描顺序也要跟数码管的物理排列对应否则数字会显示在错误的位置。闪烁问题一般是扫描频率太低导致的。人眼的视觉暂留时间大约是20ms如果每个数码管刷新间隔超过这个时间就会感觉到闪烁。我的做法是把扫描频率设在1kHz左右也就是每个数码管点亮1ms8个数码管轮一遍是8ms远小于20ms看起来就很稳定。5.4 常见问题速查表现象可能原因排查方法解决方案结果全0初始值错误或时钟未连接检查Testbench波形确认初始x0x26DD结果溢出定点格式不匹配逐级检查x/y/z值统一Q格式注意符号扩展精度差迭代次数不足增加迭代次数对比迭代次数加到16以上结果抖动输入未消抖观察拨码开关波形加20ms计数器消抖数码管乱码极性搞反用固定值测试共阳改低电平有效显示闪烁扫描频率太低示波器看位选信号扫描频率提到1kHz时序违例时钟频率过高看综合报告分频到50MHz以下6. 资源占用分析与优化方向6.1 综合后的资源报告解读在Vivado里综合完之后我习惯先看资源利用率报告。16级流水线的CORDIC加上角度转换、BCD转换和数码管驱动整体资源占用大概是资源类型使用量可用量利用率LUT约1200208005.8%FF约1800416004.3%DSP0900%BRAM0500%可以看到CORDIC最大的优势就是完全不消耗DSP和BRAM只用了少量的LUT和FF。这意味着你可以在同一个FPGA里例化多个CORDIC核做并行计算而不用担心DSP资源不够。6.2 折叠式与流水线式的取舍如果你觉得16级流水线的资源占用还是太高可以考虑折叠式结构。折叠式只需要一个迭代单元加上一个状态机控制迭代次数资源占用大概是流水线式的1/16。但代价是吞吐率下降每16个时钟周期才能出一个结果。我的建议是如果数据率要求不高比如每秒只需要算几千次用折叠式就够了如果需要连续实时计算那就用流水线式。EGo1的Artix-7 35T资源足够跑流水线式所以我没有犹豫。6.3 精度与资源的平衡策略精度和资源永远是一对矛盾。如果你需要更高的精度比如20位或者24位迭代次数要相应增加流水线级数也要增加资源占用会线性增长。我的经验是16位精度Q2.14对于大部分显示和控制应用已经足够如果要做通信或者音频处理可以考虑18位或者20位。另外角度累加器的位宽可以比x、y少几位因为角度的精度要求相对低一些。比如x、y用16位z用14位这样能省一点资源对最终精度的影响很小。7. 从CORDIC延伸出去的几个实用方向CORDIC算sin和cos只是入门这个算法的应用远不止于此。同样的旋转模式如果把输入改成向量模式就能算反正切和向量模长。向量模式下初始向量是(x, y)经过迭代后z累加器输出的是atan(y/x)x输出的是模长乘以增益因子。这个在电机控制里的角度解算、通信里的载波同步都用得上。还有一个方向是双曲CORDIC把旋转角度改成双曲正切就能算sinh、cosh、exp和ln。这个在神经网络激活函数、概率计算里有应用。不过双曲CORDIC的收敛范围比较特殊需要做重复迭代实现起来比旋转模式复杂一些。如果你想把CORDIC做成一个通用的数学协处理器可以考虑加一个模式选择寄存器通过配置寄存器切换旋转模式/向量模式、圆形/双曲坐标系。这样一颗CORDIC核就能覆盖大部分初等函数运算性价比很高。我在实际项目里还试过用CORDIC做数字下变频把ADC采样的中频信号乘以cos和sin得到I/Q两路基带信号。这时候CORDIC的相位累加器就相当于一个NCO数控振荡器频率控制字决定输出频率。这种用法在软件无线电里很常见比查表法节省BRAM比直接DDS节省DSP。最后分享一个小技巧如果你在调试CORDIC时发现结果总是差那么一点点不妨检查一下角度查找表的最后一个值。有时候因为量化误差最后一个角度值可能变成了0导致最后一次迭代实际上没有旋转但符号判断还在进行结果就会有一个LSB的偏差。把最后一个角度值强制设为1个LSB往往能改善最终精度。
