FFT(Fast Fourier Transform),快速傅立葉變換,是一種 DFT(離散傅里葉變換)的高效算法。在以時(shí)頻變換分析為基礎(chǔ)的數(shù)字處理方法中,有著不可替代的作用。
FFT 原理
◆公式推導(dǎo)
DFT 的運(yùn)算公式為:
其中,
將離散傅里葉變換公式拆分成奇偶項(xiàng),則前 N/2 個(gè)點(diǎn)可以表示為:
同理,后 N/2 個(gè)點(diǎn)可以表示為:
由此可知,后 N/2 個(gè)點(diǎn)的值完全可以通過計(jì)算前 N/2 個(gè)點(diǎn)時(shí)的中間過程值確定。對(duì) A[k] 與 B[k] 繼續(xù)進(jìn)行奇偶分解,直至變成 2 點(diǎn)的 DFT,這樣就可以避免很多的重復(fù)計(jì)算,實(shí)現(xiàn)了快速離散傅里葉變換(FFT)的過程。
◆算法結(jié)構(gòu)
8 點(diǎn) FFT 計(jì)算的結(jié)構(gòu)示意圖如下。
由圖可知,只需要簡(jiǎn)單的計(jì)算幾次乘法和加法,便可完成離散傅里葉變換過程,而不是對(duì)每個(gè)數(shù)據(jù)進(jìn)行繁瑣的相乘和累加。
◆重要特性
(1) 級(jí)的概念
每分割一次,稱為一級(jí)運(yùn)算。
設(shè) FFT 運(yùn)算點(diǎn)數(shù)為 N,共有 M 級(jí)運(yùn)算,則它們滿足:
每一級(jí)運(yùn)算的標(biāo)識(shí)為 m = 0, 1, 2, ..., M-1。
為了便于分割計(jì)算,F(xiàn)FT 點(diǎn)數(shù) N 的取值經(jīng)常為 2 的整數(shù)次冪。
(2) 蝶形單元
FFT 計(jì)算結(jié)構(gòu)由若干個(gè)蝶形運(yùn)算單元組成,每個(gè)運(yùn)算單元示意圖如下:
蝶形單元的輸入輸出滿足:
其中,。
每一個(gè)蝶形單元運(yùn)算時(shí),進(jìn)行了一次乘法和兩次加法。
每一級(jí)中,均有 N/2 個(gè)蝶形單元。
故完成一次 FFT 所需要的乘法次數(shù)和加法次數(shù)分別為:
(3) 組的概念
每一級(jí) N/2 個(gè)蝶形單元可分為若干組,每一組有著相同的結(jié)構(gòu)與因子分布。
例如 m=0 時(shí),可以分為 N/2=4 組。
m=1 時(shí),可以分為 N/4=2 組。
m=M-1 時(shí),此時(shí)只能分為 1 組。
(4) 因子分布
因子存在于 m 級(jí),其中
。
在 8 點(diǎn) FFT 第二級(jí)運(yùn)算中,即 m=1 ,蝶形運(yùn)算因子可以化簡(jiǎn)為:
(5) 碼位倒置
對(duì)于 N=8 點(diǎn)的 FFT 計(jì)算,X(0) ~ X(7) 位置對(duì)應(yīng)的 2 進(jìn)制碼為:
X(000), X(001), X(010), X(011), X(100), X(101), X(110), X(111)
將其位置的 2 進(jìn)制碼進(jìn)行翻轉(zhuǎn):
X(000), X(100), X(010), X(110), X(001), X(101), X(011), X(111)
此時(shí)位置對(duì)應(yīng)的 10 進(jìn)制為:
X(0), X(4), X(2), X(6), X(1), X(5), X(3), X(7)
恰好對(duì)應(yīng) FFT 第一級(jí)輸入數(shù)據(jù)的順序。
該特性有利于 FFT 的編程實(shí)現(xiàn)。
FFT 設(shè)計(jì)
◆設(shè)計(jì)說明
為了利用仿真簡(jiǎn)單的說明 FFT 的變換過程,數(shù)據(jù)點(diǎn)數(shù)取較小的值 8。
如果數(shù)據(jù)是串行輸入,需要先進(jìn)行緩存,所以設(shè)計(jì)時(shí)數(shù)據(jù)輸入方式為并行。
數(shù)據(jù)輸入分為實(shí)部和虛部共 2 部分,所以計(jì)算結(jié)果也分為實(shí)部和虛部。
設(shè)計(jì)采用流水結(jié)構(gòu),暫不考慮資源消耗的問題。
為了使設(shè)計(jì)結(jié)構(gòu)更加簡(jiǎn)單,這里做一步妥協(xié),乘法計(jì)算直接使用乘號(hào)。如果 FFT 設(shè)計(jì)應(yīng)用于實(shí)際,一定要將乘法結(jié)構(gòu)換成可以流水的乘法器,或使用官方提供的效率較高的乘法器 IP。
◆蝶形單元設(shè)計(jì)
蝶形單元為定點(diǎn)運(yùn)算,需要對(duì)旋轉(zhuǎn)因子進(jìn)行定點(diǎn)量化。
借助 matlab 將旋轉(zhuǎn)因子擴(kuò)大 8192 倍(左移 13 位),可參考附錄。
為了防止蝶形運(yùn)算中的乘法和加法導(dǎo)致位寬逐級(jí)增大,每一級(jí)運(yùn)算完成后,要對(duì)輸出數(shù)據(jù)進(jìn)行固定位寬的截位,也可去掉旋轉(zhuǎn)因子倍數(shù)增大而帶來的影響。
代碼如下:
`timescale 1ns/100ps
/**************** butter unit *************************
Xm(p) ------------------------ > Xm+1(p)
- - >
- -
-
- -
- - >
Xm(q) ------------------------ > Xm+1(q)
Wn -1
*//////////////////////////////////////////////////////
module butterfly
(
input clk,
input rstn,
input en,
input signed [23:0] xp_real, // Xm(p)
input signed [23:0] xp_imag,
input signed [23:0] xq_real, // Xm(q)
input signed [23:0] xq_imag,
input signed [15:0] factor_real, // Wnr
input signed [15:0] factor_imag,
output valid,
output signed [23:0] yp_real, //Xm+1(p)
output signed [23:0] yp_imag,
output signed [23:0] yq_real, //Xm+1(q)
output signed [23:0] yq_imag);
reg [4:0] en_r ;
always @(posedge clk or negedge rstn) begin
if (!rstn) begin
en_r <= 'b0 ;
end
else begin
en_r <= {en_r[3:0], en} ;
end
end
//=====================================================//
//(1.0) Xm(q) mutiply and Xm(p) delay
reg signed [39:0] xq_wnr_real0;
reg signed [39:0] xq_wnr_real1;
reg signed [39:0] xq_wnr_imag0;
reg signed [39:0] xq_wnr_imag1;
reg signed [39:0] xp_real_d;
reg signed [39:0] xp_imag_d;
always @(posedge clk or negedge rstn) begin
if (!rstn) begin
xp_real_d <= 'b0;
xp_imag_d <= 'b0;
xq_wnr_real0 <= 'b0;
xq_wnr_real1 <= 'b0;
xq_wnr_imag0 <= 'b0;
xq_wnr_imag1 <= 'b0;
end
else if (en) begin
xq_wnr_real0 <= xq_real * factor_real;
xq_wnr_real1 <= xq_imag * factor_imag;
xq_wnr_imag0 <= xq_real * factor_imag;
xq_wnr_imag1 <= xq_imag * factor_real;
//expanding 8192 times as Wnr
xp_real_d <= {{4{xp_real[23]}}, xp_real[22:0], 13'b0};
xp_imag_d <= {{4{xp_imag[23]}}, xp_imag[22:0], 13'b0};
end
end
//(1.1) get Xm(q) mutiplied-results and Xm(p) delay again
reg signed [39:0] xp_real_d1;
reg signed [39:0] xp_imag_d1;
reg signed [39:0] xq_wnr_real;
reg signed [39:0] xq_wnr_imag;
always @(posedge clk or negedge rstn) begin
if (!rstn) begin
xp_real_d1 <= 'b0;
xp_imag_d1 <= 'b0;
xq_wnr_real <= 'b0 ;
xq_wnr_imag <= 'b0 ;
end
else if (en_r[0]) begin
xp_real_d1 <= xp_real_d;
xp_imag_d1 <= xp_imag_d;
//提前設(shè)置好位寬余量,防止數(shù)據(jù)溢出
xq_wnr_real <= xq_wnr_real0 - xq_wnr_real1 ;
xq_wnr_imag <= xq_wnr_imag0 + xq_wnr_imag1 ;
end
end
//======================================================//
//(2.0) butter results
reg signed [39:0] yp_real_r;
reg signed [39:0] yp_imag_r;
reg signed [39:0] yq_real_r;
reg signed [39:0] yq_imag_r;
always @(posedge clk or negedge rstn) begin
if (!rstn) begin