ARTICLE DETAIL

资讯详情

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

手写Verilog基-2 16点FFT:从蝶形运算到FPGA实现全解析

手写Verilog基-2 16点FFT:从蝶形运算到FPGA实现全解析 一提到FFT很多FPGA工程师的习惯性动作是打开Vivado或Quartus里的FFT IP核把点数一配、流模式一选、点一下Generate完事。但今年我在一个资源预算比较紧的小项目里只需要做16点FFT输入是8bit定点数外部数据速率不高而核的Latency、DSP占用和时序行为总觉得有点“黑盒”。折腾了一圈之后我决定干脆用Verilog手写一个基-2 16点FFT。这个规模不算大但正好能把FFT硬件化的全部关键问题——蝶形运算结构、旋转因子量化、定点数位宽控制、数据调度、仿真对比——完整地走一遍。做完之后回头看这件事的价值远不止“省了一个IP核”它让我对FFT的数据流和FPGA数字信号处理的设计思路都有了更实的掌握。这篇文章就按我实际的实现路径把算法结构、RTL拆分、定位数处理、验证方法和踩过的坑一次讲清楚。1. 先算清楚这笔账16点FFT为什么值得手写1.1 资源占用与延迟的直观对比基-2 16点FFT的运算量是固定的4级蝶形每级8个蝶形总共32次复数蝶形运算。但同样的算法映射到FPGA上的架构差别非常大。全并行架构每一级8个蝶形每个蝶形一个复数乘法器总共需要32个复数乘法器等价于128个实数乘法器逻辑规模和布线压力都不小。单蝶形复用架构只实现一个蝶形运算单元用状态机控制它依次完成32次运算每次运算只用一个复数乘法器资源占用极小。IP核以Xilinx FFT IP为例配16点、Radix-2 Burst I/O模式Latency通常在百拍量级资源则取决于你选的实现方式。IP核的好处是通用性强但代价是“黑盒”综合之后未必是最优解。我把三种方案在16点这个规模下的表现粗略列了一下方案复数乘法器数量大致Latency可控性适用场景全并行/流水线32数拍到十几拍高高速流式处理单蝶形复用1约100拍很高低速、资源受限FFT IP核取决于配置数十到数百拍低通用快速集成1.2 什么情况下应该放弃IP核用手写IP核不是不能用而是有些场景下它确实不是最优选择。资源受限的小规模设计只做16点FFT用IP核有点“杀鸡用牛刀”配置页翻半天不说生成的逻辑可能比手写大不少。时序与流水需要精确对齐IP核Latency是固定的但它内部到底怎么调度、能不能和你的上下游逻辑无缝衔接很多时候要靠经验去猜。手写的话每一拍干什么都是自己控制的。需要自定义数据格式与截位策略FFT IP核的输出位宽、缩放方式是可配的但如果你有特殊的定点数格式要求或者想在中间级做自定义处理手写反而更灵活。学习与面试需求现在不少数字IC和FPGA相关的笔试面试题里都有“手撕FFT”这种题自己完整写过一遍理解深度完全不同。16点FFT的手写逻辑量并不大核心模块也就几个蝶形运算单元、旋转因子ROM、地址生成器、状态机。在资源利用率上单蝶形方案用一个DSP48就够逻辑消耗只有几百个LUT。2. 基-2 16点FFT的蝶形结构与旋转因子表2.1 DIT还是DIF选哪一个更顺手FFT的按时间抽取DIT和按频率抽取DIF本质上用同一个蝶形公式区别在于输入输出顺序。DIT要求输入按倒位序排列输出是自然序DIF则是输入自然序输出倒位序。我推荐用DIT原因很直接输入数据写入RAM时做一次bit-reverse地址映射即可后续每级蝶形运算都按自然序读数据控制逻辑简单。DIF虽然输入不需要倒序但最后输出要处理倒序如果下游还要做频谱分析或逆变换反而多一道麻烦。2.2 四级蝶形的排布规律16点基-2 FFT的蝶形间距和旋转因子指数是有严格规律的。设级数从0开始记第s级的蝶形间距是2^s每级固定8个蝶形第s级第j个蝶形的旋转因子指数由公式k j × (16 / (2^(s1))) 计算。实际算一遍就是这个表级数蝶形间距分组数旋转因子指数说明第0级18k0每组1个蝶形第1级24k0,4每组2个蝶形第2级42k0,2,4,6每组4个蝶形第3级81k0,1,...,7一组8个蝶形这里有个容易记错的地方16点FFT第二级旋转因子不是W16^0和W16^2而是W16^0和W16^4。因为第二级里每组两个蝶形合并的只是4点子序列旋转因子对应的是N/4周期上的角度。2.3 输入倒位序的Verilog处理16点需要4位地址反转。比如输入序号1二进制0001反转后地址是1000也就是8。在RTL里写一个反转函数就行function [3:0] bit_reverse4; input [3:0] addr; begin bit_reverse4 {addr[0], addr[1], addr[2], addr[3]}; end endfunction实际写入RAM时就把原始输入数据写到bit_reverse4(地址)的位置。例如第1个采样点写到地址8第2个写到地址4第3个写到地址12依次类推。这个映射可以在Testbench阶段就预先确认一遍常见的一个低级错误就是倒位序只反转了半字节或者反转方向写反导致最终输出频谱顺序是乱的。3. 定点数格式与中间位宽最容易翻车的环节3.1 Q格式定点数的基础FPGA里做不了浮点定点数格式得一开始就定好。最常用的是Q格式表示法。Q1.71位符号位7位小数范围[-1, 1-2^-7]精度约0.0078。适合表示归一化后的信号。Q1.151位符号位15位小数范围[-1, 1-2^-15]精度非常高适合存旋转因子的cos/sin值。在这个设计里输入用Q1.7格式的8bit有符号数即可旋转因子用Q1.15格式的16bit有符号数。两个定点数相乘后结果的格式是Q(1.7) × Q(1.15) ≈ Q(1.22)需要根据需求截位。3.2 两种中间位宽方案逐级右移 vs 全精度增长这块是整个设计里最容易被低估的部分。FFT每级蝶形进行一次加法和一次减法数据位宽理论上每级增长1bit4级下来最多增长4bit。16点FFT的增益正好是16也就是说一个接近满幅的8bit输入经过4级运算后中间结果的理论最大值可以达到输入的16倍。我试验过两种主流方案方案A逐级右移1位每级蝶形计算完成后把结果右移1位相当于每级做一次1/2缩放。4级之后输出等于原始FFT结果的1/16但整个运算过程中的数据位宽恒定不变不会溢出。实现简单资源占用小。代价是每级右移时都会引入量化误差不过16点FFT的级数只有4级误差累积有限实测SNR通常还有50dB以上完全够用。方案B全精度增长最后统一截位输入8bit每级把位宽扩到16bit甚至更宽中间不缩放最后输出时再截位或饱和。精度最高但存储和运算逻辑变多而且最后截位时如果处理不好照样会有溢出风险。两种方案我建议初学者直接用方案A理由很简单省心、无溢出、逻辑清晰。如果你的应用对SNR要求特别高再考虑方案B。3.3 旋转因子的量化与ROM生成旋转因子的计算公式是W16^k cos(2πk/16) − j·sin(2πk/16)在Verilog里存成两个16bit有符号数实部存cos值虚部存−sin值。基-2 16点FFT实际需要用的旋转因子是k0到7另一半可以通过对称性得到所以ROM只需要存8组。用Python生成ROM初始值最方便import math N 16 tw_r [] tw_i [] for k in range(8): angle 2 * math.pi * k / N wr int(round(math.cos(angle) * 32767)) wi -int(round(math.sin(angle) * 32767)) tw_r.append(wr 0xFFFF) tw_i.append(wi 0xFFFF) with open(twiddle_rom.mem, w) as f: for i in range(8): f.write(f{i:02X} {tw_r[i]:04X} {tw_i[i]:04X}\n)需要注意生成时要用16bit有符号数的补码形式存这样才能在Verilog里直接用signed类型读取。4. RTL实现的关键模块蝶形单元、ROM与地址生成4.1 蝶形运算单元的Verilog写法蝶形运算是FFT的核心设输入为x和y旋转因子为W则输出为x y·W和x − y·W。复数乘法展开后需要4个实数乘法和几次加法。写成Verilogmodule butterfly #( parameter DW 16 )( input wire signed [DW-1:0] xr, xi, input wire signed [DW-1:0] yr, yi, input wire signed [DW-1:0] wr, wi, output wire signed [DW-1:0] sum_r, sum_i, output wire signed [DW-1:0] dif_r, dif_i ); // y * W 的实部yr*wr - yi*wi // y * W 的虚部yr*wi yi*wr wire signed [2*DW-1:0] yr_wr yr * wr; wire signed [2*DW-1:0] yi_wi yi * wi; wire signed [2*DW-1:0] yr_wi yr * wi; wire signed [2*DW-1:0] yi_wr yi * wr; wire signed [DW-1:0] ywr (yr_wr - yi_wi) (DW-1); wire signed [DW-1:0] ywi (yr_wi yi_wr) (DW-1); assign sum_r xr ywr; assign sum_i xi ywi; assign dif_r xr - ywr; assign dif_i xi - ywi; endmodule这里有个细节乘法结果是32bit使用算数右移把它截回16bit。右移时研究一下仿真结果看看用截断还是四舍五入对SNR影响大。四舍五入可以在右移前加一个偏移量实现代价是几个加法器。4.2 旋转因子ROM的初始化ROM可以直接用Verilog数组加initial块初始化。我用上面的Python脚本生成twiddle_rom.mem然后在RTL里读入reg signed [15:0] tw_r [0:7]; reg signed [15:0] tw_i [0:7]; initial begin $readmemh(twiddle_rom.mem, tw_r); $readmemh(twiddle_rom.mem, tw_i); end不过$readmemh一次只能读一个数组实践中我通常把实部和虚部分成两个文件或者用initial块里的case语句直接查表。16点规模很小case语句反而最直观可读性也最好。4.3 地址生成与状态机控制单蝶形复用方案需要一个状态机来调度32次蝶形运算。核心控制逻辑是两级计数器stage_cnt[1:0]当前级数0到3。bfly_cnt[2:0]当前级内的蝶形序号0到7。根据这两个计数器和2.2节的规律就可以算出当前蝶形的两个输入数据地址和旋转因子地址wire [3:0] spacing 4b0001 stage_cnt; // 蝶形间距 wire [3:0] j bfly_cnt (spacing - 4b0001); // 组内蝶形序号 wire [3:0] base (bfly_cnt / spacing) * (spacing 1); // 组基地址 wire [3:0] addr0 base j; wire [3:0] addr1 base j spacing; wire [2:0] tw_addr stage_cnt 2d0 ? 3d0 : (stage_cnt 2d1 ? {j[0], 2b00} : (stage_cnt 2d2 ? {j[1:0], 1b0} : j));三级条件嵌套就能把旋转因子换算出来。状态机则负责在每个蝶形运算开始时读两个复数计算完成后写回RAM然后切换下一组地址。5. 仿真验证从激励生成到误差分析的全链路5.1 用Python生成测试向量纯随机数测试能验证功能但看不出精度损失。我习惯用已知频率的正弦信号做激励这样频谱图上一眼就能看出峰值位置对不对。import numpy as np N 16 t np.arange(N) signal np.round(100 * np.cos(2 * np.pi * 3 * t / N)).astype(int) with open(input_real.txt, w) as f: for v in signal: f.write(f{v 0xFF:02X}\n) # 虚部全0只测实信号 with open(input_imag.txt, w) as f: for _ in range(N): f.write(00\n)Testbench里用$readmemh把这两个文件读入存储数组然后启动FFTreg [7:0] mem_r [0:15]; reg [7:0] mem_i [0:15]; reg start; initial begin $readmemh(input_real.txt, mem_r); $readmemh(input_imag.txt, mem_i); #20 start 1; end5.2 结果比对与SNR计算仿真结束后把输出结果导出在Python里与numpy.fft.fft的参考结果比对。这里要特别注意幅度对齐如果中间采用的是逐级右移方案RTL输出等于真实FFT结果的1/16比较时要先把参考结果除以16。rtl_out np.loadtxt(fft_out.txt) # RTL仿真导出的复数结果 ref_out np.fft.fft(signal) # 参考FFT ref_scaled ref_out / 16 # 对齐RTL的缩放 error rtl_out - ref_scaled max_err np.max(np.abs(error)) rms_err np.sqrt(np.mean(np.abs(error)**2)) snr 20 * np.log10(np.linalg.norm(ref_scaled) / np.linalg.norm(error)) print(fmax error: {max_err:.4f}) print(frms error: {rms_err:.4f}) print(fSNR: {snr:.2f} dB)我实际做下来用8bit输入、逐级右移方案SNR大概能到50dB以上如果把截位改成四舍五入SNR还能再提高几个dB。如果你发现SNR只有30dB以下大概率不是位数不够而是旋转因子方向反了或者截位逻辑有bug。5.3 常见bug复现从仿真波形定位说两个我调试时真实踩过的坑供各位参考。坑1输出频谱顺序全乱表现是三个信号频率分量的位置完全对不上。排查时发现是倒位序地址写反了我在写RAM时用了bit_reverse4(addr[3:0])但addr在循环里已经是倒好的序等于反转了两次等于没反。修正方式是只保留第一次反转。坑2输出幅度明显偏小输入是100量级的信号输出结果却小了一截。原因在截位逻辑乘法器结果右移15位后符号位扩展被wire signed默认保留但我在加法器里用了无符号比较导致大负数被错误截断。解决方法是让蝶形运算单元的所有信号统一用signed类型并且右移用而不是。6. 实测中的几个坑与下一步扩展思考6.1 输入满摆幅时的溢出问题逐级右移方案理论上不会溢出但前提是每一级的右移都在加法和减法之后进行。如果你的数据路径是加法器和减法器共用然后单独做右移要注意先把蝶形输出完整算完再统一右移。否则高位数据在截位前就可能溢出。6.2 乘法器位宽与DSP资源配置16bit × 16bit乘法正好落在Xilinx DSP48E1的原生能力范围内一个DSP就能搞定。如果你的FPGA里DSP资源紧张也可以用移位加法和分布式算法替代乘法器但那个实现更复杂16点FFT规模不大建议还是直接用乘法器省心。综合之后查一下报告确认乘法器被正确推断成了DSP原语而不是被综合器展开成一堆LUT。6.3 向更大点数扩展的架构思路16点做完之后往32点、64点扩展的思路是现成的。基-2结构下每增加一级级数加1旋转因子表变大控制逻辑中的级数计数器和地址位宽相应扩展。从单蝶形复用切换到多蝶形流水线也只需把多个蝶形单元和对应的寄存器堆排成流水控制逻辑从状态机变成简单的地址生成器。如果你之后要处理多路并行FFT比如MIMO系统里每个通道都需要频谱分析单蝶形复用方案就力不从心了。那时可以在“空分”上做文章复制多个蝶形单元每个通道一组RAM控制逻辑共享用基地址偏移区分通道这样吞吐率可以成倍提高。我在实际项目中最终采用的是单蝶形复用逐级右移的组合综合后LUT消耗不到1000DSP使用1个Latency约100拍。在数据速率不高的场景下这个实现的功耗和面积都比同配置的IP核要清爽。每次做FFT相关设计时我都会先问自己一句到底需不需要IP核的通用性如果只是固定点数、固定位宽的单通道场景手写往往更香。
返回列表