Verilog FFT设计
fft（fast fourier transform），快速傅立叶变换，是一种 dft（离散傅里叶变换）的高效算法。 在以时频变换分析为基础的数字处理方法中，有着不可替代的作用。
fft 原理
公式推导
dft 的运算公式为：
其中，
将离散傅里叶变换公式拆分成奇偶项，则前 n/2 个点可以表示为：
同理，后 n/2 个点可以表示为：
由此可知，后 n/2 个点的值完全可以通过计算前 n/2 个点时的中间过程值确定。 对 a[k] 与 b[k] 继续进行奇偶分解，直至变成 2 点的 dft，这样就可以避免很多的重复计算，实现了快速离散傅里叶变换（fft）的过程。
算法结构
8 点 fft 计算的结构示意图如下。
由图可知，只需要简单的计算几次乘法和加法，便可完成离散傅里叶变换过程，而不是对每个数据进行繁琐的相乘和累加。
重要特性
(1) 级的概念
每分割一次，称为一级运算。
设 fft 运算点数为 n，共有 m 级运算，则它们满足：
每一级运算的标识为 m = 0, 1, 2, ..., m-1。
为了便于分割计算，fft 点数 n 的取值经常为 2 的整数次幂。
(2) 蝶形单元
fft 计算结构由若干个蝶形运算单元组成，每个运算单元示意图如下：
蝶形单元的输入输出满足：
其中，。
每一个蝶形单元运算时，进行了一次乘法和两次加法。
每一级中，均有 n/2 个蝶形单元。
故完成一次 fft 所需要的乘法次数和加法次数分别为：
(3) 组的概念
每一级 n/2 个蝶形单元可分为若干组，每一组有着相同的结构与因子分布。
例如 m=0 时，可以分为 n/2=4 组。
m=1 时，可以分为 n/4=2 组。
m=m-1 时，此时只能分为 1 组。
(4) 因子分布
因子存在于 m 级，其中 。
在 8 点 fft 第二级运算中，即 m=1 ，蝶形运算因子可以化简为：
(5) 码位倒置
对于 n=8 点的 fft 计算，x(0) ~ x(7) 位置对应的 2 进制码为：
x（000）， x（001）， x（010）， x（011）， x（100）， x（101）， x（110）， x（111）
将其位置的 2 进制码进行翻转：
x（000）， x（100）， x（010）， x（110）， x（001）， x（101）， x（011）， x（111）
此时位置对应的 10 进制为：
x（0）， x（4）， x（2）， x（6）， x（1）， x（5）， x（3）， x（7）
恰好对应 fft 第一级输入数据的顺序。
该特性有利于 fft 的编程实现。
fft 设计
设计说明
为了利用仿真简单的说明 fft 的变换过程，数据点数取较小的值 8。
如果数据是串行输入，需要先进行缓存，所以设计时数据输入方式为并行。
数据输入分为实部和虚部共 2 部分，所以计算结果也分为实部和虚部。
设计采用流水结构，暂不考虑资源消耗的问题。
为了使设计结构更加简单，这里做一步妥协，乘法计算直接使用乘号。 如果 fft 设计应用于实际，一定要将乘法结构换成可以流水的乘法器，或使用官方提供的效率较高的乘法器 ip。
蝶形单元设计
蝶形单元为定点运算，需要对旋转因子进行定点量化。
借助 matlab 将旋转因子扩大 8192 倍（左移 13 位），可参考附录。
为了防止蝶形运算中的乘法和加法导致位宽逐级增大，每一级运算完成后，要对输出数据进行固定位宽的截位，也可去掉旋转因子倍数增大而带来的影响。
代码如下：
`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; //提前设置好位宽余量，防止数据溢出 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顶层例化
根据 fft 算法结构示意图，将蝶形单元例化，完成最后的 fft 功能。
可根据每一级蝶形单元的输入输出对应关系，依次手动例化 12 次，也可利用 generate 进行例化，此时就需要非常熟悉 fft 中“组”和“级”的特点：
(1) 8 点 fft 设计，需要 3 级运算，每一级有 4 个蝶形单元，每一级的组数目分别是 4、2、1。
(2) 每一级的组内一个蝶形单元中两个输入端口的距离恒为 （m 为级标号，对应左移运算 1<(3) 每一级相邻组间的第一个蝶形单元的第一个输入端口的距离为 （对应左移运算 2<例化代码如下。
其中，矩阵信号 xm_real（xm_imag）的一维、二维地址是代表级和组的标识。
在判断信号端口之间的连接关系时，使用了看似复杂的判断逻辑，而且还带有乘号，其实最终生成的电路和手动编写代码例化 12 个蝶形单元的方式是完全相同的。 因为 generate 中的变量只是辅助生成实际的电路，相关值的计算判断都已经在编译时完成。 这些变量更不会生成实际的电路，只是为更快速的模块例化提供了一种方法。
timescale 1ns/100psmodule 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 //输入初始化，和码位有关倒置 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] ) , //是否再组内？组编号+组内编号：下组编号+新组内编号 .xp_real (xm_real[ m ] [k[m:0] < (1< 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仿真结果
大致可以看出，fft 结果可以流水输出。
用 matlab 自带的 fft 函数对相同数据进行运算，前 2 组数据 fft 结果如下。
可以看出，第一次输入的数据信号只有实部有效时，fft 结果是完全一样的。
但是第二次输入的数据复部也有信号，此时两者之间的结果开始有误差，有时误差还很大。
用 matlab 对 verilog 实现的 fft 过程进行模拟，发现此过程的 fft 结果和 verilog 实现的 fft 结果基本一致。
将有误差的两种 fft 结果取绝对值进行比较，图示如下。
可以看出，fft 结果的趋势大致相同，但在个别点有肉眼可见的误差。
设计总结：
就如设计蝶形单元时所说，旋转因子量化时，位宽的选择就会引入误差。
而且每个蝶形单元的运算结果都会进行截取，也会引入误差。
matlab 计算 fft 时不用考虑精度问题，以其最高精度对数据进行 fft 计算。
以上所述，都会导致最后两种 fft 计算方式结果的差异。
感兴趣的学者，可以将旋转因子和输入数据位宽再进行一定的增加，fft 点数也可以增加，然后再进行仿真对比，相对误差应该会减小。
附录：matlab 使用
生成旋转因子
8 点 fft 只需要用到 4 个旋转因子。 旋转因子扩大倍数为 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 ; %负数的补码 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) %实部 wnr(2*r+2,:) = dec2hex(wnr_imag) %虚部endfft 结果对比
matlab 模拟 verilog 实现 fft 的过程如下，也包括 2 种 fft 结果的对比。
clear all;close all;clc;%=======================================================% 8dots fft%=======================================================for r=0:3 wnr(r+1) = cos(pi/4*r) - j*sin(pi/4*r) ;endx = [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(:))%% stage2xm2(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(:))%% stage3xm3(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(:))%% fftfft1 = fft(x)fft2 = fft(x2)plot(1:8, abs(fft2))hold onplot(1:8, abs(xm3), 'r')
虚拟现实技术对于电子竞技产业有着什么影响
小米成绩斐然 已领衔印度智能手机市场
dfrobotSIM7000C扩展板简介
土壤含水量测定仪是如何测定土壤的含水量
工控机触摸屏一体
Verilog FFT设计
基于TDA2009的25瓦功率放大器电路
复旦大学利用柔性薄膜组装集成芯片传感器，实现多环境参数探测功能
Leader千兆六电口｜为机器视觉而生的网卡
5G如何让无人驾驶跑起来
采用新技术，短时间现场修复造纸浆泵泵体冲刷磨损
TCL Xess mini平板：多场景化、家庭娱乐“宛如生活”
亿铸“芯”势力赋能AI未来
工程师电子制作故事：锂电池充电器DIY设计
sgnx函数的概念、特征和用途
华云数据许广彬：这些企业的数字化转型故事，印证着超融合的魅力
甚薄型QFN封装技术
立得空间无人机快速测绘系统的系统架构分析
苹果正式发布了iOS 12.3的第四个测试版用户可以通过iOS的OTA更新
ITE DP转HDMI芯片IT6561
fft 原理
公式推导
dft 的运算公式为：
其中，
将离散傅里叶变换公式拆分成奇偶项，则前 n/2 个点可以表示为：
同理，后 n/2 个点可以表示为：
由此可知，后 n/2 个点的值完全可以通过计算前 n/2 个点时的中间过程值确定。 对 a[k] 与 b[k] 继续进行奇偶分解，直至变成 2 点的 dft，这样就可以避免很多的重复计算，实现了快速离散傅里叶变换（fft）的过程。
算法结构
8 点 fft 计算的结构示意图如下。
由图可知，只需要简单的计算几次乘法和加法，便可完成离散傅里叶变换过程，而不是对每个数据进行繁琐的相乘和累加。
重要特性
(1) 级的概念
每分割一次，称为一级运算。
设 fft 运算点数为 n，共有 m 级运算，则它们满足：
每一级运算的标识为 m = 0, 1, 2, ..., m-1。
为了便于分割计算，fft 点数 n 的取值经常为 2 的整数次幂。
(2) 蝶形单元
fft 计算结构由若干个蝶形运算单元组成，每个运算单元示意图如下：
蝶形单元的输入输出满足：
其中，。
每一个蝶形单元运算时，进行了一次乘法和两次加法。
每一级中，均有 n/2 个蝶形单元。
故完成一次 fft 所需要的乘法次数和加法次数分别为：
(3) 组的概念
每一级 n/2 个蝶形单元可分为若干组，每一组有着相同的结构与因子分布。
例如 m=0 时，可以分为 n/2=4 组。
m=1 时，可以分为 n/4=2 组。
m=m-1 时，此时只能分为 1 组。
(4) 因子分布
因子存在于 m 级，其中 。
在 8 点 fft 第二级运算中，即 m=1 ，蝶形运算因子可以化简为：
(5) 码位倒置
对于 n=8 点的 fft 计算，x(0) ~ x(7) 位置对应的 2 进制码为：
x（000）， x（001）， x（010）， x（011）， x（100）， x（101）， x（110）， x（111）
将其位置的 2 进制码进行翻转：
x（000）， x（100）， x（010）， x（110）， x（001）， x（101）， x（011）， x（111）
此时位置对应的 10 进制为：
x（0）， x（4）， x（2）， x（6）， x（1）， x（5）， x（3）， x（7）
恰好对应 fft 第一级输入数据的顺序。
该特性有利于 fft 的编程实现。
fft 设计
设计说明
为了利用仿真简单的说明 fft 的变换过程，数据点数取较小的值 8。
如果数据是串行输入，需要先进行缓存，所以设计时数据输入方式为并行。
数据输入分为实部和虚部共 2 部分，所以计算结果也分为实部和虚部。
设计采用流水结构，暂不考虑资源消耗的问题。
为了使设计结构更加简单，这里做一步妥协，乘法计算直接使用乘号。 如果 fft 设计应用于实际，一定要将乘法结构换成可以流水的乘法器，或使用官方提供的效率较高的乘法器 ip。
蝶形单元设计
蝶形单元为定点运算，需要对旋转因子进行定点量化。
借助 matlab 将旋转因子扩大 8192 倍（左移 13 位），可参考附录。
为了防止蝶形运算中的乘法和加法导致位宽逐级增大，每一级运算完成后，要对输出数据进行固定位宽的截位，也可去掉旋转因子倍数增大而带来的影响。
代码如下：
`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; //提前设置好位宽余量，防止数据溢出 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顶层例化
根据 fft 算法结构示意图，将蝶形单元例化，完成最后的 fft 功能。
可根据每一级蝶形单元的输入输出对应关系，依次手动例化 12 次，也可利用 generate 进行例化，此时就需要非常熟悉 fft 中“组”和“级”的特点：
(1) 8 点 fft 设计，需要 3 级运算，每一级有 4 个蝶形单元，每一级的组数目分别是 4、2、1。
(2) 每一级的组内一个蝶形单元中两个输入端口的距离恒为 （m 为级标号，对应左移运算 1<
其中，矩阵信号 xm_real（xm_imag）的一维、二维地址是代表级和组的标识。
在判断信号端口之间的连接关系时，使用了看似复杂的判断逻辑，而且还带有乘号，其实最终生成的电路和手动编写代码例化 12 个蝶形单元的方式是完全相同的。 因为 generate 中的变量只是辅助生成实际的电路，相关值的计算判断都已经在编译时完成。 这些变量更不会生成实际的电路，只是为更快速的模块例化提供了一种方法。
timescale 1ns/100psmodule 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 //输入初始化，和码位有关倒置 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] ) , //是否再组内？组编号+组内编号：下组编号+新组内编号 .xp_real (xm_real[ m ] [k[m:0] < (1<
大致可以看出，fft 结果可以流水输出。
用 matlab 自带的 fft 函数对相同数据进行运算，前 2 组数据 fft 结果如下。
可以看出，第一次输入的数据信号只有实部有效时，fft 结果是完全一样的。
但是第二次输入的数据复部也有信号，此时两者之间的结果开始有误差，有时误差还很大。
用 matlab 对 verilog 实现的 fft 过程进行模拟，发现此过程的 fft 结果和 verilog 实现的 fft 结果基本一致。
将有误差的两种 fft 结果取绝对值进行比较，图示如下。
可以看出，fft 结果的趋势大致相同，但在个别点有肉眼可见的误差。
设计总结：
就如设计蝶形单元时所说，旋转因子量化时，位宽的选择就会引入误差。
而且每个蝶形单元的运算结果都会进行截取，也会引入误差。
matlab 计算 fft 时不用考虑精度问题，以其最高精度对数据进行 fft 计算。
以上所述，都会导致最后两种 fft 计算方式结果的差异。
感兴趣的学者，可以将旋转因子和输入数据位宽再进行一定的增加，fft 点数也可以增加，然后再进行仿真对比，相对误差应该会减小。
附录：matlab 使用
生成旋转因子
8 点 fft 只需要用到 4 个旋转因子。 旋转因子扩大倍数为 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 ; %负数的补码 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) %实部 wnr(2*r+2,:) = dec2hex(wnr_imag) %虚部endfft 结果对比
matlab 模拟 verilog 实现 fft 的过程如下，也包括 2 种 fft 结果的对比。
clear all;close all;clc;%=======================================================% 8dots fft%=======================================================for r=0:3 wnr(r+1) = cos(pi/4*r) - j*sin(pi/4*r) ;endx = [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(:))%% stage2xm2(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(:))%% stage3xm3(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(:))%% fftfft1 = fft(x)fft2 = fft(x2)plot(1:8, abs(fft2))hold onplot(1:8, abs(xm3), 'r')
虚拟现实技术对于电子竞技产业有着什么影响
小米成绩斐然 已领衔印度智能手机市场
dfrobotSIM7000C扩展板简介
土壤含水量测定仪是如何测定土壤的含水量
工控机触摸屏一体
Verilog FFT设计
基于TDA2009的25瓦功率放大器电路
复旦大学利用柔性薄膜组装集成芯片传感器，实现多环境参数探测功能
Leader千兆六电口｜为机器视觉而生的网卡
5G如何让无人驾驶跑起来
采用新技术，短时间现场修复造纸浆泵泵体冲刷磨损
TCL Xess mini平板：多场景化、家庭娱乐“宛如生活”
亿铸“芯”势力赋能AI未来
工程师电子制作故事：锂电池充电器DIY设计
sgnx函数的概念、特征和用途
华云数据许广彬：这些企业的数字化转型故事，印证着超融合的魅力
甚薄型QFN封装技术
立得空间无人机快速测绘系统的系统架构分析
苹果正式发布了iOS 12.3的第四个测试版用户可以通过iOS的OTA更新
ITE DP转HDMI芯片IT6561