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)的值完全可以通過(guò)計(jì)算前 N/2 個(gè)點(diǎn)時(shí)的中間過(guò)程值確定。 對(duì) A[k] 與 B[k] 繼續(xù)進(jìn)行奇偶分解，直至變成 2 點(diǎn)的 DFT，這樣就可以避免很多的重復(fù)計(jì)算，實(shí)現(xiàn)了快速離散傅里葉變換（FFT）的過(guò)程。

算法結(jié)構(gòu)

8 點(diǎn) FFT 計(jì)算的結(jié)構(gòu)示意圖如下。

由圖可知，只需要簡(jiǎn)單的計(jì)算幾次乘法和加法，便可完成離散傅里葉變換過(guò)程，而不是對(duì)每個(gè)數(shù)據(jù)進(jìn)行繁瑣的相乘和累加。

重要特性

(1) 級(jí)的概念

每分割一次，稱(chēng)為一級(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ì)說(shuō)明

為了利用仿真簡(jiǎn)單的說(shuō)明 FFT 的變換過(guò)程，數(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)，暫不考慮資源消耗的問(wèn)題。

為了使設(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ù)增大而帶來(lái)的影響。

代碼如下：

`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
            yp_real_r      <= 'b0;
            yp_imag_r      <= 'b0;
            yq_real_r      <= 'b0;
            yq_imag_r      <= 'b0;
        end
        else if (en_r[1]) begin
            yp_real_r      <= xp_real_d1 + xq_wnr_real;
            yp_imag_r      <= xp_imag_d1 + xq_wnr_imag;
            yq_real_r      <= xp_real_d1 - xq_wnr_real;
            yq_imag_r      <= xp_imag_d1 - xq_wnr_imag;
        end
    end


    //(3) discard the low 13bits because of Wnr
    assign yp_real = {yp_real_r[39], yp_real_r[13+23:13]};
    assign yp_imag = {yp_imag_r[39], yp_imag_r[13+23:13]};
    assign yq_real = {yq_real_r[39], yq_real_r[13+23:13]};
    assign yq_imag = {yq_imag_r[39], yq_imag_r[13+23:13]};
    assign valid   = en_r[2];

endmodule

頂層例化

根據(jù) FFT 算法結(jié)構(gòu)示意圖，將蝶形單元例化，完成最后的 FFT 功能。

可根據(jù)每一級(jí)蝶形單元的輸入輸出對(duì)應(yīng)關(guān)系，依次手動(dòng)例化 12 次，也可利用 generate 進(jìn)行例化，此時(shí)就需要非常熟悉 FFT 中“組”和“級(jí)”的特點(diǎn)：

(1) 8 點(diǎn) FFT 設(shè)計(jì)，需要 3 級(jí)運(yùn)算，每一級(jí)有 4 個(gè)蝶形單元，每一級(jí)的組數(shù)目分別是 4、2、1。

(2) 每一級(jí)的組內(nèi)一個(gè)蝶形單元中兩個(gè)輸入端口的距離恒為 （m 為級(jí)標(biāo)號(hào)，對(duì)應(yīng)左移運(yùn)算 1<<< span="">m），組內(nèi)兩個(gè)蝶形單元的第一個(gè)輸入端口間的距離為 1。

(3) 每一級(jí)相鄰組間的第一個(gè)蝶形單元的第一個(gè)輸入端口的距離為 （對(duì)應(yīng)左移運(yùn)算 2<<< span="">m）。

例化代碼如下。

其中，矩陣信號(hào) xm_real（xm_imag）的一維、二維地址是代表級(jí)和組的標(biāo)識(shí)。

在判斷信號(hào)端口之間的連接關(guān)系時(shí)，使用了看似復(fù)雜的判斷邏輯，而且還帶有乘號(hào)，其實(shí)最終生成的電路和手動(dòng)編寫(xiě)代碼例化 12 個(gè)蝶形單元的方式是完全相同的。 因?yàn)?generate 中的變量只是輔助生成實(shí)際的電路，相關(guān)值的計(jì)算判斷都已經(jīng)在編譯時(shí)完成。 這些變量更不會(huì)生成實(shí)際的電路，只是為更快速的模塊例化提供了一種方法。

timescale 1ns/100ps
module fft8 (
    input                    clk,
    input                    rstn,
    input                    en,


    input signed [23:0]      x0_real,
    input signed [23:0]      x0_imag,
    input signed [23:0]      x1_real,
    input signed [23:0]      x1_imag,
    input signed [23:0]      x2_real,
    input signed [23:0]      x2_imag,
    input signed [23:0]      x3_real,
    input signed [23:0]      x3_imag,
    input signed [23:0]      x4_real,
    input signed [23:0]      x4_imag,
    input signed [23:0]      x5_real,
    input signed [23:0]      x5_imag,
    input signed [23:0]      x6_real,
    input signed [23:0]      x6_imag,
    input signed [23:0]      x7_real,
    input signed [23:0]      x7_imag,


    output                   valid,
    output signed [23:0]     y0_real,
    output signed [23:0]     y0_imag,
    output signed [23:0]     y1_real,
    output signed [23:0]     y1_imag,
    output signed [23:0]     y2_real,
    output signed [23:0]     y2_imag,
    output signed [23:0]     y3_real,
    output signed [23:0]     y3_imag,
    output signed [23:0]     y4_real,
    output signed [23:0]     y4_imag,
    output signed [23:0]     y5_real,
    output signed [23:0]     y5_imag,
    output signed [23:0]     y6_real,
    output signed [23:0]     y6_imag,
    output signed [23:0]     y7_real,
    output signed [23:0]     y7_imag
    );


    //operating data
    wire signed [23:0]             xm_real [3:0] [7:0];
    wire signed [23:0]             xm_imag [3:0] [7:0];
    wire                           en_connect [15:0] ;
    assign                         en_connect[0] = en;
    assign                         en_connect[1] = en;
    assign                         en_connect[2] = en;
    assign                         en_connect[3] = en;


    //factor, multiplied by 0x2000
    wire signed [15:0]             factor_real [3:0] ;
    wire signed [15:0]             factor_imag [3:0];
    assign factor_real[0]        = 16'h2000; //1
    assign factor_imag[0]        = 16'h0000; //0
    assign factor_real[1]        = 16'h16a0; //sqrt(2)/2
    assign factor_imag[1]        = 16'he95f; //-sqrt(2)/2
    assign factor_real[2]        = 16'h0000; //0
    assign factor_imag[2]        = 16'he000; //-1
    assign factor_real[3]        = 16'he95f; //-sqrt(2)/2
    assign factor_imag[3]        = 16'he95f; //-sqrt(2)/2


    //輸入初始化，和碼位有關(guān)倒置
    assign xm_real[0][0] = x0_real;
    assign xm_real[0][1] = x4_real;
    assign xm_real[0][2] = x2_real;
    assign xm_real[0][3] = x6_real;
    assign xm_real[0][4] = x1_real;
    assign xm_real[0][5] = x5_real;
    assign xm_real[0][6] = x3_real;
    assign xm_real[0][7] = x7_real;
    assign xm_imag[0][0] = x0_imag;
    assign xm_imag[0][1] = x4_imag;
    assign xm_imag[0][2] = x2_imag;
    assign xm_imag[0][3] = x6_imag;
    assign xm_imag[0][4] = x1_imag;
    assign xm_imag[0][5] = x5_imag;
    assign xm_imag[0][6] = x3_imag;
    assign xm_imag[0][7] = x7_imag;


    //butter instantiaiton
    //integer              index[11:0] ;
    genvar               m, k;
    generate
    //3 stage
    for(m=0; m<=2; m=m+1) begin: stage
        for (k=0; k<=3; k=k+1) begin: unit


            butterfly           u_butter(
               .clk        (clk                 ) ,
               .rstn       (rstn                ) ,
               .en         (en_connect[m*4 + k] ) ,
                       //是否再組內(nèi)？組編號(hào)+組內(nèi)編號(hào)：下組編號(hào)+新組內(nèi)編號(hào)
               .xp_real    (xm_real[ m ] [k[m:0] < (1<3:m] << (m+1)) + k[m:0] :
                           (k[3:m] << (m+1)) + (k[m:0]-(1<0] < (1<3:m] << (m+1)) + k[m:0] :
                           (k[3:m] << (m+1)) + (k[m:0]-(1<0] < (1<3:m] << (m+1)) + k[m:0] :
                           (k[3:m] << (m+1)) + (k[m:0]-(1<1<

測(cè)試平臺(tái)

testbench 編寫(xiě)如下，主要用于 16 路實(shí)、復(fù)數(shù)據(jù)的連續(xù)輸入。 因?yàn)槊看?FFT 只有 8 點(diǎn)數(shù)據(jù)，所以送入的數(shù)據(jù)比較隨意，并不是正弦波等規(guī)則的數(shù)據(jù)。

`timescale 1ns/100ps
module test ;
    reg          clk;
    reg          rstn;
    reg          en ;


    reg signed   [23:0]   x0_real;
    reg signed   [23:0]   x0_imag;
    reg signed   [23:0]   x1_real;
    reg signed   [23:0]   x1_imag;
    reg signed   [23:0]   x2_real;
    reg signed   [23:0]   x2_imag;
    reg signed   [23:0]   x3_real;
    reg signed   [23:0]   x3_imag;
    reg signed   [23:0]   x4_real;
    reg signed   [23:0]   x4_imag;
    reg signed   [23:0]   x5_real;
    reg signed   [23:0]   x5_imag;
    reg signed   [23:0]   x6_real;
    reg signed   [23:0]   x6_imag;
    reg signed   [23:0]   x7_real;
    reg signed   [23:0]   x7_imag;


    wire                  valid;
    wire signed  [23:0]   y0_real;
    wire signed  [23:0]   y0_imag;
    wire signed  [23:0]   y1_real;
    wire signed  [23:0]   y1_imag;
    wire signed  [23:0]   y2_real;
    wire signed  [23:0]   y2_imag;
    wire signed  [23:0]   y3_real;
    wire signed  [23:0]   y3_imag;
    wire signed  [23:0]   y4_real;
    wire signed  [23:0]   y4_imag;
    wire signed  [23:0]   y5_real;
    wire signed  [23:0]   y5_imag;
    wire signed  [23:0]   y6_real;
    wire signed  [23:0]   y6_imag;
    wire signed  [23:0]   y7_real;
    wire signed  [23:0]   y7_imag;


    initial begin
        clk = 0; //50MHz
        rstn = 0 ;
        #10 rstn = 1;
        forever begin
            #10 clk = ~clk; //50MHz
        end
    end


    fft8 u_fft (
      .clk        (clk    ),
      .rstn       (rstn    ),
      .en         (en     ),
      .x0_real    (x0_real),
      .x0_imag    (x0_imag),
      .x1_real    (x1_real),
      .x1_imag    (x1_imag),
      .x2_real    (x2_real),
      .x2_imag    (x2_imag),
      .x3_real    (x3_real),
      .x3_imag    (x3_imag),
      .x4_real    (x4_real),
      .x4_imag    (x4_imag),
      .x5_real    (x5_real),
      .x5_imag    (x5_imag),
      .x6_real    (x6_real),
      .x6_imag    (x6_imag),
      .x7_real    (x7_real),
      .x7_imag    (x7_imag),


      .valid      (valid),
      .y0_real    (y0_real),
      .y0_imag    (y0_imag),
      .y1_real    (y1_real),
      .y1_imag    (y1_imag),
      .y2_real    (y2_real),
      .y2_imag    (y2_imag),
      .y3_real    (y3_real),
      .y3_imag    (y3_imag),
      .y4_real    (y4_real),
      .y4_imag    (y4_imag),
      .y5_real    (y5_real),
      .y5_imag    (y5_imag),
      .y6_real    (y6_real),
      .y6_imag    (y6_imag),
      .y7_real    (y7_real),
      .y7_imag    (y7_imag));


    //data input
    initial      begin
        en = 0 ;
        x0_real = 24'd10;
        x1_real = 24'd20;
        x2_real = 24'd30;
        x3_real = 24'd40;
        x4_real = 24'd10;
        x5_real = 24'd20;
        x6_real = 24'd30;
        x7_real = 24'd40;


        x0_imag = 24'd0;
        x1_imag = 24'd0;
        x2_imag = 24'd0;
        x3_imag = 24'd0;
        x4_imag = 24'd0;
        x5_imag = 24'd0;
        x6_imag = 24'd0;
        x7_imag = 24'd0;
        @(negedge clk) ;
        en = 1 ;
        forever begin
            @(negedge clk) ;
            x0_real = (x0_real > 22'h3F_ffff) ? 'b0 : x0_real + 1 ;
            x1_real = (x1_real > 22'h3F_ffff) ? 'b0 : x1_real + 1 ;
            x2_real = (x2_real > 22'h3F_ffff) ? 'b0 : x2_real + 31 ;
            x3_real = (x3_real > 22'h3F_ffff) ? 'b0 : x3_real + 1 ;
            x4_real = (x4_real > 22'h3F_ffff) ? 'b0 : x4_real + 23 ;
            x5_real = (x5_real > 22'h3F_ffff) ? 'b0 : x5_real + 1 ;
            x6_real = (x6_real > 22'h3F_ffff) ? 'b0 : x6_real + 6 ;
            x7_real = (x7_real > 22'h3F_ffff) ? 'b0 : x7_real + 1 ;


            x0_imag = (x0_imag > 22'h3F_ffff) ? 'b0 : x0_imag + 2 ;
            x1_imag = (x1_imag > 22'h3F_ffff) ? 'b0 : x1_imag + 5 ;
            x2_imag = (x2_imag > 22'h3F_ffff) ? 'b0 : x2_imag + 3 ;
            x3_imag = (x3_imag > 22'h3F_ffff) ? 'b0 : x3_imag + 6 ;
            x4_imag = (x4_imag > 22'h3F_ffff) ? 'b0 : x4_imag + 4 ;
            x5_imag = (x5_imag > 22'h3F_ffff) ? 'b0 : x5_imag + 8 ;
            x6_imag = (x6_imag > 22'h3F_ffff) ? 'b0 : x6_imag + 11 ;
            x7_imag = (x7_imag > 22'h3F_ffff) ? 'b0 : x7_imag + 7 ;
        end
    end


   //finish simulation
   initial #1000       $finish ;
endmodule

仿真結(jié)果

大致可以看出，F(xiàn)FT 結(jié)果可以流水輸出。

用 matlab 自帶的 FFT 函數(shù)對(duì)相同數(shù)據(jù)進(jìn)行運(yùn)算，前 2 組數(shù)據(jù) FFT 結(jié)果如下。

可以看出，第一次輸入的數(shù)據(jù)信號(hào)只有實(shí)部有效時(shí)，F(xiàn)FT 結(jié)果是完全一樣的。

但是第二次輸入的數(shù)據(jù)復(fù)部也有信號(hào)，此時(shí)兩者之間的結(jié)果開(kāi)始有誤差，有時(shí)誤差還很大。

用 matlab 對(duì) Verilog 實(shí)現(xiàn)的 FFT 過(guò)程進(jìn)行模擬，發(fā)現(xiàn)此過(guò)程的 FFT 結(jié)果和 Verilog 實(shí)現(xiàn)的 FFT 結(jié)果基本一致。

將有誤差的兩種 FFT 結(jié)果取絕對(duì)值進(jìn)行比較，圖示如下。

可以看出，F(xiàn)FT 結(jié)果的趨勢(shì)大致相同，但在個(gè)別點(diǎn)有肉眼可見(jiàn)的誤差。

設(shè)計(jì)總結(jié)：

就如設(shè)計(jì)蝶形單元時(shí)所說(shuō)，旋轉(zhuǎn)因子量化時(shí)，位寬的選擇就會(huì)引入誤差。

而且每個(gè)蝶形單元的運(yùn)算結(jié)果都會(huì)進(jìn)行截取，也會(huì)引入誤差。

matlab 計(jì)算 FFT 時(shí)不用考慮精度問(wèn)題，以其最高精度對(duì)數(shù)據(jù)進(jìn)行 FFT 計(jì)算。

以上所述，都會(huì)導(dǎo)致最后兩種 FFT 計(jì)算方式結(jié)果的差異。

感興趣的學(xué)者，可以將旋轉(zhuǎn)因子和輸入數(shù)據(jù)位寬再進(jìn)行一定的增加，F(xiàn)FT 點(diǎn)數(shù)也可以增加，然后再進(jìn)行仿真對(duì)比，相對(duì)誤差應(yīng)該會(huì)減小。

附錄：matlab 使用

生成旋轉(zhuǎn)因子

8 點(diǎn) FFT 只需要用到 4 個(gè)旋轉(zhuǎn)因子。 旋轉(zhuǎn)因子擴(kuò)大倍數(shù)為 8192。

clear all;close all;clc;
%=======================================================
% Wnr calcuting
%=======================================================
for r = 0:3 
    Wnr_factor  = cos(pi/4*r) - j*sin(pi/4*r) ;
    Wnr_integer = floor(Wnr_factor * 2^13) ;

    if (real(Wnr_integer)<0) 
        Wnr_real    = real(Wnr_integer) + 2^16 ; %負(fù)數(shù)的補(bǔ)碼
    else
        Wnr_real    = real(Wnr_integer) ;
    end
    if (imag(Wnr_integer)<0) 
        Wnr_imag    = imag(Wnr_integer) + 2^16 ; 
    else
        Wnr_imag    = imag(Wnr_integer);
    end

    Wnr(2*r+1,:)  =  dec2hex(Wnr_real)   %實(shí)部
    Wnr(2*r+2,:)  =  dec2hex(Wnr_imag)   %虛部
end

FFT 結(jié)果對(duì)比

matlab 模擬 Verilog 實(shí)現(xiàn) FFT 的過(guò)程如下，也包括 2 種 FFT 結(jié)果的對(duì)比。

clear all;close all;clc;
%=======================================================
% 8dots fft
%=======================================================
for r=0:3
    Wnr(r+1)  = cos(pi/4*r) - j*sin(pi/4*r) ;
end
x       = [10, 20, 30, 40, 10, 20 ,30 ,40];
step    = [1+2j, 1+5j, 31+3j, 1+6j, 23+4j, 1+8j, 6+11j, 1+7j];
x2      = x + step;
xm0     = [x2(0+1), x2(4+1), x2(2+1), x2(6+1), x2(1+1), x2(5+1),         x2(3+1), x2(7+1)] ;


%% stage1 
xm1(1) = xm0(1) + xm0(2)*Wnr(1) ;
xm1(2) = xm0(1) - xm0(2)*Wnr(1) ;
xm1(3) = xm0(3) + xm0(4)*Wnr(1) ;
xm1(4) = xm0(3) - xm0(4)*Wnr(1) ;
xm1(5) = xm0(5) + xm0(6)*Wnr(1) ;
xm1(6) = xm0(5) - xm0(6)*Wnr(1) ;
xm1(7) = xm0(7) + xm0(8)*Wnr(1) ;
xm1(8) = xm0(7) - xm0(8)*Wnr(1) ;
floor(xm1(:))


%% stage2
xm2(1) = xm1(1) + xm1(3)*Wnr(1) ;
xm2(3) = xm1(1) - xm1(3)*Wnr(1) ;
xm2(2) = xm1(2) + xm1(4)*Wnr(2) ;
xm2(4) = xm1(2) - xm1(4)*Wnr(2) ;
xm2(5) = xm1(5) + xm1(7)*Wnr(1) ;
xm2(7) = xm1(5) - xm1(7)*Wnr(1) ;
xm2(6) = xm1(6) + xm1(8)*Wnr(2) ;
xm2(8) = xm1(6) - xm1(8)*Wnr(2) ;
floor(xm2(:))


%% stage3
xm3(1) = xm2(1) + xm2(5)*Wnr(1) ;
xm3(5) = xm2(1) - xm2(5)*Wnr(1) ;
xm3(2) = xm2(2) + xm2(6)*Wnr(2) ;
xm3(6) = xm2(2) - xm2(6)*Wnr(2) ;
xm3(3) = xm2(3) + xm2(7)*Wnr(3) ;
xm3(7) = xm2(3) - xm2(7)*Wnr(3) ;
xm3(4) = xm2(4) + xm2(8)*Wnr(4) ;
xm3(8) = xm2(4) - xm2(8)*Wnr(4) ;
floor(xm3(:))


%% fft
fft1 = fft(x)
fft2 = fft(x2)
plot(1:8, abs(fft2))
hold on
plot(1:8, abs(xm3), 'r')