ARTICLE DETAIL

资讯详情

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

手写16点基-2 FFT的Verilog实现:从蝶形网络到RTL全流程

手写16点基-2 FFT的Verilog实现:从蝶形网络到RTL全流程 写Verilog不是会写个计数器、状态机就算入门。想把FFT这种在FPGA里真正高频出现的DSP算法跑通需要在算法理解和RTL实现之间反复横跳。16点的基-2 FFT就是极好的切入点——点数不大但蝶形网络、旋转因子、位反转、定点化、乒乓存储、状态机控制一样不少把这一套结构吃透后面做OFDM解调、雷达脉冲压缩、音频频谱分析都能直接迁移思路。这篇文章是我手写16点基-2 FFT的经验记录从算法推导到Verilog代码再到仿真验证和踩坑全流程讲一遍。适合已经有Verilog基础、正在学数字信号处理或者准备做FPGA信号处理的同学参考。如果你只会调IP核想弄明白FFT内部到底发生了什么这篇也值得读完。1. 项目概述与整体设计思路1.1 16点基-2 FFT到底在算什么FFT不是玄学它只是DFT的快速算法。16点DFT的定义式是X(k) Σ x(n) · W_16^(nk)n从0到15k也是0到15其中W_16 e^(-j2π/16)。直接按这个公式算16×16256次复数乘法。基-2 FFT利用了旋转因子的对称性和周期性把计算变成4级蝶形网络每级8个蝶形总共32个蝶形运算每个蝶形只做1次复数乘法计算量直接降到32次。16点规模虽然不算大但省下的这个8倍差距到1024点就会变成几百倍的差距这就是FFT存在的意义。我们这次要做的就是把这4级蝶形网络用Verilog搭出来。输入16个复数样本经过运算后输出16个频域复数结果。在FPGA里复数用实部加虚部表示每个部分一个16位有符号数占32位。整个设计规模不大非常适合手写学习和验证算法正确性。1.2 串行蝶形方案为什么适合入门FPGA上实现FFT常见有三条路直接用厂商IP核、用HLS高级综合、手写RTL。IP核效率最高但内部是黑盒FIFO深度、块浮点、缩放策略这些参数不调过几次根本理解不了HLS写起来像C但生成的逻辑难控制也难调试。手写RTL这条路上又有全并行和串行两种典型架构。全并行方案是把32个蝶形单元全部铺开4级流水线数据像流水线一样涌进去吞吐率很高但代价是消耗大量DSP乘法器16点可能就要上百个一般入门开发板扛不住代码也不好看。我选的是单蝶形复用方案核心就一个蝶形运算单元先把输入数据存进一组寄存器然后依次取出需要运算的A、B两点数据和对应的旋转因子算完写回另一组缓冲区。一个蝶形算完再算下一个4级网络顺序完成。这样做逻辑层次非常清楚方便对照算法一步步看波形、验证中间结果而且资源占用极低普通板子轻松跑。1.3 数据流与存储结构规划整体数据流是这样设计的输入端16个复数样本按顺序进来先做位反转重排写入buffer04级蝶形运算逐级处理每完成一级读写buffer交换一次最终结果放在最后写入的那个buffer里。因为整个过程只有两个buffer轮流读写所以叫乒乓结构后面代码里会详细展开。存储这块我用寄存器组实现而不是RAM原语原因很简单16点规模太小用RAM反而浪费寄存器组读写逻辑直观也便于仿真时直接查看每个点的中间值。真实工程里点数如果到1024以上这里就要换Block RAM了这个我们在最后扩展部分会说。1.4 数据格式与定点化考虑FPGA内部没有浮点单元所以输入数据和旋转因子都要定点化。输入实部虚部各16位有符号数取值区间大致在[-32768, 32767]可以理解为Q1.15格式整数位1位小数位15位。旋转因子的实部虚部绝对值都不超过1直接用16位有符号数表示相当于一个小数点在符号位后的定点数范围约[-1, 1)。这里有个新人比较容易忽略的点FFT每一级蝶形运算都会让数据幅度最多翻倍4级下来最大增长16倍。如果输入接近满幅中间结果必定溢出。解决办法是在每级蝶形输出后做一次右移1位的缩放让数据幅度稳定在合理范围。代价是最终输出整体缩小了16倍对比算法参考结果时记得除回来。具体怎么在代码里处理第三章会给出实现。2. 算法原理与关键参数设计2.1 蝶形运算的数学推导与符号约定基-2 DIT蝶形运算的基本公式是Y_m A W × BY_{mspan} A - W × BA和B是输入的两点复数数据W是相应的旋转因子span是当前级数据跨距。一个蝶形做一次复数乘法和两次复数加减法。这里W的符号方向极其重要我用的是工程上最常见的DFT符号约定W e^(-j2π/N)也就是旋转因子的虚部是负的。如果符号写反输出频谱会左右颠倒这个问题我见过的踩坑率极高。复数乘法展开来是real(W·B) B_re × W_re - B_im × W_imimag(W·B) B_re × W_im B_im × W_re在Verilog里实现就是4个整数乘法器加2个加法器、2个减法器。16位乘16位得到32位结果再通过右移15位截断回16位完成一次定点复乘。2.2 旋转因子的生成与定点化存储16点基-2 FFT一共只需要8个不同的旋转因子从W^0到W^7。它们的实际值如下表所示我把16位定点近似值也一并列出来方便直接抄进ROM。索引三角函数实际值16位定点表示实部虚部01 0j0x7FFF, 0x000010.9239 - 0.3827j0x7641, 0xCF0420.7071 - 0.7071j0x5A82, 0xA57E30.3827 - 0.9239j0x30FC, 0x89BF40 - 1j0x0000, 0x80005-0.3827 - 0.9239j0xCF04, 0x89BF6-0.7071 - 0.7071j0xA57E, 0xA57E7-0.9239 - 0.3827j0x89BF, 0xCF04可以看到0x8000对应-1.00x7FFF对应约1.0这就是Q1.15格式的全部秘密。实际设计里用Verilog的case语句按索引输出对应的实部虚部值就是只读ROM16点规模根本不需要真的调用Block RAM一个组合逻辑就搞定了。2.3 四级蝶形网络的地址规律分析16点基-2 DIT总共4级每一级的跨距span和使用的旋转因子索引都不同。这是整个RTL实现里最重要的一张表级数跨距span蝶形数量旋转因子索引stage 018W^0stage 128W^0, W^4stage 248W^0, W^2, W^4, W^6stage 388W^0, W^1, ..., W^7对应的地址生成规律是当前级共有N/(2×span)组每组有span只蝶形数据地址A为组号×2×span 组内偏移地址B为地址A加span。旋转因子索引为组内偏移左移(3 - 当前级号)位。这个规律在代码里可以直接写成位运算很优雅后面贴代码大家一看就懂。2.4 位反转与输入重排在按时间抽取的DIT算法中输入序列需要按二进制序号反转后的地址重排16点就是4位地址的bit0和bit3交换、bit1和bit2交换。比如输入第2个样本序号0001位反转后是1000即地址8所以要放到buffer的第8个位置。位反转处理我建议放在数据加载阶段一边写入一边重排一劳永逸。否则在每一级运算时都要处理混乱的地址映射代码会变得很难看。加载完成后4级蝶形的地址计算就完全按照2.3节的规律顺序访问非常干净。3. Verilog核心实现与代码详解3.1 顶层模块划分与端口定义先看顶层模块的端口定义。我习惯把数据位宽做成参数这样将来扩到24位或32位精度时只改一个地方。module fft16_top #( parameter DW 16 )( input wire clk, input wire rst_n, input wire start, input wire signed [DW-1:0] din_re, input wire signed [DW-1:0] din_im, input wire din_valid, output reg busy, output reg done, output reg signed [DW-1:0] dout_re, output reg signed [DW-1:0] dout_im, output reg dout_valid );端口不算多。start拉高后开始接收16点输入数据din_valid每一拍的有效数据都要能采到。busy信号告诉外部当前正在运算done信号表示FFT计算完毕此时从dout_re和dout_im读取频点0的结果。我简化了一下完整16点输出其实就存放在内部buffer里按地址0到15依次读取即可。如果想做一个连续输出接口加个计数器把16个点流水读出就行逻辑不复杂。3.2 双buffer的乒乓切换双buffer在代码里就是两组16深度的寄存器数组一组实部一组虚部再加一个写buffer选择信号和一个读buffer选择信号。reg signed [DW-1:0] buf_re [0:1][0:15]; reg signed [DW-1:0] buf_im [0:1][0:15]; reg rd_buf; // 当前读buffer reg wr_buf; // 当前写buffer加载阶段写wr_buf为0的buffer运算开始后从rd_buf为0的buffer读计算完的结果写入wr_buf为1的buffer。这一级8个蝶形全部算完将rd_buf和wr_buf都取反下一级从Buffer1读、往Buffer0写。这样读写物理上分离不会出现同一时刻读和写同一块存储的冲突问题。这个乒乓结构在流式数据处理里非常常见很多同学学到这里只记得“双buffer能提高速度”其实它更本质的作用是彻底解决读写冲突让每一级运算可以在确定的时间窗口内完成。3.3 组合蝶形单元与复数乘法蝶形单元我直接写成组合逻辑先取数再复乘再做加减和右移最后在时钟沿写回。单周期完成一只蝶形对于16点规模的教学设计完全够用。在高时钟频率场景下这段组合路径可能需要插入流水寄存器后面第4章会给出优化提示。// 数据读取 wire signed [DW-1:0] a_re buf_re[rd_buf][addr_a]; wire signed [DW-1:0] a_im buf_im[rd_buf][addr_a]; wire signed [DW-1:0] b_re buf_re[rd_buf][addr_b]; wire signed [DW-1:0] b_im buf_im[rd_buf][addr_b]; // 复乘 W*B wire signed [2*DW-1:0] b_w_re b_re * w_re - b_im * w_im; wire signed [2*DW-1:0] b_w_im b_re * w_im b_im * w_re; // 截断回16位, 先偶数右移15位 wire signed [DW-1:0] b_w_re_t b_w_re[2*DW-2:DW-1]; wire signed [DW-1:0] b_w_im_t b_w_im[2*DW-2:DW-1]; // 蝶形输出并右移1位做防溢出缩放 wire signed [DW-1:0] y0_re (a_re b_w_re_t) 1; wire signed [DW-1:0] y0_im (a_im b_w_im_t) 1; wire signed [DW-1:0] y1_re (a_re - b_w_re_t) 1; wire signed [DW-1:0] y1_im (a_im - b_w_im_t) 1;这里两个细节值得展开。第一32位乘法结果取[30:15]这段等价于算术右移15位把Q1.15格式的旋转因子乘法结果缩回到Q1.15。第二蝶形加减后我用无符号右移运算符做算术右移1位这是有符号数的正确右移方式能保持符号位。很多初学者在这用最终负数会变成很大的正数频谱直接是错的。3.4 旋转因子表用case实现ROM旋转因子表用组合逻辑case实现这里需要特别小心的是虚部符号我在注释里标清楚防止以后自己看代码时被绕晕。reg signed [DW-1:0] w_re; reg signed [DW-1:0] w_im; always (*) begin case (w_addr) 3d0: begin w_re 16sh7FFF; w_im 16sh0000; end // 1.0000 0.0000j 3d1: begin w_re 16sh7641; w_im 16shCF04; end // 0.9239 - 0.3827j 3d2: begin w_re 16sh5A82; w_im 16shA57E; end // 0.7071 - 0.7071j 3d3: begin w_re 16sh30FC; w_im 16sh89BF; end // 0.3827 - 0.9239j 3d4: begin w_re 16sh0000; w_im 16sh8000; end // 0.0000 - 1.0000j 3d5: begin w_re 16shCF04; w_im 16sh89BF; end // -0.3827 - 0.9239j 3d6: begin w_re 16shA57E; w_im 16shA57E; end // -0.7071 - 0.7071j 3d7: begin w_re 16sh89BF; w_im 16shCF04; end // -0.9239 - 0.3827j default: begin w_re 16sh7FFF; w_im 16sh0000; end endcase end注意负数在Verilog里用s符号标记16shCF04表示有符号数-0.3827的补码表示直接参与有符号乘法时解释正确。如果忘了加s标记后面对它做乘法时会被当作正数处理结果就全乱了。3.5 控制状态机的时序设计状态机只有4个状态IDLE、LOAD、PROC、DONE。加载16点数据需要16拍运算32只蝶形需要32拍总共不到60拍就能完成一轮FFT。localparam IDLE 3d0; localparam LOAD 3d1; localparam PROC 3d2; localparam DONE 3d3; reg [2:0] state; reg [1:0] stage_cnt; // 0~3 reg [3:0] bf_cnt; // 当前蝶形编号0~7 reg [3:0] load_cnt; // 加载计数LOAD阶段的工作是采样din_valid有效的数据写入当前写buffer的位置地址做位反转。加载完16点后自动跳到PROCstage_cnt清零bf_cnt清零。PROC阶段每次处理一个蝶形计算地址、取旋转因子、执行蝶形运算、写回目标buffer。bf_cnt到7表示本级8只蝶形全部完成此时切换buffer、stage_cnt加1如果已经是第3级则进入DONE状态。这个状态机的核心设计思路是“浅状态深计数”状态本身只要4个真正的复杂度都放在stage_cnt和bf_cnt这两个计数器的配合上调试起来非常直观用modelsim看波形时一目了然。3.6 地址与旋转因子索引生成地址生成和旋转因子索引是整个代码里最精华的部分用位运算直接映射2.3节的规律。wire [3:0] span 4d1 stage_cnt; // 1, 2, 4, 8 wire [3:0] group_id bf_cnt stage_cnt; // 组号 wire [3:0] offset bf_cnt (span - 4d1); // 组内偏移 wire [3:0] addr_a (group_id (stage_cnt 1b1)) | offset; wire [3:0] addr_b addr_a span; wire [2:0] w_addr offset (3d3 - stage_cnt);拿stage 2验证一下这时stage_cnt等于2span4bf_cnt从0循环到7。bf_cnt为5时group_id521offset531所以addr_a(13)|19addr_b9413正好是第二组里第2号蝶形取地址9和13正确。旋转因子索引w_addr1(3-2)2查表得W^2和2.3节舞台参数表一致。这段组合逻辑我在仿真里反复对照过公式确认在4级所有情况下都成立。扩展到大点数时只需要把stage_cnt的位宽从2位改成log2(N)对应的位数公式完全不用动这也是这个写法比较高级的地方。4. 仿真验证、问题排查与扩展建议4.1 测试平台与参考数据对比仿真我推荐用Icarus Verilog加GTKWave轻量免费适合学习调试。编译运行命令如下iverilog -o fft_sim.vvp fft16_top.v fft16_tb.v vvp fft_sim.vvp gtkwave fft_sim.vcd测试向量我用一个最简单也最能说明问题的信号单频余弦x[n] cos(2π·2n/16)n0到15。这个信号理论上在频点2和频点14处各有一条谱线幅度为N/28。用Python的numpy.fft.fft先算一遍参考结果然后和FPGA仿真输出做比对。实际验证时输入数据要定点化把cos值乘以32767取整作为16位有符号数喂给模块。FPGA输出由于每级右移1位整体是参考结果的1/16所以对比时把Python结果除以16再比较。我实测下来频点2和14的幅度和参考值误差在1%以内虚部是一个很小的残差这就是定点化引入的量化噪声。4.2 精度评估与误差控制用16位定点做FFT误差来源主要有三个。第一个是输入数据的量化误差正弦值取整到整数本身就有±0.5的量化噪声。第二个是旋转因子存储误差我用的16位近似值和真实三角函数值有约零点几个百分点的偏差。第三个是每级右移1位时的舍入误差直接截断会引入直流偏置更精细的做法是加0.5后取整也就是四舍五入。我实测这个实现对比浮点FFT信噪比大约能做到60到70dB对绝大多数信号处理场景已经够用。如果你想要更高精度把数据位宽从16位升到24位或者用块浮点策略SNR会迅速改善但如果只是想学习FFT体系16位这个精度水平已经能完美暴露所有算法问题是最合适的教学配置。4.3 常见问题排查实录我在调试过程中遇到的坑这里整理成速查表希望你们能少走弯路。现象可能原因解决办法输出频谱左右颠倒旋转因子虚部符号写反确认W e^(-j2π/N)虚部为负结果全是很大的正数有符号数用了逻辑右移用算术右移点数错乱、地址越界加载时没做位反转或位反转写错检查4位地址bit0/bit3、bit1/bit2交换中间结果溢出变成负大数级间没做缩放每级蝶形输出右移1位状态机卡死在busy级结束标志bf_cnt判断时机不对检查bf_cnt在8只蝶形后的清零逻辑结果整体偏小级间缩放太多确认每级只右移1位总缩放为N其中“输出频谱左右颠倒”这个坑我印象最深。第一次跑出结果时频点2的能量跑到了频点14上折腾半天以为是地址错乱最后发现是旋转因子的虚部符号写反了。后来我把符号约定直接写进代码注释再没犯过这个错。4.4 性能评估与资源优化建议以Cyclone IV这类入门级FPGA为例单蝶形复用方案占用大约1000个左右的寄存器主要就是两组16×16×2的寄存器数组DSP乘法器只用4个因为复乘只需要4个整数乘法器而全并行方案需要128个。时序方面16位乘16位加两个加法器的组合路径在50MHz下完全稳定100MHz以上建议在蝶形单元中间插入流水寄存器。如果追求吞吐率建议改成流水线架构每一级配一个独立的蝶形运算单元级和级之间用寄存器打拍。这样数据可以连续不断地流进第一级完成第一个点的时间是4拍左右之后每个时钟周期出一个频点吞吐率是单蝶形方案的8倍。代价是面积增加约4倍这就是典型的面积换速度。4.5 从16点到更大点数的扩展路径把16点改成64点、256点需要改动的部分其实非常少。第一buffer深度从16改到N第二stage_cnt位宽从2位改成log2(N)对应的位数第三旋转因子表从8组扩展到N/2组如果N很大就改用CORDIC算法实时计算。地址生成公式和状态机结构完全不用动这也是当初选这个高度规整的基-2结构的原因。我自己做完16点之后顺手扩到了64点改代码花了不到半小时仿真一次通过。那种感觉就是之前花在理解地址规律、调试位反转上的时间全部回本了。最后再分享一个小习惯调试FFT时我最常用的方法是把每一级处理完的中间buffer用$display打印出来和Python里按级手算的结果逐级对比。这样一旦出错能立刻定位到是哪一级的哪只蝶形出了问题而不是在最终输出面前瞎猜。你如果刚写完这个模块不妨也试试这个手段会让你对整个蝶形网络的理解深一大截。
返回列表