ARTICLE DETAIL

资讯详情

深耕编程入门与网站建设的一线实战洞察。

FPGA实战:CORDIC算法实现sin/cos,从原理到EGo1上板验证

FPGA实战:CORDIC算法实现sin/cos,从原理到EGo1上板验证 CORDIC这个算法我第一次接触的时候觉得它挺反直觉的——明明是要算三角函数它却全程只做加减和移位连一个乘法器都不用。但恰恰是这一点让它在FPGA上特别吃香。你想想FPGA里最宝贵的资源是什么是DSP Slice和Block RAM。如果为了算个sin/cos就把DSP占满那其他逻辑就别想跑了。CORDIC的精髓就在于把复杂的三角函数运算拆解成一连串固定角度的旋转逼近每次旋转只需要判断方向、做加减、再移位硬件开销小到令人发指。这篇内容我打算把基于FPGA的CORDIC旋转模式实现sin和cos运算这件事从头到尾讲透包括算法原理怎么理解、Verilog代码怎么写、流水线怎么设计、EGo1板卡上怎么验证、以及我在实际调试中踩过的那些坑。不管你是刚学FPGA的学生还是想找个靠谱的CORDIC实现方案拿来用的工程师这篇应该都能给你一些直接能抄的东西。1. 为什么CORDIC在FPGA上比查表法和泰勒展开更值得选1.1 三种常见sin/cos实现方案的硬碰硬对比在FPGA里算sin和cos能走的路其实就那么几条我把它们拉出来做个对比你就明白CORDIC的位置了。方案核心思路资源消耗精度速度适用场景查表法预存sin值到ROM大量Block RAM取决于表深度单周期读出精度要求低、角度范围小泰勒展开多项式逼近多个DSP乘法器高多周期软件实现或DSP资源充裕CORDIC迭代旋转逼近仅加减/移位/寄存器中等可调流水线后每周期一个结果资源受限、需要连续输出查表法的问题很直接你要16位精度角度分辨率假设0.01度那ROM深度就是36000个点每个点存16位差不多要吃掉576Kb的Block RAM。EGo1上那颗Artix-7的BRAM总共才多少这么一搞别的模块就别想用BRAM了。而且查表法还有个致命伤——它只能输出离散角度对应的值中间角度得靠插值插值又得加乘法器绕来绕去反而不划算。泰勒展开在数学上很优雅sin(x) x - x³/3! x⁵/5! - ...但你在FPGA上实现试试每一项都要乘法阶乘还得预计算角度大了收敛还慢。一个16位精度的sin运算少说也得三四个DSP乘法器级联延迟还不好控制。CORDIC就不一样了。它的核心迭代公式只有加法和移位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))其中d_i是旋转方向1或-12^(-i)就是右移i位。你看乘法变成了移位角度累加变成了加减法。整个运算下来一个乘法器都不需要全是LUT和寄存器就能搞定。这就是CORDIC在FPGA上的核心竞争力——用逻辑资源换DSP资源而逻辑资源在FPGA里通常比DSP富裕得多。1.2 旋转模式的几何直觉把目标角度转到零CORDIC有两种工作模式旋转模式Rotation Mode和向量模式Vectoring Mode。我们这里用的是旋转模式它的几何意义特别直观。想象你有一个初始向量(1, 0)对应角度0度。现在你想让它旋转到目标角度θ。CORDIC的做法不是一次性转过去而是一步一步地转每次转一个固定的角度arctan(2^(-i))但方向可以选——如果当前累积角度还没到θ就正转如果超过了就反转。这就像你开车找路不知道具体该打多少方向盘但你知道往左打一点和往右打一点哪个更接近目标。每次调整的角度越来越小最后无限逼近目标角度。迭代N次之后累积角度误差小于arctan(2^(-N))对于16位精度N取16就足够了。旋转模式的关键在于初始向量(1, 0)经过N次旋转后x分量就是cos(θ)y分量就是sin(θ)。但这里有个细节——每次旋转都会让向量长度略微增加因为旋转的模长因子是sqrt(1 2^(-2i))。所有迭代的模长因子乘积趋近于一个常数K ≈ 1.646760258。所以最终结果需要乘以1/K ≈ 0.607252935来进行幅度校正。这个校正可以在初始化时就把x设为1/K而不是1这样迭代完直接就是正确幅度。我在代码里就是这么干的省掉了一次乘法。1.3 角度累加表那些必须记住的固定角度CORDIC迭代中每次旋转的角度是固定的就是arctan(2^(-i))。这些角度值需要预先算好存在一个查找表里。对于16位迭代角度表如下以弧度为单位用定点数表示迭代次数iarctan(2^(-i)) 弧度值定点表示Q15格式00.78539816342573610.46364760901519320.2449786631802730.1243549945407440.0624188100204550.0312398334102460.015623728651270.007812341125680.003906230112890.001953122564100.000976562232110.000488281216120.00024414068130.00012207034140.00006103522150.00003051761注意看最后几行角度值已经小到定点表示只有1了。这意味着16位迭代的精度极限大概在2^(-15)弧度左右。如果你需要更高精度就得增加迭代次数同时角度表的位宽也要跟着增加。这里有个容易踩的坑角度表的定点格式必须和你的角度累加器位宽匹配。我一开始用Q15格式存角度表但角度累加器用了18位结果每次累加都丢精度最后输出误差大得离谱。后来统一成Q17格式才解决。2. Verilog实现从单周期迭代到全流水线架构2.1 单周期迭代版本先跑通再优化我建议你第一次实现CORDIC的时候先写一个单周期迭代的版本。就是一个时钟周期完成一次迭代16次迭代就是16个时钟周期出一个结果。虽然速度慢但逻辑简单容易调试。module cordic_single ( input wire clk, input wire rst_n, input wire start, input wire [17:0] angle_in, // Q17格式目标角度 output reg [17:0] cos_out, output reg [17:0] sin_out, output reg done ); // 角度查找表Q17格式 reg [17:0] atan_table [0:15]; initial begin atan_table[0] 18d102944; // arctan(1) * 2^17 atan_table[1] 18d60775; atan_table[2] 18d32109; atan_table[3] 18d16297; atan_table[4] 18d8180; atan_table[5] 18d4096; atan_table[6] 18d2048; atan_table[7] 18d1024; atan_table[8] 18d512; atan_table[9] 18d256; atan_table[10] 18d128; atan_table[11] 18d64; atan_table[12] 18d32; atan_table[13] 18d16; atan_table[14] 18d8; atan_table[15] 18d4; end reg [4:0] iter_cnt; reg [17:0] x_reg, y_reg, z_reg; reg busy; // 1/K的Q17表示K ≈ 1.646760258 localparam [17:0] INV_K 18d79623; // 0.607252935 * 2^17 always (posedge clk or negedge rst_n) begin if (!rst_n) begin iter_cnt 5d0; x_reg 18d0; y_reg 18d0; z_reg 18d0; busy 1b0; done 1b0; end else if (start !busy) begin // 初始化x 1/K, y 0, z 目标角度 x_reg INV_K; y_reg 18d0; z_reg angle_in; iter_cnt 5d0; busy 1b1; done 1b0; end else if (busy) begin if (iter_cnt 5d15) begin busy 1b0; done 1b1; cos_out x_reg; sin_out y_reg; end else begin iter_cnt iter_cnt 1b1; if (z_reg[17]) begin // z为负需要正转 x_reg x_reg - (y_reg iter_cnt); y_reg y_reg (x_reg iter_cnt); z_reg z_reg atan_table[iter_cnt]; end else begin // z为正需要反转 x_reg x_reg (y_reg iter_cnt); y_reg y_reg - (x_reg iter_cnt); z_reg z_reg - atan_table[iter_cnt]; end end end end endmodule这段代码有几个关键点需要解释。第一z_reg[17]是符号位判断因为角度用Q17格式最高位是符号位为1表示负数。第二是算术右移对于有符号数来说算术右移会保留符号位这正是我们需要的。第三初始化时x设为INV_K而不是1这样迭代完直接就是cos值不需要额外校正。但单周期版本有个问题每次迭代都要等一个时钟周期16次迭代就是16个周期。如果你需要连续计算多个角度吞吐率就很低。这时候就需要流水线架构。2.2 全流水线设计每个时钟周期出一个结果流水线的思路是把16次迭代展开成16级每级之间用寄存器隔开。这样虽然延迟还是16个周期但吞吐率变成了每个时钟周期都能接收新数据、输出新结果。module cordic_pipeline ( input wire clk, input wire rst_n, input wire [17:0] angle_in, output wire [17:0] cos_out, output wire [17:0] sin_out ); // 16级流水线寄存器 reg signed [17:0] x [0:16]; reg signed [17:0] y [0:16]; reg signed [17:0] z [0:16]; // 角度表 reg signed [17:0] atan_table [0:15]; initial begin atan_table[0] 18sd102944; atan_table[1] 18sd60775; atan_table[2] 18sd32109; atan_table[3] 18sd16297; atan_table[4] 18sd8180; atan_table[5] 18sd4096; atan_table[6] 18sd2048; atan_table[7] 18sd1024; atan_table[8] 18sd512; atan_table[9] 18sd256; atan_table[10] 18sd128; atan_table[11] 18sd64; atan_table[12] 18sd32; atan_table[13] 18sd16; atan_table[14] 18sd8; atan_table[15] 18sd4; end localparam signed [17:0] INV_K 18sd79623; integer i; always (posedge clk or negedge rst_n) begin if (!rst_n) begin for (i 0; i 16; i i 1) begin x[i] 18sd0; y[i] 18sd0; z[i] 18sd0; end end else begin // 第一级初始化 x[0] INV_K; y[0] 18sd0; z[0] angle_in; // 后续15级迭代 for (i 0; i 15; i i 1) begin if (z[i][17]) begin x[i1] x[i] - (y[i] i); y[i1] y[i] (x[i] i); z[i1] z[i] atan_table[i]; end else begin x[i1] x[i] (y[i] i); y[i1] y[i] - (x[i] i); z[i1] z[i] - atan_table[i]; end end end end assign cos_out x[15]; assign sin_out y[15]; endmodule流水线版本看起来代码量差不多但结构完全不同。每一级迭代的结果都锁存在寄存器里下一级直接读上一级的结果。这样综合出来的电路是一条长长的流水线每个时钟周期都能吞进一个新的角度值同时吐出一个16周期前输入的角度对应的sin/cos值。注意流水线版本里我用了reg signed因为算术右移对无符号数和对有符号数的行为不同。CORDIC迭代中x和y可能为负必须用有符号数才能保证右移正确。2.3 位宽选择与精度权衡18位到底够不够位宽选择是CORDIC实现中最需要动脑子的事情之一。位宽太窄精度不够位宽太宽资源浪费。我一般从以下几个维度来考虑角度精度如果你需要0.01度的角度分辨率那角度累加器的位宽至少要能表示2π/0.01 ≈ 628个刻度也就是10位。但CORDIC的精度不仅取决于角度位宽还取决于迭代次数。16次迭代的理论角度误差是arctan(2^(-16)) ≈ 0.000015弧度 ≈ 0.00087度远高于0.01度的要求。数据精度x和y的位宽决定了sin/cos的输出精度。每增加一位精度提升约6dB。18位有符号数的动态范围是-131072到131071对应Q17格式就是-1到1分辨率是2^(-17) ≈ 0.0000076。这个精度对于大多数应用足够了。内部增长CORDIC迭代过程中x和y的幅度会增长约1.6467倍。如果你初始x设为1/K那迭代过程中x和y的幅度始终在1附近不会溢出。但如果你初始x设为1那迭代完x和y的幅度会到1.6467需要额外留出2位保护位。我在EGo1上实测过18位位宽、16次迭代的配置sin/cos输出误差在±2个LSB以内也就是±0.000015左右。这个精度对于电机控制、信号生成、坐标变换这些应用完全够用。3. EGo1板卡上板验证从仿真到硬件的完整链路3.1 EGo1的硬件资源与引脚分配要点EGo1是Xilinx Artix-7系列的一款教学板卡具体型号是XC7A35T-1CSG324C。它的资源对于CORDIC验证来说绰绰有余33280个逻辑单元、1800Kb Block RAM、90个DSP Slice。我们只用LUT和寄存器就能实现CORDICDSP一个都不用占。上板验证需要用到几个外设时钟源、复位按键、以及输出显示。EGo1上有一个100MHz的晶振可以直接作为系统时钟。复位可以用板上的按键但按键有机械抖动需要做消抖处理。输出显示我建议用两种方式一是通过UART把sin/cos值传到电脑上看波形二是用板上的LED或数码管显示几个特定角度的结果。引脚分配是个细致活。EGo1的引脚约束文件需要根据原理图来写。时钟引脚是E3复位按键是C12UART的TX和RX分别是M5和N5。这些在EGo1的用户手册里都有但手册有时候版本对不上我建议你直接看板子背面的丝印或者用Vivado的引脚规划器对照原理图确认。# EGo1约束文件片段 set_property PACKAGE_PIN E3 [get_ports clk] set_property IOSTANDARD LVCMOS33 [get_ports clk] set_property PACKAGE_PIN C12 [get_ports rst_n] set_property IOSTANDARD LVCMOS33 [get_ports rst_n] set_property PACKAGE_PIN M5 [get_ports uart_tx] set_property IOSTANDARD LVCMOS33 [get_ports uart_tx]有个细节容易忽略EGo1的时钟是100MHz但CORDIC流水线在100MHz下可能时序紧张。如果你综合后发现时序不满足可以先用50MHz跑或者把流水线级数增加。我在100MHz下跑16级流水线时序余量大概还有0.5ns算是比较稳的。3.2 测试数据生成怎么验证输出是对的上板之前仿真验证是必须的。但仿真激励怎么写直接决定了你能不能发现潜在问题。我一般会构造三类测试数据第一类是特殊角度0度、30度、45度、60度、90度、180度、270度。这些角度的sin/cos值有解析解容易对照。比如45度时sin和cos都应该是0.7071对应Q17格式就是92682。第二类是边界角度接近0度和接近360度的值。这些角度容易暴露符号处理和溢出问题。比如角度为1度时sin值很小cos值接近1如果位宽不够sin值可能被截断成0。第三类是随机角度用MATLAB或Python生成一批随机角度算出理论sin/cos值和仿真输出对比。我一般会跑1000个随机点统计最大误差和平均误差。# Python生成测试激励 import numpy as np angles np.random.uniform(0, 2*np.pi, 1000) with open(test_angles.txt, w) as f: for a in angles: # 转换为Q17定点格式 fixed int(a / (2*np.pi) * 2**17) 0x3FFFF f.write(f{fixed:018b}\n)仿真的时候我用了一个简单的testbench每个时钟周期喂一个角度然后收集输出。对比的时候要注意流水线有16个周期的延迟所以输出要延迟16个周期再和输入对齐。3.3 实测波形分析那些仿真发现不了的问题仿真过了不代表上板就能跑。我在EGo1上实测的时候遇到了几个仿真阶段完全没暴露的问题。第一个问题是复位释放的时机。EGo1的复位按键没有专门的复位芯片释放时有个缓慢的上升沿。如果CORDIC状态机在复位释放的瞬间就开始工作可能因为复位信号还在阈值附近而导致状态机进入非法状态。我的解决办法是加一个复位同步器用两级触发器对复位信号做同步并且延迟几个周期再释放。第二个问题是时钟质量。EGo1的100MHz晶振是普通的有源晶振抖动比较大。CORDIC流水线对时钟抖动不敏感但如果你后面接了DAC输出模拟波形抖动就会体现在输出上。我建议在时钟输入后加一个BUFG并且如果条件允许用MMCM把时钟抖动滤一下。第三个问题是UART传输的误码。我一开始用115200波特率传sin/cos值发现偶尔会有几个数据点跳变。后来查出来是UART的波特率发生器在100MHz时钟下分频系数不是整数累积误差导致采样点偏移。改成921600波特率后问题消失因为分频系数更接近整数。4. 调试中踩过的坑与性能优化经验4.1 角度累加溢出的隐蔽问题CORDIC迭代中角度累加器z的位宽必须足够容纳目标角度加上所有atan值的和。最坏情况下如果目标角度接近π/2而每次迭代都往同一个方向累加z的绝对值可能超过π。如果你用Q17格式π对应的是411774而18位有符号数的范围是-131072到131071根本装不下。我一开始就踩了这个坑。角度输入范围是0到2πQ17格式下最大是823549但18位有符号数最大才131071。结果就是角度输入直接被截断输出完全不对。解决办法有两个一是增加角度累加器的位宽至少要到20位二是把角度输入范围限制在-π/2到π/2利用三角函数的周期性做预处理。我选择了第二种方案因为CORDIC本身在-π/2到π/2范围内收敛最好超出这个范围精度会下降。// 角度预处理把任意角度映射到[-π/2, π/2] always (*) begin case (angle_in[19:18]) // 假设输入是20位高2位表示象限 2b00: begin // 第一象限 angle_mapped angle_in[17:0]; cos_sign 1b0; sin_sign 1b0; end 2b01: begin // 第二象限 angle_mapped 18sd102944 - angle_in[17:0]; // π/2 - θ cos_sign 1b1; sin_sign 1b0; end // ... 第三、四象限类似 endcase end4.2 流水线时序收敛的几种手段16级流水线在100MHz下时序收敛我试了好几种方法才搞定。最直接的是在每级迭代之间插入额外的寄存器把组合逻辑路径打断。但这样会增加延迟从16周期变成32周期。另一种方法是优化移位器的实现。在Verilog里综合出来是一个桶形移位器对于可变移位量它会综合成多级MUX。但CORDIC里每级的移位量是固定的第i级移i位所以综合工具应该能优化成简单的连线。如果你发现时序报告里移位器是关键路径可以手动把它写成固定移位// 不要这样写综合器可能不知道移位量是固定的 y_shifted y i; // 可以这样写明确告诉综合器移位量 case (i) 0: y_shifted y; 1: y_shifted {y[17], y[17:1]}; 2: y_shifted {{2{y[17]}}, y[17:2]}; // ... endcase还有一种方法是降低时钟频率。EGo1的100MHz对于16级流水线来说确实有点紧张如果你不追求极致速度用50MHz跑会更稳。我实测50MHz下时序余量有2ns以上非常安全。4.3 资源占用实测LUT和寄存器的真实消耗我在Vivado里综合了流水线版本的CORDIC目标器件是XC7A35T-1CSG324C综合策略是默认的Vivado Synthesis。资源占用如下资源类型使用量可用量占比LUT1247208006.0%FF918416002.2%DSP0900%BRAM0500%你看DSP和BRAM一个都没用全是LUT和FF。1247个LUT对于XC7A35T来说只占了6%剩下的资源足够你跑其他逻辑。这就是CORDIC的魅力——用最廉价的逻辑资源实现三角函数运算。如果你把迭代次数从16增加到20LUT会增加到约1600个FF增加到约1150个精度提升到2^(-19)弧度。这个代价完全可以接受。5. 从CORDIC到实际应用几个可以直接复用的场景5.1 数字信号发生器中的正交本振生成CORDIC最直接的应用就是生成正交本振信号。在通信系统里你需要一对频率相同、相位差90度的正弦波来做IQ调制解调。用CORDIC的话你只需要一个相位累加器每个时钟周期累加频率控制字然后把累加结果作为角度输入CORDIC输出的cos和sin就是正交本振。这个方案的优点是频率切换非常快因为相位累加器是纯数字的改变频率控制字就能瞬间切换频率。而且CORDIC的流水线结构保证了每个时钟周期都能输出一对新的IQ样本采样率就是时钟频率。我在EGo1上试过用CORDIC生成1kHz到1MHz的正弦波通过UART传到电脑上用Python画波形频谱非常干净谐波抑制在-60dB以下。这个性能对于一般的信号发生器应用完全够用。5.2 电机控制中的Park变换与Clarke变换在电机控制里Park变换和Clarke变换都需要用到sin和cos。传统的做法是用查表法但查表法的精度受限于表深度而且占用大量ROM。用CORDIC的话你可以实时计算任意角度的sin/cos精度只取决于迭代次数。具体来说Park变换的公式是Id Iα * cos(θ) Iβ * sin(θ) Iq -Iα * sin(θ) Iβ * cos(θ)你只需要把CORDIC输出的cos和sin接到两个乘法器上就能完成Park变换。整个链路延迟只有CORDIC的流水线延迟加上乘法器的延迟对于电机控制来说完全在可接受范围内。5.3 坐标旋转与相位检测CORDIC的向量模式可以用来做相位检测但旋转模式也可以间接实现。如果你有两个正交信号I和Q想知道它们的相位角可以用CORDIC的向量模式。但如果你已经知道角度想生成对应的I和Q那就是旋转模式。我在一个项目里用CORDIC做过相位锁定环PLL的鉴相器。把输入信号和本地振荡信号做乘法得到I和Q然后用CORDIC向量模式算出相位误差再反馈给相位累加器。整个环路全部在FPGA里实现锁定时间小于1ms相位误差小于0.1度。6. 写在最后一些个人体会CORDIC这个算法看起来简单但真正写好、调通、上板稳定运行还是需要不少细节上的打磨。我最大的体会是不要一上来就追求流水线和高精度先用单周期版本把功能跑通用仿真确认算法正确再逐步优化架构。很多问题在单周期版本里很容易定位到了流水线里就变得很难调试。另一个体会是位宽的选择要留有余量。我一开始用16位发现精度不够改成18位后精度够了但角度累加器又溢出了最后角度用20位、数据用18位才稳定。每次改位宽都要重新算角度表、重新验证挺折腾的。所以如果你时间充裕建议一开始就把位宽定得宽一点比如数据20位、角度22位这样后面调整的空间更大。还有一点EGo1这块板子虽然资源不多但跑CORDIC绰绰有余。如果你手头有其他的Artix-7或Zynq板卡代码几乎不用改只需要调整引脚约束就行。CORDIC的可移植性非常好这也是它经久不衰的原因之一。
返回列表