ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

FPGA上CORDIC算法实现sin/cos:从数学推导到EGo1上板验证

FPGA上CORDIC算法实现sin/cos:从数学推导到EGo1上板验证 1. 为什么要在FPGA里用CORDIC算三角函数很多人第一次听到用FPGA算sin和cos脑子里蹦出来的第一个方案是查表法把0到90度的正弦值预先算好存进一块ROM里用的时候按角度查表再配合象限判断和线性插值。这个思路没错在精度要求不高、资源又比较紧张的场合查表法确实够用。但它有两个绕不开的硬伤一是精度和存储深度直接挂钩想要16位精度ROM的规模会迅速膨胀二是它本质上是个死方法角度分辨率被表长锁死想提高分辨率就得重新生成表、重新综合迭代成本高。CORDICCoordinate Rotation Digital Computer坐标旋转数字计算机走的是另一条路。它不存表而是通过一系列固定角度的旋转逼近把目标角度转出来。核心思想非常朴素把任意角度分解成一堆预先定好的、大小递减的基准角之和每次只判断往左转还是往右转转完累加。因为每次旋转的角度是固定的旋转矩阵里的cos项可以提出来当常数增益剩下的乘法全部退化成移位和加减法。这意味着整个运算只需要加法器、减法器、移位寄存器和一张极小的角度查找表非常适合在FPGA这种逻辑资源有限、又极度依赖并行和流水线的器件上实现。我这次做的是CORDIC的旋转模式Rotation Mode目标是输入一个角度θ输出对应的cosθ和sinθ。旋转模式和向量模式容易搞混这里先厘清旋转模式是已知角度求坐标向量模式是已知坐标求角度和模长。做sin/cos用的是旋转模式做arctan和求模用的是向量模式两者迭代公式一样只是初始条件和判断方向不同。选EGo1这块板子来上板验证理由很实际。EGo1是Xilinx Artix-7系列的教学板逻辑资源够跑一个中等规模的CORDIC核板载的数码管、LED、按键资源齐全特别适合把输入角度、输出结果这种交互式验证做出来。相比在仿真里看波形上板能看到真实的数码管显示对理解定点数、角度映射、象限处理这些细节帮助极大。这篇内容适合已经写过基本Verilog、懂一点时序逻辑、想入门FPGA数字信号处理的同学也适合做过查表法、想换个思路理解CORDIC的老手。2. CORDIC旋转模式的数学推导与定点化取舍2.1 从旋转矩阵到迭代公式二维平面里把一个点(x, y)绕原点旋转角度φ标准旋转矩阵是x x·cosφ - y·sinφ y x·sinφ y·cosφCORDIC的巧妙之处在于它不一次性转φ而是拆成n次小旋转每次转的角度是atan(2^(-i))i从0开始递增。为什么选这个角度序列因为tan(atan(2^(-i))) 2^(-i)正好是2的负整数次幂乘以它等价于右移i位硬件上零成本。把cosφ提出来第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或-1由剩余角度z_i的符号决定z_i大于0就往正方向转d_i1小于0就往负方向转d_i-1。z_i是还差多少角度没转完的累加器初始值就是目标角度θ。迭代足够多次后z_i趋近于0此时(x_n, y_n)就是旋转后的坐标。这里有个关键点必须讲清楚每次迭代都提了一个cos(atan(2^(-i)))出来n次迭代累积的增益是K ∏ cos(atan(2^(-i))) ≈ 0.607252935也就是说如果初始输入是(x_0, y_0) (1/K, 0)那么迭代结束后x_n cosθy_n sinθ。这个1/K ≈ 1.646760258是CORDIC里最重要的一个常数通常直接预置成初始x值省掉一次乘法。2.2 迭代次数、精度与位宽的三角关系迭代次数n决定了角度分辨率。第i次迭代能分辨的最小角度是atan(2^(-i))当i增大这个角度迅速变小。n次迭代后角度误差量级约为atan(2^(-(n-1)))。想要16位精度n取16基本够用想要更高精度n要相应增加但每增加一次迭代就多一级流水线资源和延迟都上升。位宽的选择更考验经验。定点数做CORDIC通常用Q格式比如Q1.15表示1位符号位加15位小数位。但CORDIC迭代过程中x和y会先被增益放大到约1.6467倍如果输入是满量程的1.0中间值会超过1定点数就溢出了。所以实际工程里有两个常见做法一是把初始x设成1/K让最终结果落在[-1,1]内中间过程最大值约1.0刚好不溢出二是留出保护位比如内部用Q2.16甚至Q3.17输出时再截断。我这次用的是16次迭代、内部18位定点1位符号2位整数15位小数输出取高16位。这个配置在Artix-7上综合出来的资源占用很舒服精度实测能到小数点后4位稳定数码管显示足够。2.3 角度预处理把任意角度压进收敛区间CORDIC旋转模式有个收敛范围限制只有当目标角度落在[-99.7°, 99.7°]约±π/2内迭代才能收敛。因为所有atan(2^(-i))之和的极限就是99.7°。超过这个范围z_i永远归不了零。解决办法是做象限折叠。把输入角度θ先规约到第一象限用三角函数的对称性还原原始角度范围规约操作cos还原sin还原[0, π/2]直接用θcossin[π/2, π]θ π - θ-cossin[π, 3π/2]θ θ - π-cos-sin[3π/2, 2π]θ 2π - θcos-sin这一步在硬件上就是几个比较器和加减法代价很小但少了它整个设计就只能在第一象限工作。很多网上的CORDIC代码只做第一象限上板一测就发现角度一大结果全错问题就出在这。3. Verilog实现从角度累加器到流水线输出3.1 顶层模块的接口设计先定接口。输入是系统时钟、复位、一个角度值用16位无符号表示0到2π即0到65535对应0到360度输出是16位有符号的cos和sin。为了让上板验证直观我额外加了一个角度步进逻辑用按键控制角度每次加固定值数码管轮流显示角度、cos、sin。module cordic_sincos ( input wire clk, input wire rst_n, input wire [15:0] angle_in, // 0~65535 映射 0~360度 output reg signed [15:0] cos_out, output reg signed [15:0] sin_out, output reg valid );角度用16位无符号表示整圈这个映射很实用360度对应655361度约等于182分辨率约0.0055度对数码管显示绰绰有余。用无符号整圈表示的好处是象限判断直接看高两位不用做浮点比较。3.2 象限折叠的硬件实现把16位角度的高两位拿出来判断象限低14位作为第一象限内的角度。这里有个细节第一象限的角度范围是0到π/2对应低14位的0到16383。但CORDIC的atan表是按弧度或按归一化角度建的需要统一量纲。我的做法是把角度统一用归一化到π/2的整数表示即16384对应90度。这样atan(2^(-i))对应的值可以预先算好存成一张16项的查找表// atan(2^-i) 归一化到 16384 90度 // atan(2^-0)45度 - 8192 // atan(2^-1)26.565度 - 4836 // ... 依次递减这张表是CORDIC唯一的存储开销16个16位数连一块小ROM都算不上直接写成case或者常量数组综合器会自动优化成组合逻辑。象限折叠的代码逻辑always (*) begin case (angle_in[15:14]) 2b00: begin // 第一象限 z_init {2b00, angle_in[13:0]}; cos_sign 1b0; sin_sign 1b0; end 2b01: begin // 第二象限 z_init {2b00, 14d16384 - angle_in[13:0]}; cos_sign 1b1; sin_sign 1b0; end 2b10: begin // 第三象限 z_init {2b00, angle_in[13:0] - 14d16384}; cos_sign 1b1; sin_sign 1b1; end 2b11: begin // 第四象限 z_init {2b00, 14d32768 - angle_in[13:0]}; cos_sign 1b0; sin_sign 1b1; end endcase end注意第三象限的规约是θ - π对应低14位减去16384第四象限是2π - θ对应32768减去低14位。这些数值关系在纸上画个单位圆推一遍就清楚了别死记。3.3 迭代核心16级流水线的写法CORDIC的迭代天然适合流水线。每一级做三件事根据z的符号决定旋转方向、更新x和y、更新z。16级就是16个always块或者一个generate循环。genvar i; generate for (i 0; i 16; i i 1) begin : cordic_stage always (posedge clk or negedge rst_n) begin if (!rst_n) begin x[i1] 0; y[i1] 0; z[i1] 0; end else begin if (z[i][17]) begin // z为负顺时针转 x[i1] x[i] (y[i] i); y[i1] y[i] - (x[i] i); z[i1] z[i] atan_table[i]; end else begin // z为正逆时针转 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 endgenerate这里有几个坑必须提醒。第一是算术右移对有符号数才能保持符号位千万别用。第二x和y的位宽要留够我内部用18位第i级右移i位后还要参与加减位宽不够会截断出错。第三z的符号判断看最高位18位有符号数最高位是符号位z[i][17]为1表示负数。初始值设置x[0] 18d39797即1/K × 2^15 ≈ 1.6467 × 32768y[0] 0z[0] 象限折叠后的角度。这个39797是CORDIC的魔法数字记不住就现场算1.646760258 × 32768 53961等等这里要小心量纲。如果内部是Q2.151.0对应32768那1/K对应53961。但如果x最终要落在[-1,1]初始x设成1/K中间过程最大值约1.0用Q2.15刚好。我实际用的是Q2.15初始x 53961。3.4 输出截断与符号还原16级流水线走完x[16]和y[16]就是第一象限的cos和sin但它们是Q2.15格式且被增益放大过——不对因为初始x已经除了增益所以最终结果就是标准值。取高16位作为输出再根据象限折叠时记录的cos_sign和sin_sign还原符号always (posedge clk or negedge rst_n) begin if (!rst_n) begin cos_out 0; sin_out 0; valid 0; end else begin cos_out cos_sign ? -x[16][17:2] : x[16][17:2]; sin_out sin_sign ? -y[16][17:2] : y[16][17:2]; valid 1b1; end end注意流水线延迟是16个时钟周期valid信号要相应打拍否则上层模块会在结果还没出来时就读走旧数据。这个延迟在仿真里一眼能看出来上板时如果数码管刷新逻辑没对齐会看到显示跳变别慌是时序没对齐。4. EGo1上板验证从约束文件到数码管显示4.1 引脚约束与时钟配置EGo1板载的晶振是100MHz但CORDIC流水线跑100MHz完全没问题甚至能跑到150MHz以上。不过数码管刷新和按键消抖不需要这么快我习惯分频出一个1kHz左右的慢时钟给显示逻辑CORDIC核还是用100MHz主时钟两者用valid信号握手。约束文件里最关键的是时钟引脚和数码管、按键的引脚。EGo1的数码管是共阳极段选低电平点亮位选也是低有效。按键需要消抖我用一个20位计数器做20ms消抖实测很稳。# 时钟约束 create_clock -period 10.000 -name sys_clk [get_ports clk] # 数码管段选a~g, dp set_property PACKAGE_PIN ... [get_ports {seg[0]}] # 位选 set_property PACKAGE_PIN ... [get_ports {an[0]}]引脚号这里不一一列了每块板子的约束文件不一样照着自己板子的原理图填。但有个经验数码管的段选线最好加限流电阻EGo1板载已经带了直接连就行。4.2 角度输入与显示切换逻辑上板验证的核心是能改角度、能看结果。我用两个按键一个控制角度加1度一个控制角度减1度。角度寄存器16位加1度就是加18265536/360。显示用4位数码管前两位显示角度整数部分后两位轮流显示cos和sin的小数部分用一个拨码开关切换显示cos还是sin。这里有个显示技巧cos和sin是Q1.15格式转成十进制显示需要做定点转BCD。我图省事直接用乘以10000再取整的方式显示4位小数虽然占点资源但直观。比如cos_out 32767乘以10000除以32768得到9999显示0.9999。// 定点转显示简化版 wire signed [31:0] cos_scaled cos_out * 10000 / 32768;综合器对除以2的幂会优化成移位但32768不是2的幂的简单形式32768 2^15是2的幂所以除法就是右移15位。乘以10000会消耗一个DSPArtix-7的DSP够用不心疼。4.3 实测结果与误差分析上板跑起来后我测了几个关键点输入角度理论cos实测cos理论sin实测sin0°1.00000.99990.00000.000130°0.86600.86610.50000.499945°0.70710.70700.70710.707290°0.00000.00011.00000.9999180°-1.0000-0.99990.00000.0001误差在±0.0002以内主要来源是16次迭代的截断误差和输出截断。想再提高精度把迭代次数加到18次、内部位宽加到20位误差能压到±0.00005但资源占用会上升约30%。对数码管显示来说16次迭代完全够用。有个现象值得说在90°附近cos的理论值是0实测是0.0001这是正常的。因为CORDIC在接近坐标轴时z累加器很难精确归零最后几次迭代会在±1个LSB之间抖动。这不是bug是算法固有特性接受它就行。5. 踩过的坑与调试经验5.1 符号位处理不当导致全盘错误第一次综合上板数码管显示的结果完全乱套正负号全反。查了半天发现是和用混了。Verilog里对有符号数做右移必须用否则高位补0负数右移后变成大正数整个迭代就崩了。这个坑很隐蔽因为仿真时如果testbench里的激励恰好都是正数根本发现不了。我的建议是所有参与CORDIC迭代的寄存器都声明成signed右移一律用并且在testbench里必须加负角度和跨象限的测试用例。5.2 象限折叠的边界角度处理在90°、180°、270°这些边界上象限判断的高两位会跳变如果折叠逻辑没处理好会出现结果突变。比如角度从89.99°到90.01°高两位从00变成01规约后的角度从16383变成16384-163840理论上cos从0变到0、sin从1变到1是连续的。但如果折叠公式写错比如第二象限用了θ-π而不是π-θ结果就会从1突变成-1。我的调试方法是在testbench里扫一遍0到360度每0.1度打一个点把cos和sin的曲线画出来看有没有跳变。这个方法比单点测试有效得多能一次性暴露所有边界问题。5.3 流水线valid信号的对齐CORDIC核有16级流水线意味着输入角度后要等16个时钟才有输出。如果上层模块用组合逻辑直接读cos_out会读到中间态的垃圾值。我一开始就是按键一按就读结果数码管闪得没法看。后来加了一个16级的移位寄存器跟踪valid只有valid拉高时才更新显示寄存器问题解决。这个经验适用于所有流水线设计延迟必须显式跟踪不能假设输出随时有效。很多新手在写FIR、FFT、CORDIC这类流水线模块时都会栽在这上面。5.4 资源占用与优化方向综合报告显示这个CORDIC核在Artix-7上用了约320个LUT、280个FF、1个DSP。LUT主要消耗在16级迭代的加减法和角度查找表上DSP用在显示转换的乘法。如果想省资源可以把迭代次数降到12次精度损失到小数点后3位LUT能降到200左右。如果想提频率可以在每级迭代之间插入寄存器做成全流水但延迟会增加。一个容易被忽略的优化点角度查找表不要用ROM直接用case语句或者常量数组综合器会把它优化成组合逻辑比Block RAM更省资源而且没有读延迟。16项的表用LUT实现绰绰有余。6. 从sin/cos延伸到其他CORDIC应用把旋转模式跑通之后CORDIC的向量模式就很容易理解了。向量模式是把旋转模式的条件反过来不是根据z的符号决定方向而是根据y的符号决定方向目标是让y归零。迭代结束后z累加器里就是atan(y/x)x累加器里是模长乘以增益。用这个可以算arctan、求复数模长、做坐标变换。再进一步CORDIC还能做双曲函数、指数、对数、平方根只要把旋转角度序列换成atanh(2^(-i))迭代公式稍作调整即可。这些在通信里的载波恢复、锁相环、矩阵求逆等场景都有应用。我个人的体会是CORDIC的价值不在于它算得有多快而在于它用极少的硬件资源实现了超越函数的计算这在资源受限的FPGA和ASIC里是无可替代的。如果你已经跑通了sin/cos下一步建议试试用CORDIC做IQ信号的幅度和相位提取这是数字通信里非常实用的一个功能。把I路和Q路分别接到x和y输入向量模式跑一遍模长就是幅度z就是相位一个核同时出两个结果效率很高。
返回列表