前面我们讲了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输出结果一致。文章来源:https://www.toymoban.com/news/detail-742259.html
参考材料
1024点fft原理及fpga实现_雷凌峻毅的博客-CSDN博客文章来源地址https://www.toymoban.com/news/detail-742259.html
到了这里,关于FPGA实现8点FFT的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!