一提到FFT很多FPGA工程師的習(xí)慣性動(dòng)作是打開Vivado或Quartus里的FFT IP核把點(diǎn)數(shù)一配、流模式一選、點(diǎn)一下Generate完事。但今年我在一個(gè)資源預(yù)算比較緊的小項(xiàng)目里只需要做16點(diǎn)FFT輸入是8bit定點(diǎn)數(shù)外部數(shù)據(jù)速率不高而核的Latency、DSP占用和時(shí)序行為總覺得有點(diǎn)“黑盒”。折騰了一圈之后我決定干脆用Verilog手寫一個(gè)基-2 16點(diǎn)FFT。這個(gè)規(guī)模不算大但正好能把FFT硬件化的全部關(guān)鍵問題——蝶形運(yùn)算結(jié)構(gòu)、旋轉(zhuǎn)因子量化、定點(diǎn)數(shù)位寬控制、數(shù)據(jù)調(diào)度、仿真對(duì)比——完整地走一遍。做完之后回頭看這件事的價(jià)值遠(yuǎn)不止“省了一個(gè)IP核”它讓我對(duì)FFT的數(shù)據(jù)流和FPGA數(shù)字信號(hào)處理的設(shè)計(jì)思路都有了更實(shí)的掌握。這篇文章就按我實(shí)際的實(shí)現(xiàn)路徑把算法結(jié)構(gòu)、RTL拆分、定位數(shù)處理、驗(yàn)證方法和踩過的坑一次講清楚。1. 先算清楚這筆賬16點(diǎn)FFT為什么值得手寫1.1 資源占用與延遲的直觀對(duì)比基-2 16點(diǎn)FFT的運(yùn)算量是固定的4級(jí)蝶形每級(jí)8個(gè)蝶形總共32次復(fù)數(shù)蝶形運(yùn)算。但同樣的算法映射到FPGA上的架構(gòu)差別非常大。全并行架構(gòu)每一級(jí)8個(gè)蝶形每個(gè)蝶形一個(gè)復(fù)數(shù)乘法器總共需要32個(gè)復(fù)數(shù)乘法器等價(jià)于128個(gè)實(shí)數(shù)乘法器邏輯規(guī)模和布線壓力都不小。單蝶形復(fù)用架構(gòu)只實(shí)現(xiàn)一個(gè)蝶形運(yùn)算單元用狀態(tài)機(jī)控制它依次完成32次運(yùn)算每次運(yùn)算只用一個(gè)復(fù)數(shù)乘法器資源占用極小。IP核以Xilinx FFT IP為例配16點(diǎn)、Radix-2 Burst I/O模式Latency通常在百拍量級(jí)資源則取決于你選的實(shí)現(xiàn)方式。IP核的好處是通用性強(qiáng)但代價(jià)是“黑盒”綜合之后未必是最優(yōu)解。我把三種方案在16點(diǎn)這個(gè)規(guī)模下的表現(xiàn)粗略列了一下方案復(fù)數(shù)乘法器數(shù)量大致Latency可控性適用場(chǎng)景全并行/流水線32數(shù)拍到十幾拍高高速流式處理單蝶形復(fù)用1約100拍很高低速、資源受限FFT IP核取決于配置數(shù)十到數(shù)百拍低通用快速集成1.2 什么情況下應(yīng)該放棄IP核用手寫IP核不是不能用而是有些場(chǎng)景下它確實(shí)不是最優(yōu)選擇。資源受限的小規(guī)模設(shè)計(jì)只做16點(diǎn)FFT用IP核有點(diǎn)“殺雞用牛刀”配置頁翻半天不說生成的邏輯可能比手寫大不少。時(shí)序與流水需要精確對(duì)齊IP核Latency是固定的但它內(nèi)部到底怎么調(diào)度、能不能和你的上下游邏輯無縫銜接很多時(shí)候要靠經(jīng)驗(yàn)去猜。手寫的話每一拍干什么都是自己控制的。需要自定義數(shù)據(jù)格式與截位策略FFT IP核的輸出位寬、縮放方式是可配的但如果你有特殊的定點(diǎn)數(shù)格式要求或者想在中間級(jí)做自定義處理手寫反而更靈活。學(xué)習(xí)與面試需求現(xiàn)在不少數(shù)字IC和FPGA相關(guān)的筆試面試題里都有“手撕FFT”這種題自己完整寫過一遍理解深度完全不同。16點(diǎn)FFT的手寫邏輯量并不大核心模塊也就幾個(gè)蝶形運(yùn)算單元、旋轉(zhuǎn)因子ROM、地址生成器、狀態(tài)機(jī)。在資源利用率上單蝶形方案用一個(gè)DSP48就夠邏輯消耗只有幾百個(gè)LUT。2. 基-2 16點(diǎn)FFT的蝶形結(jié)構(gòu)與旋轉(zhuǎn)因子表2.1 DIT還是DIF選哪一個(gè)更順手FFT的按時(shí)間抽取DIT和按頻率抽取DIF本質(zhì)上用同一個(gè)蝶形公式區(qū)別在于輸入輸出順序。DIT要求輸入按倒位序排列輸出是自然序DIF則是輸入自然序輸出倒位序。我推薦用DIT原因很直接輸入數(shù)據(jù)寫入RAM時(shí)做一次bit-reverse地址映射即可后續(xù)每級(jí)蝶形運(yùn)算都按自然序讀數(shù)據(jù)控制邏輯簡(jiǎn)單。DIF雖然輸入不需要倒序但最后輸出要處理倒序如果下游還要做頻譜分析或逆變換反而多一道麻煩。2.2 四級(jí)蝶形的排布規(guī)律16點(diǎn)基-2 FFT的蝶形間距和旋轉(zhuǎn)因子指數(shù)是有嚴(yán)格規(guī)律的。設(shè)級(jí)數(shù)從0開始記第s級(jí)的蝶形間距是2^s每級(jí)固定8個(gè)蝶形第s級(jí)第j個(gè)蝶形的旋轉(zhuǎn)因子指數(shù)由公式k j × (16 / (2^(s1))) 計(jì)算。實(shí)際算一遍就是這個(gè)表級(jí)數(shù)蝶形間距分組數(shù)旋轉(zhuǎn)因子指數(shù)說明第0級(jí)18k0每組1個(gè)蝶形第1級(jí)24k0,4每組2個(gè)蝶形第2級(jí)42k0,2,4,6每組4個(gè)蝶形第3級(jí)81k0,1,...,7一組8個(gè)蝶形這里有個(gè)容易記錯(cuò)的地方16點(diǎn)FFT第二級(jí)旋轉(zhuǎn)因子不是W16^0和W16^2而是W16^0和W16^4。因?yàn)榈诙?jí)里每組兩個(gè)蝶形合并的只是4點(diǎn)子序列旋轉(zhuǎn)因子對(duì)應(yīng)的是N/4周期上的角度。2.3 輸入倒位序的Verilog處理16點(diǎn)需要4位地址反轉(zhuǎn)。比如輸入序號(hào)1二進(jìn)制0001反轉(zhuǎn)后地址是1000也就是8。在RTL里寫一個(gè)反轉(zhuǎn)函數(shù)就行function [3:0] bit_reverse4; input [3:0] addr; begin bit_reverse4 {addr[0], addr[1], addr[2], addr[3]}; end endfunction實(shí)際寫入RAM時(shí)就把原始輸入數(shù)據(jù)寫到bit_reverse4(地址)的位置。例如第1個(gè)采樣點(diǎn)寫到地址8第2個(gè)寫到地址4第3個(gè)寫到地址12依次類推。這個(gè)映射可以在Testbench階段就預(yù)先確認(rèn)一遍常見的一個(gè)低級(jí)錯(cuò)誤就是倒位序只反轉(zhuǎn)了半字節(jié)或者反轉(zhuǎn)方向?qū)懛磳?dǎo)致最終輸出頻譜順序是亂的。3. 定點(diǎn)數(shù)格式與中間位寬最容易翻車的環(huán)節(jié)3.1 Q格式定點(diǎn)數(shù)的基礎(chǔ)FPGA里做不了浮點(diǎn)定點(diǎn)數(shù)格式得一開始就定好。最常用的是Q格式表示法。Q1.71位符號(hào)位7位小數(shù)范圍[-1, 1-2^-7]精度約0.0078。適合表示歸一化后的信號(hào)。Q1.151位符號(hào)位15位小數(shù)范圍[-1, 1-2^-15]精度非常高適合存旋轉(zhuǎn)因子的cos/sin值。在這個(gè)設(shè)計(jì)里輸入用Q1.7格式的8bit有符號(hào)數(shù)即可旋轉(zhuǎn)因子用Q1.15格式的16bit有符號(hào)數(shù)。兩個(gè)定點(diǎn)數(shù)相乘后結(jié)果的格式是Q(1.7) × Q(1.15) ≈ Q(1.22)需要根據(jù)需求截位。3.2 兩種中間位寬方案逐級(jí)右移 vs 全精度增長(zhǎng)這塊是整個(gè)設(shè)計(jì)里最容易被低估的部分。FFT每級(jí)蝶形進(jìn)行一次加法和一次減法數(shù)據(jù)位寬理論上每級(jí)增長(zhǎng)1bit4級(jí)下來最多增長(zhǎng)4bit。16點(diǎn)FFT的增益正好是16也就是說一個(gè)接近滿幅的8bit輸入經(jīng)過4級(jí)運(yùn)算后中間結(jié)果的理論最大值可以達(dá)到輸入的16倍。我試驗(yàn)過兩種主流方案方案A逐級(jí)右移1位每級(jí)蝶形計(jì)算完成后把結(jié)果右移1位相當(dāng)于每級(jí)做一次1/2縮放。4級(jí)之后輸出等于原始FFT結(jié)果的1/16但整個(gè)運(yùn)算過程中的數(shù)據(jù)位寬恒定不變不會(huì)溢出。實(shí)現(xiàn)簡(jiǎn)單資源占用小。代價(jià)是每級(jí)右移時(shí)都會(huì)引入量化誤差不過16點(diǎn)FFT的級(jí)數(shù)只有4級(jí)誤差累積有限實(shí)測(cè)SNR通常還有50dB以上完全夠用。方案B全精度增長(zhǎng)最后統(tǒng)一截位輸入8bit每級(jí)把位寬擴(kuò)到16bit甚至更寬中間不縮放最后輸出時(shí)再截位或飽和。精度最高但存儲(chǔ)和運(yùn)算邏輯變多而且最后截位時(shí)如果處理不好照樣會(huì)有溢出風(fēng)險(xiǎn)。兩種方案我建議初學(xué)者直接用方案A理由很簡(jiǎn)單省心、無溢出、邏輯清晰。如果你的應(yīng)用對(duì)SNR要求特別高再考慮方案B。3.3 旋轉(zhuǎn)因子的量化與ROM生成旋轉(zhuǎn)因子的計(jì)算公式是W16^k cos(2πk/16) ? j·sin(2πk/16)在Verilog里存成兩個(gè)16bit有符號(hào)數(shù)實(shí)部存cos值虛部存?sin值?；?2 16點(diǎn)FFT實(shí)際需要用的旋轉(zhuǎn)因子是k0到7另一半可以通過對(duì)稱性得到所以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)需要注意生成時(shí)要用16bit有符號(hào)數(shù)的補(bǔ)碼形式存這樣才能在Verilog里直接用signed類型讀取。4. RTL實(shí)現(xiàn)的關(guān)鍵模塊蝶形單元、ROM與地址生成4.1 蝶形運(yùn)算單元的Verilog寫法蝶形運(yùn)算是FFT的核心設(shè)輸入為x和y旋轉(zhuǎn)因子為W則輸出為x y·W和x ? y·W。復(fù)數(shù)乘法展開后需要4個(gè)實(shí)數(shù)乘法和幾次加法。寫成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 的實(shí)部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這里有個(gè)細(xì)節(jié)乘法結(jié)果是32bit使用算數(shù)右移把它截回16bit。右移時(shí)研究一下仿真結(jié)果看看用截?cái)噙€是四舍五入對(duì)SNR影響大。四舍五入可以在右移前加一個(gè)偏移量實(shí)現(xiàn)代價(jià)是幾個(gè)加法器。4.2 旋轉(zhuǎn)因子ROM的初始化ROM可以直接用Verilog數(shù)組加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一次只能讀一個(gè)數(shù)組實(shí)踐中我通常把實(shí)部和虛部分成兩個(gè)文件或者用initial塊里的case語句直接查表。16點(diǎn)規(guī)模很小case語句反而最直觀可讀性也最好。4.3 地址生成與狀態(tài)機(jī)控制單蝶形復(fù)用方案需要一個(gè)狀態(tài)機(jī)來調(diào)度32次蝶形運(yùn)算。核心控制邏輯是兩級(jí)計(jì)數(shù)器stage_cnt[1:0]當(dāng)前級(jí)數(shù)0到3。bfly_cnt[2:0]當(dāng)前級(jí)內(nèi)的蝶形序號(hào)0到7。根據(jù)這兩個(gè)計(jì)數(shù)器和2.2節(jié)的規(guī)律就可以算出當(dāng)前蝶形的兩個(gè)輸入數(shù)據(jù)地址和旋轉(zhuǎn)因子地址wire [3:0] spacing 4b0001 stage_cnt; // 蝶形間距 wire [3:0] j bfly_cnt (spacing - 4b0001); // 組內(nèi)蝶形序號(hào) 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));三級(jí)條件嵌套就能把旋轉(zhuǎn)因子換算出來。狀態(tài)機(jī)則負(fù)責(zé)在每個(gè)蝶形運(yùn)算開始時(shí)讀兩個(gè)復(fù)數(shù)計(jì)算完成后寫回RAM然后切換下一組地址。5. 仿真驗(yàn)證從激勵(lì)生成到誤差分析的全鏈路5.1 用Python生成測(cè)試向量純隨機(jī)數(shù)測(cè)試能驗(yàn)證功能但看不出精度損失。我習(xí)慣用已知頻率的正弦信號(hào)做激勵(lì)這樣頻譜圖上一眼就能看出峰值位置對(duì)不對(duì)。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只測(cè)實(shí)信號(hào) with open(input_imag.txt, w) as f: for _ in range(N): f.write(00\n)Testbench里用$readmemh把這兩個(gè)文件讀入存儲(chǔ)數(shù)組然后啟動(dòng)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 結(jié)果比對(duì)與SNR計(jì)算仿真結(jié)束后把輸出結(jié)果導(dǎo)出在Python里與numpy.fft.fft的參考結(jié)果比對(duì)。這里要特別注意幅度對(duì)齊如果中間采用的是逐級(jí)右移方案RTL輸出等于真實(shí)FFT結(jié)果的1/16比較時(shí)要先把參考結(jié)果除以16。rtl_out np.loadtxt(fft_out.txt) # RTL仿真導(dǎo)出的復(fù)數(shù)結(jié)果 ref_out np.fft.fft(signal) # 參考FFT ref_scaled ref_out / 16 # 對(duì)齊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)我實(shí)際做下來用8bit輸入、逐級(jí)右移方案SNR大概能到50dB以上如果把截位改成四舍五入SNR還能再提高幾個(gè)dB。如果你發(fā)現(xiàn)SNR只有30dB以下大概率不是位數(shù)不夠而是旋轉(zhuǎn)因子方向反了或者截位邏輯有bug。5.3 常見bug復(fù)現(xiàn)從仿真波形定位說兩個(gè)我調(diào)試時(shí)真實(shí)踩過的坑供各位參考。坑1輸出頻譜順序全亂表現(xiàn)是三個(gè)信號(hào)頻率分量的位置完全對(duì)不上。排查時(shí)發(fā)現(xiàn)是倒位序地址寫反了我在寫RAM時(shí)用了bit_reverse4(addr[3:0])但addr在循環(huán)里已經(jīng)是倒好的序等于反轉(zhuǎn)了兩次等于沒反。修正方式是只保留第一次反轉(zhuǎn)。坑2輸出幅度明顯偏小輸入是100量級(jí)的信號(hào)輸出結(jié)果卻小了一截。原因在截位邏輯乘法器結(jié)果右移15位后符號(hào)位擴(kuò)展被wire signed默認(rèn)保留但我在加法器里用了無符號(hào)比較導(dǎo)致大負(fù)數(shù)被錯(cuò)誤截?cái)?。解決方法是讓蝶形運(yùn)算單元的所有信號(hào)統(tǒng)一用signed類型并且右移用而不是。6. 實(shí)測(cè)中的幾個(gè)坑與下一步擴(kuò)展思考6.1 輸入滿擺幅時(shí)的溢出問題逐級(jí)右移方案理論上不會(huì)溢出但前提是每一級(jí)的右移都在加法和減法之后進(jìn)行。如果你的數(shù)據(jù)路徑是加法器和減法器共用然后單獨(dú)做右移要注意先把蝶形輸出完整算完再統(tǒng)一右移。否則高位數(shù)據(jù)在截位前就可能溢出。6.2 乘法器位寬與DSP資源配置16bit × 16bit乘法正好落在Xilinx DSP48E1的原生能力范圍內(nèi)一個(gè)DSP就能搞定。如果你的FPGA里DSP資源緊張也可以用移位加法和分布式算法替代乘法器但那個(gè)實(shí)現(xiàn)更復(fù)雜16點(diǎn)FFT規(guī)模不大建議還是直接用乘法器省心。綜合之后查一下報(bào)告確認(rèn)乘法器被正確推斷成了DSP原語而不是被綜合器展開成一堆LUT。6.3 向更大點(diǎn)數(shù)擴(kuò)展的架構(gòu)思路16點(diǎn)做完之后往32點(diǎn)、64點(diǎn)擴(kuò)展的思路是現(xiàn)成的?；?2結(jié)構(gòu)下每增加一級(jí)級(jí)數(shù)加1旋轉(zhuǎn)因子表變大控制邏輯中的級(jí)數(shù)計(jì)數(shù)器和地址位寬相應(yīng)擴(kuò)展。從單蝶形復(fù)用切換到多蝶形流水線也只需把多個(gè)蝶形單元和對(duì)應(yīng)的寄存器堆排成流水控制邏輯從狀態(tài)機(jī)變成簡(jiǎn)單的地址生成器。如果你之后要處理多路并行FFT比如MIMO系統(tǒng)里每個(gè)通道都需要頻譜分析單蝶形復(fù)用方案就力不從心了。那時(shí)可以在“空分”上做文章復(fù)制多個(gè)蝶形單元每個(gè)通道一組RAM控制邏輯共享用基地址偏移區(qū)分通道這樣吞吐率可以成倍提高。我在實(shí)際項(xiàng)目中最終采用的是單蝶形復(fù)用逐級(jí)右移的組合綜合后LUT消耗不到1000DSP使用1個(gè)Latency約100拍。在數(shù)據(jù)速率不高的場(chǎng)景下這個(gè)實(shí)現(xiàn)的功耗和面積都比同配置的IP核要清爽。每次做FFT相關(guān)設(shè)計(jì)時(shí)我都會(huì)先問自己一句到底需不需要IP核的通用性如果只是固定點(diǎn)數(shù)、固定位寬的單通道場(chǎng)景手寫往往更香。