阶段一：信号基础

🦋 第05课：FFT蝶形运算

快速傅里叶变换的核心——蝶形运算单元的原理与实现

1. 为什么需要FFT？

直接计算N点DFT需要N²次复数乘法，当N很大时（如N=1024），计算量达百万级。FFT利用旋转因子的周期性和对称性，将计算量降至(N/2)log₂N，效果显著：

N	DFT乘法次数	FFT乘法次数	加速比
64	4,096	192	21×
256	65,536	1,024	64×
1024	1,048,576	5,120	205×
4096	16,777,216	24,576	683×

2. 按时间抽取（DIT）的FFT算法

DIT-FFT的核心思想是将N点DFT分解为两个N/2点DFT：

X[k] = Σ x[n]·W_N^(nk) = Σ x[2r]·W_N^(2rk) + Σ x[2r+1]·W_N^((2r+1)k)
= G[k] + W_N^k · H[k]，其中G[k]为偶数点DFT，H[k]为奇数点DFT

2.1 蝶形运算

FFT的基本运算单元是蝶形运算（Butterfly），每个蝶形执行以下操作：

A' = A + W·B
B' = A - W·B

其中A、B为输入，W为旋转因子，A'、B'为输出。一个蝶形需要1次复数乘法和2次复数加法。

2.2 8点FFT信号流图

8点FFT需要log₂8=3级蝶形，每级4个蝶形，共12个蝶形：

第1级（2点DFT）  第2级（4点DFT）  第3级（8点DFT）
x[0]─┬─A0        ─┬───────┬──X[0]
x[4]─┘           │  ┌────┘──X[4]
x[2]─┬─A2        ─┤──┬─────X[2]
x[6]─┘           │  └─────X[6]
x[1]─┬─A1        ─┤──┬────X[1]
x[5]─┘           │  └────X[5]
x[3]─┬─A3        ─┤──┬────X[3]
x[7]─┘           └─────X[7]

输入：位反转顺序 {0,4,2,6,1,5,3,7}
输出：自然顺序 {0,1,2,3,4,5,6,7}

3. 旋转因子的性质

旋转因子 W_N^k = e^(-j2πk/N) 具有三大关键性质：

🔑 旋转因子三大性质

周期性：W_N^(k+N) = W_N^k
对称性：W_N^(k+N/2) = -W_N^k（简化半数计算）
可约性：W_N^(k) = W_(N/2)^(k/2)（分解的基础）

4. Verilog实现：蝶形运算单元

//=============================================
// dsp_butterfly.v
// FFT蝶形运算单元
// A' = A + W * B
// B' = A - W * B
// W = Wr + j*Wi (旋转因子)
//=============================================
module dsp_butterfly #(
    parameter DATA_WIDTH = 16,
    parameter TWIDDLE_WIDTH = 16,
    parameter OUT_WIDTH = 32
)(
    input  wire signed [DATA_WIDTH-1:0]    ar_in,     // A实部
    input  wire signed [DATA_WIDTH-1:0]    ai_in,     // A虚部
    input  wire signed [DATA_WIDTH-1:0]    br_in,     // B实部
    input  wire signed [DATA_WIDTH-1:0]    bi_in,     // B虚部
    input  wire signed [TWIDDLE_WIDTH-1:0] wr,        // 旋转因子实部
    input  wire signed [TWIDDLE_WIDTH-1:0] wi,        // 旋转因子虚部
    output wire signed [OUT_WIDTH-1:0]     ar_out,    // A'实部
    output wire signed [OUT_WIDTH-1:0]     ai_out,    // A'虚部
    output wire signed [OUT_WIDTH-1:0]     br_out,    // B'实部
    output wire signed [OUT_WIDTH-1:0]     bi_out     // B'虚部
);

    // W * B = (Wr + j*Wi)(Br + j*Bi)
    //       = (Wr*Br - Wi*Bi) + j*(Wr*Bi + Wi*Br)
    wire signed [DATA_WIDTH+TWIDDLE_WIDTH-1:0] wr_br = wr * br_in;
    wire signed [DATA_WIDTH+TWIDDLE_WIDTH-1:0] wi_bi = wi * bi_in;
    wire signed [DATA_WIDTH+TWIDDLE_WIDTH-1:0] wr_bi = wr * bi_in;
    wire signed [DATA_WIDTH+TWIDDLE_WIDTH-1:0] wi_br = wi * br_in;

    // 旋转后B的实部和虚部
    wire signed [DATA_WIDTH+TWIDDLE_WIDTH-1:0] wb_re = wr_br - wi_bi;
    wire signed [DATA_WIDTH+TWIDDLE_WIDTH-1:0] wb_im = wr_bi + wi_br;

    // A' = A + W*B
    assign ar_out = {(OUT_WIDTH-DATA_WIDTH){ar_in[DATA_WIDTH-1]}, ar_in} + wb_re;
    assign ai_out = {(OUT_WIDTH-DATA_WIDTH){ai_in[DATA_WIDTH-1]}, ai_in} + wb_im;

    // B' = A - W*B
    assign br_out = {(OUT_WIDTH-DATA_WIDTH){ar_in[DATA_WIDTH-1]}, ar_in} - wb_re;
    assign bi_out = {(OUT_WIDTH-DATA_WIDTH){ai_in[DATA_WIDTH-1]}, ai_in} - wb_im;

endmodule

5. Verilog实现：旋转因子生成器

//=============================================
// dsp_twiddle_generator.v
// FFT旋转因子生成器
// W_N^k = cos(2πk/N) - j*sin(2πk/N)
//=============================================
module dsp_twiddle_generator #(
    parameter N = 256,              // FFT点数
    parameter DATA_WIDTH = 16
)(
    input  wire                        clk,
    input  wire                        rst_n,
    input  wire                        en,
    input  wire [$clog2(N)-1:0]        k_index,    // 旋转因子索引
    output reg  signed [DATA_WIDTH-1:0] cos_out,    // cos(2πk/N)
    output reg  signed [DATA_WIDTH-1:0] sin_out     // -sin(2πk/N)
);

    // 旋转因子查找表
    reg signed [DATA_WIDTH-1:0] cos_table [0:N/2-1];  // 利用对称性只需半表
    reg signed [DATA_WIDTH-1:0] sin_table [0:N/2-1];

    initial begin
        integer i;
        for (i = 0; i < N/2; i = i + 1) begin
            cos_table[i] = $rtoi(
                (2.0**(DATA_WIDTH-1)-1) * $cos(2.0*3.14159265*i/N)
            );
            sin_table[i] = $rtoi(
                -(2.0**(DATA_WIDTH-1)-1) * $sin(2.0*3.14159265*i/N)
            );
        end
    end

    // 利用对称性：W^(k) = -W^(k+N/2)
    wire use_negate = (k_index >= N/2);
    wire [$clog2(N/2)-1:0] table_idx = use_negate ? 
        k_index[$clog2(N/2)-1:0] : k_index[$clog2(N/2)-1:0];

    always @(posedge clk) begin
        if (en) begin
            if (use_negate) begin
                cos_out <= -cos_table[table_idx];
                sin_out <= -sin_table[table_idx];
            end else begin
                cos_out <= cos_table[table_idx];
                sin_out <= sin_table[table_idx];
            end
        end
    end

endmodule

6. Verilog实现：位反转地址生成

//=============================================
// dsp_bit_reverse.v
// 位反转地址生成器
// FFT输入需要位反转顺序重排
//=============================================
module dsp_bit_reverse #(
    parameter ADDR_WIDTH = 8     // log2(N)
)(
    input  wire [ADDR_WIDTH-1:0] addr_in,
    output wire [ADDR_WIDTH-1:0] addr_out
);

    // 位反转：最高位变最低位，次高位变次低位...
    genvar i;
    generate
        for (i = 0; i < ADDR_WIDTH; i = i + 1) begin : gen_bit_rev
            assign addr_out[i] = addr_in[ADDR_WIDTH-1-i];
        end
    endgenerate

endmodule

7. Verilog实现：单级蝶形处理

//=============================================
// dsp_fft_stage.v
// FFT单级蝶形处理器
// 处理N点FFT中的某一级
//=============================================
module dsp_fft_stage #(
    parameter N = 16,
    parameter STAGE = 0,            // 当前级号（0=第一级）
    parameter DATA_WIDTH = 16,
    parameter TWIDDLE_WIDTH = 16,
    parameter OUT_WIDTH = 32
)(
    input  wire                              clk,
    input  wire                              rst_n,
    input  wire                              data_valid,
    input  wire signed [DATA_WIDTH-1:0]      data_re_in [0:N-1],
    input  wire signed [DATA_WIDTH-1:0]      data_im_in [0:N-1],
    output reg  signed [OUT_WIDTH-1:0]       data_re_out [0:N-1],
    output reg  signed [OUT_WIDTH-1:0]       data_im_out [0:N-1],
    output reg                               out_valid
);

    localparam LOG_N = $clog2(N);
    localparam BUTTERFLIES_PER_GROUP = 2**STAGE;
    localparam GROUPS = N / (2 * BUTTERFLIES_PER_GROUP);
    localparam STRIDE = BUTTERFLIES_PER_GROUP;

    integer i;
    reg [$clog2(N)-1:0] top_idx, bot_idx;
    reg [$clog2(N/2)-1:0] twiddle_idx;
    reg signed [DATA_WIDTH-1:0] a_re, a_im, b_re, b_im;
    reg signed [TWIDDLE_WIDTH-1:0] w_re, w_im;
    reg signed [OUT_WIDTH-1:0] ar_out, ai_out, br_out, bi_out;

    // 旋转因子查找表（本级所需）
    reg signed [TWIDDLE_WIDTH-1:0] stage_cos [0:N/2-1];
    reg signed [TWIDDLE_WIDTH-1:0] stage_sin [0:N/2-1];

    initial begin
        integer j;
        for (j = 0; j < N/2; j = j + 1) begin
            stage_cos[j] = $rtoi(
                (2.0**(TWIDDLE_WIDTH-1)-1) * $cos(2.0*3.14159265*j*GROUPS/N)
            );
            stage_sin[j] = $rtoi(
                -(2.0**(TWIDDLE_WIDTH-1)-1) * $sin(2.0*3.14159265*j*GROUPS/N)
            );
        end
    end

    // 处理状态机
    localparam IDLE = 0, PROCESS = 1, DONE_ST = 2;
    reg [1:0] state;
    reg [$clog2(N)-1:0] proc_cnt;

    always @(posedge clk or negedge rst_n) begin
        if (!rst_n) begin
            state <= IDLE;
            out_valid <= 1'b0;
            proc_cnt <= 0;
            for (i = 0; i < N; i = i + 1) begin
                data_re_out[i] <= 0;
                data_im_out[i] <= 0;
            end
        end else begin
            case (state)
                IDLE: begin
                    out_valid <= 1'b0;
                    if (data_valid) begin
                        state <= PROCESS;
                        proc_cnt <= 0;
                        // 先直通所有数据
                        for (i = 0; i < N; i = i + 1) begin
                            data_re_out[i] <= {(OUT_WIDTH-DATA_WIDTH){data_re_in[i][DATA_WIDTH-1]}, data_re_in[i]};
                            data_im_out[i] <= {(OUT_WIDTH-DATA_WIDTH){data_im_in[i][DATA_WIDTH-1]}, data_im_in[i]};
                        end
                    end
                end
                PROCESS: begin
                    // 计算当前蝶形对
                    top_idx = (proc_cnt / STRIDE) * 2 * STRIDE + (proc_cnt % STRIDE);
                    bot_idx = top_idx + STRIDE;
                    twiddle_idx = (proc_cnt % STRIDE) * GROUPS;

                    a_re = data_re_in[top_idx];
                    a_im = data_im_in[top_idx];
                    b_re = data_re_in[bot_idx];
                    b_im = data_im_in[bot_idx];
                    w_re = stage_cos[twiddle_idx];
                    w_im = stage_sin[twiddle_idx];

                    // 蝶形计算
                    data_re_out[top_idx] <= {(OUT_WIDTH-DATA_WIDTH){a_re[DATA_WIDTH-1]}, a_re}
                                            + (w_re * b_re - w_im * b_im);
                    data_im_out[top_idx] <= {(OUT_WIDTH-DATA_WIDTH){a_im[DATA_WIDTH-1]}, a_im}
                                            + (w_re * b_im + w_im * b_re);
                    data_re_out[bot_idx] <= {(OUT_WIDTH-DATA_WIDTH){a_re[DATA_WIDTH-1]}, a_re}
                                            - (w_re * b_re - w_im * b_im);
                    data_im_out[bot_idx] <= {(OUT_WIDTH-DATA_WIDTH){a_im[DATA_WIDTH-1]}, a_im}
                                            - (w_re * b_im + w_im * b_re);

                    if (proc_cnt == N/2 - 1) begin
                        state <= DONE_ST;
                    end else begin
                        proc_cnt <= proc_cnt + 1'b1;
                    end
                end
                DONE_ST: begin
                    out_valid <= 1'b1;
                    state <= IDLE;
                end
            endcase
        end
    end

endmodule

8. DIT与DIF对比

特性	DIT（按时间抽取）	DIF（按频率抽取）
输入顺序	位反转	自然顺序
输出顺序	自然顺序	位反转
蝶形运算	先乘旋转因子，后加减	先加减，后乘旋转因子
适用场景	数据需要预处理	数据按自然顺序输入

9. 练习

📝 练习1：手工蝶形

对4点序列 x[n]={1,2,3,4}，画出DIT-FFT的完整信号流图，标注每个蝶形的中间结果和旋转因子。

📝 练习2：位反转

对16点FFT，写出输入数据的位反转顺序排列。

📝 练习3：旋转因子优化

修改 dsp_twiddle_generator.v，利用旋转因子的对称性（W^(k+N/4)=jW^k），将查找表压缩到N/4大小。

📝 练习4：定点精度

分析蝶形运算在Q1.15格式下的量化误差积累。对于16级FFT（65536点），最坏情况下误差放大多少倍？

🏆 成就解锁：蝶形大师

✅ 理解了FFT相比DFT的复杂度优势

✅ 掌握了DIT-FFT的分解原理

✅ 实现了蝶形运算单元（Verilog）

✅ 实现了旋转因子生成器

✅ 实现了位反转地址和单级蝶形处理

11. 补充：FFT蝶形运算的进阶主题

FFT蝶形运算在进阶应用中还有以下重要主题值得深入研究：

📚 进阶研究方向

多通道实现：时分复用单个处理核心服务多个通道，资源利用率最大化。关键挑战是上下文切换的开销和系数存储管理。
动态重配置：运行时切换滤波器系数，支持多标准兼容。需要双缓冲系数存储，避免切换过程中的毛刺。
低功耗设计：时钟门控、操作数隔离、精度缩放、近似计算等技术组合使用，可将功耗降低50-80%。
容错设计：冗余计算、错误检测与纠正、 graceful degradation等技术在安全关键应用中的实现。
验证方法学：形式化验证、覆盖率驱动的随机测试、断言验证在DSP模块中的实践。

12. 与前后课程的关联

本课内容在整个DSP课程体系中处于承上启下的位置：

🔗 课程关联图

前置知识：离散信号与系统（第02课）、Z变换（第03课）是理解本课的理论基础
并行知识：DFT/FFT（第04-06课）提供频域分析工具，可验证本课设计的频率响应
后续应用：本课的设计方法将在实战项目（第25-30课）中得到综合运用
硬件实现：直接型/转置型结构（第09-10课）的选择直接影响本课的硬件架构

13. 设计经验总结

在实际工程中，FFT蝶形运算的设计经验总结如下：（1）始终从系统级需求出发，不要过早陷入实现细节；（2）先用浮点仿真验证算法正确性，再转为定点；（3）系数量化后必须重新验证频率响应；（4）硬件实现前用C/MATLAB模型作为参考；（5）综合后检查时序和资源是否符合预期；（6）板级验证时用已知的测试信号对比参考输出。每一步验证都是必要的，跳过任何一步都可能导致最终的调试时间成倍增加。

📌 关键设计检查清单

□ 算法仿真验证通过
□ 定点量化分析完成
□ Verilog代码通过lint检查
□ 仿真波形与参考模型一致
□ 综合时序满足要求
□ 资源使用在预算内
□ 板级测试通过

← 上一课：DFT 📚 目录下一课：FFT完整引擎 →