FPGA实现8点FFT

前面我们讲了FFT的原理以及其在C++上的实现，可以参考我的博客：

快速傅里叶变换学习（超详细，附代码实现）_Patarw_Li的博客-CSDN博客

C++实现FFT算法（迭代版本）_Patarw_Li的博客-CSDN博客

下面我们会在FPGA上用Verilog实现8点FFT，下面是需要注意的几点：

1. 旋转因子

在FPGA中直接计算旋转因子是一件比较麻烦的事，因此我们使用MATLAB将旋转因子计算好后直接在verilog中赋值即可：

 生成旋转因子实部和虚部的matlab代码：

%fft旋转因子生成表
%w代表返回值,n代表运算点数
%这里将w放大，是因为浮点运算比较消耗时间，因此将其化为整数
clear all;
clc;

n = 8;
for i=0:n/2
    wr=cos(-2*pi*i/n)*256 %实部
    wi=sin(-2*pi*i/n)*256 %虚部
end

之所以要乘以256是因为浮点运算比较消耗时间，因此将其化为整数。

在verilog中的赋值代码：

wire signed [23:0] wndatareal[0:3];    //旋转因子实部数据
wire signed [23:0] wndataimg[0:3];     //旋转因子虚部数据
//通过matlab计算出旋转因子w0、w1、w2、w3，再乘以256，然后直接赋值
assign wndatareal[0] = 24'd256;
assign wndataimg[0] = 24'd0;
assign wndatareal[1] = 24'd181;
assign wndataimg[1] = -24'd181;
assign wndatareal[2] = 24'd0;
assign wndataimg[2] = -24'd256;
assign wndatareal[3] = -24'd181;
assign wndataimg[3] = -24'd181;

2. 码位倒置

代码实现：

5'd1:begin  //码位倒置
            dft_oridata[data_cnt] <= input_data[{data_cnt[0],data_cnt[1],data_cnt[2]}];
            
            data_cnt <= data_cnt+1;
            if(data_cnt == N-1) begin
                state <= state+1;
                data_cnt <= 0;
            end
end

3. 蝶形运算

设z1=a+bi，z2=c+di(a、b、c、d∈R)是任意两个复数，那么它们的积(a+bi)(c+di)=(ac-bd)+(bc+ad)i。 

代码实现：

//(a+bi)*(c+di)=a*c-b*d + (a*d+b*c)i
cache_real = dft_firoutreal[data_cnt+cal_stage]*wndatareal[wndata_cnt]-
             dft_firoutimg[data_cnt+cal_stage]*wndataimg[wndata_cnt];
cache_img = dft_firoutreal[data_cnt+cal_stage]*wndataimg[wndata_cnt]+
            dft_firoutimg[data_cnt+cal_stage]*wndatareal[wndata_cnt];
            
cache_realres[31:0] = (dft_firoutreal[data_cnt]<<8) + cache_real;
cache_imgres[31:0] = (dft_firoutimg[data_cnt]<<8) + cache_img;
dft_secoutreal[data_cnt] = cache_realres[31:8];
dft_secoutimg[data_cnt] = cache_imgres[31:8];
            
cache_realres[31:0] = (dft_firoutreal[data_cnt]<<8) - cache_real;
cache_imgres[31:0] = (dft_firoutimg[data_cnt]<<8) - cache_img;
dft_secoutreal[data_cnt+cal_stage] = cache_realres[31:8];
dft_secoutimg[data_cnt+cal_stage] = cache_imgres[31:8];

代码的上面4行是乘以旋转因子，下面8行是两个序列相加减，之所以dft_firoutreal要往左移8位是因为旋转因子被放大了256倍，后面在加减完成后，再缩小256倍。

总体Verilog代码（包含testbench）：

//2023/4/18 lzp 8点fft
`timescale 1ns/1ps


module fft
(
    input clk_100m,             //时钟
    input rst_n,                //复位
    input signed [23:0] data,   //输入数据
    output [23:0] data_r,  //输出数据的实部
    output [23:0] data_i   //输出数据的虚部
);

//常数定义
parameter N = 4'd8;

//数组和标志位定义
wire signed [23:0] wndatareal[0:3];    //旋转因子实部数据
wire signed [23:0] wndataimg[0:3];     //旋转因子虚部数据
//通过matlab计算出旋转因子w0、w1、w2、w3，再乘以256，然后直接赋值
assign wndatareal[0] = 24'd256;
assign wndataimg[0] = 24'd0;
assign wndatareal[1] = 24'd181;
assign wndataimg[1] = -24'd181;
assign wndatareal[2] = 24'd0;
assign wndataimg[2] = -24'd256;
assign wndatareal[3] = -24'd181;
assign wndataimg[3] = -24'd181;

reg signed [23:0] input_data [0:N-1];    //原始输入数据,最高位为符号位
reg signed [23:0] dft_oridata [0:N-1];   //码位倒置后的输入数据，最高位为符号位

//第一级输出
reg signed [23:0] dft_firoutreal [0:N-1];  //第一级DFT输出数据实部，最高位为符号位
reg signed [23:0] dft_firoutimg [0:N-1];   //第一级DFT输出数据虚部，最高位为符号位
//第二级输出
reg signed [23:0] dft_secoutreal [0:N-1];  //第二级DFT输出数据实部，最高位为符号位
reg signed [23:0] dft_secoutimg [0:N-1];   //第二级DFT输出数据虚部，最高位为符号位
//第三级输出
reg signed [23:0] dft_trdoutreal [0:N-1];  //第三级DFT输出数据实部，最高位为符号位
reg signed [23:0] dft_trdoutimg [0:N-1];   //第三级DFT输出数据虚部，最高位为符号位

//fft流水线
reg [4:0]   state;        //状态机
reg [7:0]   fft_stage;    //fft运算阶数
reg [3:0]  cal_stage;    //每组做几次蝶形运算
reg [3:0]   data_cnt;     //数据计数
reg [3:0]  wndata_cnt;   //旋转因子计数
reg [3:0]  group_cnt;    //每一阶计算时都会数据进行不同的分组，该寄存器计数已计算的分组

reg [31:0]  cache_real;   //数据实部缓存（用于计算）
reg [31:0]  cache_img;    //数据虚部缓存（用于计算）
reg [31:0]  cache_realres;//实部计算结果
reg [31:0]  cache_imgres; //虚部计算结果

reg [23:0] fft_data_r;
reg [23:0] fft_data_i;

always@(posedge clk_100m or negedge rst_n)
if(!rst_n) begin
    state <= 0;
    fft_stage <= 1;
    cal_stage <= 1;
    data_cnt <= 0;
    wndata_cnt <= 0;
    group_cnt <= 0;
    
    cache_real <= 0;
    cache_img <= 0;
    cache_realres <= 0;
    cache_imgres <= 0;
    
    fft_data_r <= 0;
    fft_data_i <= 0;
end
else begin
    case(state)
        5'd0:begin  //将输入的data填入input_data中
            input_data[data_cnt] <= data;
            
            data_cnt <= data_cnt+1;
            if(data_cnt == N-1) begin
                state <= state+1;
                data_cnt <= 0;
            end
        end
        5'd1:begin  //码位倒置
            dft_oridata[data_cnt] <= input_data[{data_cnt[0],data_cnt[1],data_cnt[2]}];
            
            data_cnt <= data_cnt+1;
            if(data_cnt == N-1) begin
                state <= state+1;
                data_cnt <= 0;
            end
        end
        5'd2:begin  //第一级蝶形运算，N/2点DFT
            cache_real = dft_oridata[data_cnt+cal_stage]*wndatareal[wndata_cnt];
            cache_img = dft_oridata[data_cnt+cal_stage]*wndataimg[wndata_cnt];
            
            cache_realres[31:0] = (dft_oridata[data_cnt]<<8) + cache_real;
            cache_imgres[31:0] = cache_img;
            dft_firoutreal[data_cnt] = cache_realres[31:8];
            dft_firoutimg[data_cnt] = cache_imgres[31:8];
            
            cache_realres[31:0] = (dft_oridata[data_cnt]<<8) - cache_real;
            cache_imgres[31:0] = 0-cache_img;
            dft_firoutreal[data_cnt+cal_stage] = cache_realres[31:8];
            dft_firoutimg[data_cnt+cal_stage] = cache_imgres[31:8];
            
            data_cnt = data_cnt+1;
            wndata_cnt = wndata_cnt+1;
            if(wndata_cnt == cal_stage) begin //代表这组计算完毕，开始下一组计算
                group_cnt = group_cnt+1;
                data_cnt = data_cnt+cal_stage; //将data_cnt指向下一组的起始位置
                wndata_cnt = 0;
                if(group_cnt == (N>>fft_stage)) begin //代表所有组都计算完毕，开始下一级蝶形运算
                    state <= state+1;
                    group_cnt <= 0;
                    data_cnt <= 0;
                    fft_stage <= fft_stage+1;
                    cal_stage <= cal_stage<<1;
                end
            end
        end
        5'd3:begin  //第二级蝶形运算，N/4点DFT
            //(a+bi)*(c+di)=a*c-b*d + (a*d+b*c)i
            cache_real = dft_firoutreal[data_cnt+cal_stage]*wndatareal[wndata_cnt]-
                         dft_firoutimg[data_cnt+cal_stage]*wndataimg[wndata_cnt];
            cache_img = dft_firoutreal[data_cnt+cal_stage]*wndataimg[wndata_cnt]+
                        dft_firoutimg[data_cnt+cal_stage]*wndatareal[wndata_cnt];
            
            cache_realres[31:0] = (dft_firoutreal[data_cnt]<<8) + cache_real;
            cache_imgres[31:0] = (dft_firoutimg[data_cnt]<<8) + cache_img;
            dft_secoutreal[data_cnt] = cache_realres[31:8];
            dft_secoutimg[data_cnt] = cache_imgres[31:8];
            
            cache_realres[31:0] = (dft_firoutreal[data_cnt]<<8) - cache_real;
            cache_imgres[31:0] = (dft_firoutimg[data_cnt]<<8) - cache_img;
            dft_secoutreal[data_cnt+cal_stage] = cache_realres[31:8];
            dft_secoutimg[data_cnt+cal_stage] = cache_imgres[31:8];
            
            data_cnt = data_cnt+1;
            wndata_cnt = wndata_cnt+1;
            if(wndata_cnt == cal_stage) begin //代表这组计算完毕，开始下一组计算
                group_cnt = group_cnt+1;
                data_cnt = data_cnt+cal_stage;//将data_cnt指向下一组的起始位置
                wndata_cnt = 0;
                if(group_cnt == (N>>fft_stage)) begin //代表所有组都计算完毕，开始下一级蝶形运算
                    state <= state+1;
                    group_cnt <= 0;
                    data_cnt <= 0;
                    fft_stage <= fft_stage+1;
                    cal_stage <= cal_stage<<1;
                end
            end
        end
        5'd4:begin  //第三级蝶形运算，N/8点DFT
            //(a+bi)*(c+di)=a*c-b*d + (a*d+b*c)i
            cache_real = dft_secoutreal[data_cnt+cal_stage]*wndatareal[wndata_cnt]-
                         dft_secoutimg[data_cnt+cal_stage]*wndataimg[wndata_cnt];
            cache_img = dft_secoutreal[data_cnt+cal_stage]*wndataimg[wndata_cnt]+
                        dft_secoutimg[data_cnt+cal_stage]*wndatareal[wndata_cnt];
            
            cache_realres[31:0] = (dft_secoutreal[data_cnt]<<8) + cache_real;
            cache_imgres[31:0] = (dft_secoutimg[data_cnt]<<8) + cache_img;
            dft_trdoutreal[data_cnt] = cache_realres[31:8];
            dft_trdoutimg[data_cnt] = cache_imgres[31:8];
            
            cache_realres[31:0] = (dft_secoutreal[data_cnt]<<8) - cache_real;
            cache_imgres[31:0] = (dft_secoutimg[data_cnt]<<8) - cache_img;
            dft_trdoutreal[data_cnt+cal_stage] = cache_realres[31:8];
            dft_trdoutimg[data_cnt+cal_stage] = cache_imgres[31:8];
            
            data_cnt = data_cnt+1;
            wndata_cnt = wndata_cnt+1;
            if(wndata_cnt == cal_stage) begin //代表这组计算完毕，开始下一组计算
                group_cnt = group_cnt+1;
                data_cnt = data_cnt+cal_stage;//将data_cnt指向下一组的起始位置
                wndata_cnt = 0;
                if(group_cnt == (N>>fft_stage)) begin //最后一级蝶形运算
                    state <= state+1;
                    group_cnt <= 0;
                    data_cnt <= 0;
                    fft_stage <= fft_stage+1;
                    cal_stage <= cal_stage<<1;
                end
            end
        end
        5'd5:begin  
            fft_data_r <= dft_trdoutreal[data_cnt];
            fft_data_i <= dft_trdoutimg[data_cnt];
            
            data_cnt <= data_cnt+1; 
            if(data_cnt == N-1) begin
                state <= state+1; 
                data_cnt <= 0;
            end
        end
        default:begin 
            cal_stage <= 1;
            fft_stage <= 1; 
            data_cnt <= 0; 
        end
    endcase
end

//连接到模块输出
assign data_r = fft_data_r;
assign data_i = fft_data_i;

endmodule

//testbench
module fft_tb;

parameter N = 4'd8;

reg clk_100m;
reg rst_n;
reg signed [23:0] data;
wire [23:0] data_r;
wire [23:0] data_i;

reg [23:0] input_data [0:N-1];
reg [3:0]   data_cnt;     //数据计数

initial begin
    clk_100m=0;
    rst_n<=0;
    data<=0;
    data_cnt<=0;
    //用于输入的数据1 0 4 3 1 0 4 3
    input_data[0]<=1;
    input_data[1]<=0;
    input_data[2]<=4;
    input_data[3]<=3;
    input_data[4]<=1;
    input_data[5]<=0;
    input_data[6]<=4;
    input_data[7]<=3;
    
    #15 rst_n<=1;
end

//每隔一个时钟周期填入一个数据
always#10 begin
    data <= input_data[data_cnt];
    data_cnt <= data_cnt+1;
    if(data_cnt == N-1) begin
        data_cnt <= 0;
    end
end

//同名例化
fft fft1(
    clk_100m,
    rst_n,
    data,
    data_r,
    data_i
);

always#5 clk_100m<=~clk_100m;

endmodule

测试

输入：

matlab中fft输出结果：

fpga仿真：

将如下三个参数加入观察：input_data是输入数据，dft_trdoutreal和dft_trdoutimg对应输出的实部和虚部：

仿真波形图：

可以看到和matlab输出结果一致。

参考材料

1024点fft原理及fpga实现_雷凌峻毅的博客-CSDN博客