ARTICLE DETAIL

资讯详情

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

FPGA中CORDIC+LUT实现sin/cos:从迭代原理到Verilog状态机

FPGA中CORDIC+LUT实现sin/cos:从迭代原理到Verilog状态机 简介面向FPGA初、中级学习者与高校教研人群资源包内含基于LUT查找表方法的Cordic算法Verilog实现可输出正弦和余弦波形并配套完整testbench测试程序适用于课程设计、算法验证或项目参考。资源包为rar压缩格式共496个文件约31.03MB包含Vivado工程文件xpr、v、vcd、仿真与编译脚本bat、tcl、综合与仿真日志log、rpt、wdb及演示操作视频avi目录结构清晰便于按需检索。目前已有1776人学习下载经作者实测需使用Vivado2019.2及以上版本打开工程后参考录像操作即可快速复现。资源亮点在于形成从代码、测试激励到操作演示的完整闭环尤其借助testbench和自动化脚本能直观理解Cordic算法在不同角度下正余弦计算过程同时工程路径必须为英文这一注意事项也帮助读者避开常见的中文路径报错问题。1. 先想清楚CORDIC配LUT要解决什么问题把相位转成正弦、余弦FPGA 里最常见的两条路是查表和迭代。直接查表把整个周期的波形存成 ROM相位宽度每加 1 位表容量就翻一倍到 20 位相位时单张表已经 2MBI/Q 两路就是 4MB中小型器件很难接受。CORDIC 换了个思路不存波形只存一串“旋转角常数”每次迭代靠加法和移位逼近目标角度资源消耗几乎只跟迭代次数有关。13 次迭代能够把 sin/cos 做到千分之一量级这串旋转角常数本身就是一张很小的 LUT。把 LUT 当系数存储器用是数字下变频、电机矢量控制、相噪分析里最稳的做法。这篇按“迭代公式 → 定点化 → Verilog 状态机 → testbench 验证 → 上板调试”的顺序把整条链路讲清楚代码可以直接建工程跑仿真。2. CORDIC旋转迭代与LUT表项怎么对齐2.1 旋转模式的基本递推式CORDIC 的核心思想是把一个向量不断旋转每次旋转的角度越来越小最终让向量转到目标角度上。旋转模式Rotation Mode要求把角度寄存器 z 逐步压到 0递推式如下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 由 z 的符号决定z ≥ 0 时取 1z 0 时取 -1。这里的 arctan(2^{-i}) 就是每次旋转的固定步长i0 时是 45°i1 时约 26.565°越往后步长越小。硬件上 2^{-i} 的乘法就是右移 i 位这是 CORDIC 能在 FPGA 上低成本实现的关键。所谓 LUT 查找表方法就是指把这一串 arctan 值预先算好存进一张常数表迭代到第 i 轮时直接查表取角度不需要现场用除法或者泰勒展开去算。旋转的向量公式展开后第 i 轮等价于对 (x, y) 做一次 Givens 旋转旋转角为 d_i * arctan(2^{-i})。连续做完 N 轮总旋转角就是各轮旋转角的代数和当 z 收敛到 0 时这个代数和恰好等于输入角度。2.2 增益补偿与初值选取每次旋转除了转动向量还会把模长放大 sqrt(1 2^{-2i})。N 轮迭代的累计增益为A ∏ sqrt(1 2^{-2i}) ≈ 1.64676 (N → ∞)如果直接把 (1, 0) 作为初始向量迭代结束后得到的是放大过的 A 倍 sin/cos不是真实值。解决办法是把 x 初始化为 1/A 的定点值让后续放大刚好抵消x_0 1 / A ≈ 0.6073 y_0 0 z_0 目标角度迭代完成后 x_N 就是 cos(z0)y_N 就是 sin(z0)。1/A 这个常数在 16 位 Q1.14 定点格式下为 0.6073 * 16384 ≈ 9950也就是 16 进制 0x26DE。注意这个补偿常数在迭代次数变化时差不太多N ≥ 10 后 A 基本收敛到 1.64676直接用 9950 不会有可感知误差。2.3 象限归约与 atan 表生成CORDIC 的收敛域大约在 [-99.88°, 99.88°]超过这个范围直接输入角度会不收敛。处理办法是先做象限归约把任意角度折到 [-90°, 90°] 内。FPGA 里做一个 16 位相位字让 2π 对应 16 位无符号整数 65536取角度的高 2 位判断象限剩下的 14 位送到迭代器。下面是 13 轮迭代需要的 atan 表角度值按 65536 对应 2π 换算例如 45° 对应 65536 * 45 / 360 8192。我给出的数值是四舍五入后的定点整数索引 i旋转角度16位定点值说明045.0004096arctan(1)126.5652419arctan(0.5)214.0361278arctan(0.25)37.125649arctan(0.125)43.576326arctan(0.0625)51.790163arctan(1/32)60.89582arctan(1/64)70.44841arctan(1/128)80.22420arctan(1/256)90.11210arctan(1/512)100.0565arctan(1/1024)110.0283arctan(1/2048)120.0141arctan(1/4096)生成规则很简单角度值 atan(2^{-i}) / (2 * π) * 65536取整即可。i11 和 i12 已经小于 3 个 LSB再往后即使继续迭代对 16 位输出也几乎没有贡献。这张表放进 Verilog 的 function 里综合时会映射成常数多路选择器也就是纯组合逻辑的 LUT不会占 BRAM。3. Verilog实现定点格式、数据通路与状态机3.1 定点格式与顶层端口x、y 方向分量用有符号 16 位 Q1.14 定点数1.0 对应 16384最小分辨率为 2^{-14}约 6.1e-5。角度输入用无符号 16 位2π 对应 65536这样相位累加器可以直接给到 angle_in不需要额外转换。输出同样是有符号 16 位 Q1.14sin_out/cos_out 除以 16384 就是真实的浮点数值。顶层端口定义如下端口方向位宽说明clkinput1系统时钟rst_ninput1异步复位低有效startinput1输入有效标志高电平启动一轮计算angle_ininput16无符号相位2π 65536busyoutput1计算中拉高禁止再次 startdoneoutput1计算结果有效持续一个时钟周期sin_outoutput16正弦输出Q1.14 有符号cos_outoutput16余弦输出Q1.14 有符号迭代次数通过参数 ITER_NUM 控制默认 13。状态机只有 IDLE、RUN、DONE 三个状态每次 start 后连续迭代 ITER_NUM 轮结束后在 DONE 状态输出象限校正结果。用状态机而不是流水线好处是代码直观、资源极少适合先把算法跑通吞吐量要求高的场景再改成流水线结构。3.2 用 function 生成 atan 查找表atan 常数表直接用 Verilog function 实现综合工具会把它优化成常数 LUT。查找索引是迭代计数 iter_cnt每个时钟周期查一次延迟为零符合状态机逐轮迭代的节奏。function signed [15:0] atan_lut; input integer idx; begin case (idx) 0: atan_lut 16sd4096; 1: atan_lut 16sd2419; 2: atan_lut 16sd1278; 3: atan_lut 16sd649; 4: atan_lut 16sd326; 5: atan_lut 16sd163; 6: atan_lut 16sd82; 7: atan_lut 16sd41; 8: atan_lut 16sd20; 9: atan_lut 16sd10; 10: atan_lut 16sd5; 11: atan_lut 16sd3; 12: atan_lut 16sd1; default: atan_lut 16sd0; endcase end endfunction这里所有常量都写成 16sd 形式明确告诉综合器是有符号数。如果写成 16d4096在后续与 signed 寄存器做加减法时工具会按无符号处理容易产生位宽扩展警告严重时迭代结果完全错误。这个细节对可综合性很重要。default 分支保底防止综合时生成锁存器。3.3 迭代状态机与象限校正主模块采用三段式框架组合逻辑算下一拍数据时序逻辑做状态跳转。核心迭代在 RUN 状态完成x/y/z 三个寄存器的更新逻辑是这套实现的心脏。module cordic_sincos_lut #( parameter DATA_WIDTH 16, parameter ITER_NUM 13 )( input wire clk, input wire rst_n, input wire start, input wire [DATA_WIDTH-1:0] angle_in, output reg busy, output reg done, output reg signed [DATA_WIDTH-1:0] sin_out, output reg signed [DATA_WIDTH-1:0] cos_out ); localparam signed [DATA_WIDTH-1:0] K_CONST 16sd9950; localparam IDLE 2d0; localparam RUN 2d1; localparam DONE 2d2; reg [1:0] state; reg [3:0] iter_cnt; reg [1:0] quad_reg; reg signed [DATA_WIDTH-1:0] x_reg, y_reg, z_reg; reg signed [DATA_WIDTH-1:0] x_next, y_next, z_next; wire signed [DATA_WIDTH-1:0] x_shift x_reg iter_cnt; wire signed [DATA_WIDTH-1:0] y_shift y_reg iter_cnt; // 象限判定与z初值组合逻辑 reg [1:0] quad_sel; reg signed [DATA_WIDTH-1:0] z_src; always (*) begin case (angle_in[DATA_WIDTH-1:DATA_WIDTH-2]) 2b00: begin quad_sel 2b00; z_src $signed(angle_in); end 2b01: begin quad_sel 2b01; z_src 16sd32768 - $signed(angle_in); end 2b10: begin quad_sel 2b10; z_src $signed(angle_in) - 16sd32768; end 2b11: begin quad_sel 2b11; z_src -$signed(angle_in); end endcase end // 迭代组合逻辑根据z符号决定旋转方向 always (*) begin if (z_reg[15] 1b0) begin x_next x_reg - y_shift; y_next y_reg x_shift; z_next z_reg - atan_lut(iter_cnt); end else begin x_next x_reg y_shift; y_next y_reg - x_shift; z_next z_reg atan_lut(iter_cnt); end end always (posedge clk or negedge rst_n) begin if (!rst_n) begin state IDLE; busy 1b0; done 1b0; x_reg 16sd0; y_reg 16sd0; z_reg 16sd0; iter_cnt 4d0; quad_reg 2d0; sin_out 16sd0; cos_out 16sd0; end else begin case (state) IDLE: begin done 1b0; if (start) begin busy 1b1; x_reg K_CONST; y_reg 16sd0; z_reg z_src; quad_reg quad_sel; iter_cnt 4d0; state RUN; end end RUN: begin x_reg x_next; y_reg y_next; z_reg z_next; if (iter_cnt ITER_NUM - 1) begin state DONE; end else begin iter_cnt iter_cnt 1b1; end end DONE: begin case (quad_reg) 2b00: begin sin_out y_reg; cos_out x_reg; end 2b01: begin sin_out y_reg; cos_out -x_reg; end 2b10: begin sin_out -y_reg; cos_out -x_reg; end 2b11: begin sin_out -y_reg; cos_out x_reg; end endcase busy 1b0; done 1b1; state IDLE; end endcase end end endmodule代码里几个关键点需要注意。x_shift 和 y_shift 用连续赋值 wire由 x_reg 和 iter_cnt 实时决定组合逻辑块根据 z_reg 的符号位选择旋转方向这决定了 x/y/z_next 的加减方式。z_reg[15] 为 1 表示负数此时执行反向旋转把 z 往 0 方向拉。RUN 状态里每次时钟上升沿把 x_next/y_next/z_next 锁存iter_cnt 加 1最后一次迭代完成后直接跳 DONE下一拍 x_reg 就是最终迭代结果DONE 中读取 x_reg/y_reg 并做象限校正。象限校正的 4 种情况对应最初折叠时丢失的符号信息quad_reg原始角度范围校正公式000° ~ 90°sin y, cos x0190° ~ 180°sin y, cos -x10180° ~ 270°sin -y, cos -x11270° ~ 360°sin -y, cos x例如 135° 折叠后 z_src 32768 - 24576 8192对应 45°迭代得到 sin45/cos45再按 quad_reg01 校正sin 保持正、cos 取负正好是 sin135° 和 cos135°。4. Testbench设计与仿真结果分析4.1 激励组织定点角度与期望值testbench 要覆盖四个象限的关键角度同时包含一个非特殊角验证普适性。角度值按 LSB 360° / 65536 计算例如 30° 就是 65536 * 30 / 360 ≈ 5461。期望值直接用浮点常数写进数组不依赖仿真器的数学库这样 Modelsim 和 Vivado xsim 都能原样运行。测试点选择0°、30°、45°、60°、90° 覆盖第一象限135°、180° 覆盖第二、三象限边界225°、270°、315° 覆盖第四象限最后加一个 125° 的非特殊角。125° 的 sin 参考值是 sin(125°) 0.819152cos 参考值是 -0.573576人为设一个不对称角度比只测特殊角更能暴露象限校正逻辑的问题。4.2 误差计算与检查逻辑timescale 1ns/1ps module tb_cordic_sincos; reg clk 0; reg rst_n 0; reg start 0; reg [15:0] angle_in 0; wire busy; wire done; wire signed [15:0] sin_out; wire signed [15:0] cos_out; cordic_sincos_lut #( .DATA_WIDTH(16), .ITER_NUM(13) ) dut ( .clk (clk), .rst_n (rst_n), .start (start), .angle_in(angle_in), .busy (busy), .done (done), .sin_out (sin_out), .cos_out (cos_out) ); always #5 clk ~clk; reg [15:0] angles [0:10]; real sin_exp [0:10]; real cos_exp [0:10]; integer i; real sin_got, cos_got; real err_sin, err_cos, max_err_sin, max_err_cos; initial begin angles[0] 16d0; sin_exp[0] 0.0; cos_exp[0] 1.0; angles[1] 16d5461; sin_exp[1] 0.5; cos_exp[1] 0.866025403; angles[2] 16d8192; sin_exp[2] 0.707106781; cos_exp[2] 0.707106781; angles[3] 16d10923; sin_exp[3] 0.866025403; cos_exp[3] 0.5; angles[4] 16d16384; sin_exp[4] 1.0; cos_exp[4] 0.0; angles[5] 16d24576; sin_exp[5] 0.707106781; cos_exp[5] -0.707106781; angles[6] 16d32768; sin_exp[6] 0.0; cos_exp[6] -1.0; angles[7] 16d40960; sin_exp[7] -0.707106781; cos_exp[7] -0.707106781; angles[8] 16d49152; sin_exp[8] -1.0; cos_exp[8] 0.0; angles[9] 16d57344; sin_exp[9] -0.707106781; cos_exp[9] 0.707106781; angles[10] 16d22756; sin_exp[10] 0.819152044; cos_exp[10] -0.573576436; max_err_sin 0.0; max_err_cos 0.0; #100; rst_n 1; for (i 0; i 11; i i 1) begin (posedge clk); angle_in angles[i]; start 1b1; (posedge clk); start 1b0; wait (done 1b1); (posedge clk); sin_got $signed(sin_out) / 16384.0; cos_got $signed(cos_out) / 16384.0; err_sin sin_got - sin_exp[i]; err_cos cos_got - cos_exp[i]; if (err_sin 0.0) err_sin -err_sin; if (err_cos 0.0) err_cos -err_cos; if (err_sin max_err_sin) max_err_sin err_sin; if (err_cos max_err_cos) max_err_cos err_cos; $display(angle%0d sin%0.6f exp%0.6f cos%0.6f exp%0.6f err_sin%0.6f err_cos%0.6f, angles[i], sin_got, sin_exp[i], cos_got, cos_exp[i], err_sin, err_cos); end $display(max_err_sin%0.6f max_err_cos%0.6f, max_err_sin, max_err_cos); #100; $finish; end endmoduletestbench 的流程是每个测试点先等时钟上升沿拉高 start 一拍然后 wait(done) 等待计算完成done 拉高后再等一个时钟沿确保输出寄存器稳定。$signed(sin_out) / 16384.0 把定点输出还原成浮点。误差累加器记录整组测试的最大误差最后打印。wait(done) 在信号已经为高时会立即通过不会造成死锁但如果 DONE 状态持续多拍wait 可能连续通过所以 testbench 里让 done 只在 DONE 状态输出一个周期是很必要的。另外要注意 start 只能拉高一拍。如果在 busy1 期间再次 start状态机还停在 RUN并不会重新加载初值但可能造成迭代计数错乱。工程上可以在 IDLE 外再把 start 打一拍做握手或者干脆用 valid/ready 信号替代简单 start。4.3 仿真输出与误差量化跑完 11 个测试点后$display 会打印每组角度对应的 sin/cos 实测值和期望值。以 30° 为例angle_in5461 时sin_out 应该在 8192 附近即 0.5 的定点表示cos_out 在 14189 附近0.866025 * 16384 ≈ 14188.7。125° 那组能验证象限校正的符号逻辑sin 为正、cos 为负。13 次迭代 16 位数据位宽的典型误差来源有三个角度量化误差约 0.005°迭代残差约 2^{-12}以及定点舍入累积。综合下来最大误差在 1e-3 量级多数角度在 5e-4 以内。看打印结果时优先关注最大误差的测试点通常是 30° 或 60° 这种非整数相位字的角度因为角度量化本身就是误差的一部分。如果 max_err 超过 0.005先查 z_src 的折叠逻辑再查 x/y 是否在迭代中出现符号位扩展错误。5. 精度、资源占用与三个常见坑5.1 迭代次数与位宽匹配迭代次数不是越多越好。第 i 轮旋转角是 arctan(2^{-i})到第 13 轮已经只有 0.014°16 位定点数据的 LSB 是 0.0055°180° / 32768两者已经相当接近。继续增加迭代旋转角小于数据分辨率角度寄存器 z 根本感知不到变化白白增加时延。下面是不同迭代次数对精度的粗略影响ITER_NUM最小旋转角度16位定点下可感知建议场景80.224可以粗粒度相位误差约0.3°100.056勉强电机控制误差约0.08°120.028边缘通信解调误差约0.04°130.014低于LSB分辨率推荐误差小于0.01°数据位宽如果从 16 位扩到 20 位迭代次数可以相应增加到 1516 次此时 LSB 更小小角度迭代才有意义。反过来12 位位宽配 13 次迭代就是浪费最后几轮 z 的变化在 x/y 的量化噪声里完全消失。迭代次数和位宽的关系是迭代次数 ≈ 位宽 * 0.8这是快速估算的常用经验值。5.2 有符号移位 与 的区别Verilog 里 是无符号逻辑右移 是算术右移。x_reg 和 y_reg 是有符号 reg如果误写为 负数右移时高位补 0符号就丢了。比如 y_reg -163840xC000 1 得到 0x600024576而 1 得到 0xE000-8192。CORDIC 迭代中 x/y 会频繁出现负值位移结果直接影响加减方向一个符号错误会导致输出完全错乱且很难排查。如果工具链对 支持不佳可以手动拼接符号位扩展wire signed [15:0] y_shift { {16{y_reg[15]}}, y_reg } (16 iter_cnt);这段代码先把 y_reg 扩展到 32 位并复制符号位再整体右移取低 16 位就是算术右移效果。实际使用中 Vivado、ModelSim、Quartus 都支持 只有很老的 95 风格代码才需要这种替换。排查仿真波形时如果 x/y 出现非对称的异常跳变优先检查右移是否变成了逻辑移位。5.3 象限边界与角度格式陷阱角度输入是无符号 16 位2π 对应 65536也就是 0 和 65536 是同一个相位。但 65536 在 16 位里存不下实际输入只能是 0 到 65535。所以 360° 等价于 0° 输入而 65535 表示 359.99°接近 0° 而非 360°。这在象限判断上没有影响因为 65535 的高 2 位是 11走第 4 象限折叠得到 z_src 1约 0.0055°输出 cos 接近 1、sin 接近 0结果是正确的。真正的坑在 180° 边界。angle_in 0x8000 32768高 2 位是 10z_src 32768 - 32768 0迭代输出 cos1、sin0象限校正后 cos-1、sin0正确。但如果在 90° 边界angle_in 0x4000 16384高 2 位是 01z_src 32768 - 16384 16384也就是 90° 折叠后依然是 90°这超出了 CORDIC 的收敛域上限。所以必须注意折叠后的 z_src 应该落在 [0, 16384) 区间90° 边界在 Q1 和 Q2 的共用点上Q1 折叠保持 16384 不变Q2 折叠得到 16384两个结果都指向 90°而 CORDIC 在 90° 处收敛正常实际误差小于 0.02°。如果迭代次数较少这个边界误差会偏大建议 ITER_NUM 不要低于 10。6. 上板后的验证与调试技巧6.1 用 Modelsim 命令行做自动化回归写完代码后不要在 GUI 里手动点按钮加信号效率太低。用命令行方式把仿真流程固化下来改一次参数重跑一遍只需要几秒。在工程目录下创建 run.do 文件vlib work vlog -sv cordic_sincos_lut.v tb_cordic_sincos.v vsim -c -do run -all; quit work.tb_cordic_sincosvsim 的 -c 参数表示 console 模式不弹图形界面。run -all 跑完整个测试$display 打印的结果直接输出到终端。如果想把结果存成文件在 testbench 里用 $fopen/$fdisplay 写日志或者命令行加-l sim.log把全部输出重定向到日志文件。这种方式适合反复调整 ITER_NUM 和 DATA_WIDTH 时对比误差变化。检查波形时重点看三组信号z_reg 是否单调逼近 0x_reg 是否始终在正负 16384 之间iter_cnt 是否按顺序递增。z_reg 如果在某轮迭代后远离 0说明 atan_lut 的索引或旋转方向反了。把这三组信号加到波形窗口一眼就能定位是折叠逻辑还是迭代逻辑出了问题。6.2 用 ILA 观察迭代内部状态纯仿真通过后上板前最好在 Vivado 里加 ILA 核观察实际运行时的迭代过程。ILA 和 testbench 看到的都是 x_reg、y_reg、z_reg、iter_cnt、state 这些信号区别是 ILA 抓的是真实硬件时序能暴露时钟域、复位释放时间这类仿真里发现不了的问题。在 Vivado 中把 cordic_sincos_lut 的 x_reg、y_reg、z_reg、iter_cnt、quad_reg、state 加入 debug 列表采样深度选 1024触发条件设为 done 上升沿。实际抓波形时会看到每轮迭代 x/y/z 的变化轨迹。如果 done 触发后发现 z_reg 在迭代前几轮就归零说明 atan_lut 的数值给大了z 被过冲后再拉回来如果 z_reg 一直不归零说明角度折叠后的 z_src 超出了收敛域检查角度输入的高 2 位映射是否和预期一致。上板调试还有一个讨巧的做法把 sin_out 直接接到 Pmod 引脚上外接示波器看波形但这样只能看个大概正弦的毛刺和符号翻转不好分辨。更实用的做法是先把 testbench 里 10 个测试点的误差全部跑通再用 ILA 抽查一个 30° 和一个 125° 的内部数据这样既能验证算法正确性又不会因为反复上板浪费时间。如果后续要做连续相位输入比如 DDS 累加器驱动把 start 改成每 N 个时钟周期自动拉高一次即可。Xilinx 也有现成的 CORDIC IP 核可以直接用参数化界面里选 sin/cos 模式和流水线结构它在内部做的就是这篇文章里同样的迭代过程只是把 x/y 数据通路做成了全流水。自己写 Verilog 版本的价值在于能精确控制角度格式、迭代次数和输出时延在非标准数据接口的项目里反而比 IP 核更灵活。本文还有配套的精品资源点击获取
返回列表