64点FFT处理器（含verilog源码）（上）

2021-07-22 10:55:25

欢迎大家关注我的微信公众号：

原文链接：64点FFT处理器（上）

前言

截止2022年2月15日，中国科学院大学《高等数字集成电路分析及设计》课程终于完结，所以我计划分享几个自己完成的实践作业，供大家交流学习。

设计收获

对FFT/IFF算法有了清晰的理解
因为本设计为结课大作业，所以我进行了比较详细的文档介绍，并在源码中增加了自动化测试脚本，方便读者快速复现。

64点FFT处理器设计报告正文

一、设计内容

设计一个FFT处理器时序逻辑电路，计算64点FFT和IFFT（N = 64）。模块整体采用流水线结构实现，能够处理连续多组输入数据。
顶层模块名为fft_64，输入输出功能定义：

名称	方向	位宽	描述
clk	I	1	系统时钟
rst_n	I	1	系统异步复位，低电平有效
inv	I	1	模式控制，0表示FFT运算，1表示IFFT运算
valid_in	I	1	输入数据有效指示，高电平有效
sop_in	I	1	每组输入数据（64个数）第一个有效数据指示，高电平有效
x_re	I	16	输入数据实部，二进制补码定点格式
x_im	I	16	输入数据虚部，二进制补码定点格式
valid_out	O	1	输出数据有效指示，高电平有效
sop_out	O	1	每组输出数据（64个数）第一个有效数据指示，高电平有效
y_re	O	16	输出数据实部，二进制补码定点格式
y_im	O	16	输出数据虚部，二进制补码定点格式

设计要求：

Verilog实现代码可综合，给出详细设计文档、综合以及仿真结果。
在每组数据处理过程中，inv信号值保持不变。
计算过程进行适当精度控制，保证输出结果精确度，输出定点格式（精度范围）可以根据需要进行调整，需要对计算结果进行误差分析。

二、FFT原理及设计

2.1 FFT算法原理

FFT（Fast Fourier Transform），快速傅立叶变换，是一种DFT（离散傅里叶变换）的高效算法。在以时频变换分析为基础的数字处理方法中，有着不可替代的作用。
设x(n)为N点有限长序列，则其DFT的运算公式为

X ( k ) = ∑ n = 1 N − 1 x ( n ) W N n k , k = 0 , 1 , . . . , N − 1 (2-1) X(k)=\sum_{n=1}^{N-1}x(n)W_N^{nk}, k=0,1,...,N-1 \tag{2-1} X(k)=n=1∑N−1x(n)WNnk,k=0,1,...,N−1(2-1)

其中： W N = e − j 2 π / N W_N=e^{-j2π/N} WN=e−j2π/N
将离散傅里叶变换公式拆分成奇偶项，则前 N/2 个点可以表示为：
X ( k ) = ∑ r = 0 N / 2 − 1 X ( 2 r ) W N 2 r k + ∑ r = 0 N / 2 − 1 X ( 2 r + 1 ) W N ( 2 r + 1 ) k = ∑ r = 0 N / 2 − 1 X ( 2 r ) W N / 2 r k + W N k ∑ r = 0 N / 2 − 1 X ( 2 r + 1 ) W N / 2 r k = A ( k ) + W N k B ( k ) , k = 0 , 1 , … , N / 2 − 1 X(k)=\sum_{r=0}^{N/2-1}X(2r)W_N^{2rk} +\sum_{r=0}^{N/2-1}X(2r+1)W_N^{(2r+1)k} \\ = \sum_{r=0}^{N/2-1}X(2r) W_{N/2}^{rk}+W_N^k \sum_{r=0}^{N/2-1}X(2r+1) W_{N/2}^{rk} \\ = A(k)+W_N^k B(k), k=0,1,…,N/2-1 X(k)=r=0∑N/2−1X(2r)WN2rk+r=0∑N/2−1X(2r+1)WN(2r+1)k=r=0∑N/2−1X(2r)WN/2rk+WNkr=0∑N/2−1X(2r+1)WN/2rk=A(k)+WNkB(k),k=0,1,…,N/2−1

同理，后N/2个点可以表示为：

X ( k + N / 2 ) = A ( k ) − W N k B ( k ) , k = 0 , 1 , … , N / 2 − 1 X(k+N/2)=A(k)-W_N^k B(k), k=0,1,…,N/2-1 X(k+N/2)=A(k)−WNkB(k),k=0,1,…,N/2−1

由此可知，后 N/2 个点的值完全可以通过计算前 N/2 个点时的中间过程值确定。对 A(k) 与 B(k) 继续进行奇偶分解，直至变成 2 点的 DFT，这样就可以避免很多的重复计算，实现了快速离散傅里叶变换（FFT）的过程。
考虑到64点数据量很大，不方便绘图解释原理，这里以8点FFT为例介绍基2FFT算法结构。

说明：图中 W ( r ) = W N r , r = 0 , 1 , 2 , 3 , 4 , 5 , 6 , 7 W(r)=W_N^r, r=0,1,2,3,4,5,6,7 W(r)=WNr,r=0,1,2,3,4,5,6,7

要介绍上面这一计算流图，首先要了解结构中的一些概念：

（1）级的概念
如图2-1中所以，每一列即为一级运算，设FFT运算点数为N，则需要被分割成 M = l o g 2 ⁡ N M=log_2⁡N M=log2⁡N级运算，每一级的标识为 m = 0 , 1 , 2 , … , M − 1 m=0,1,2,…,M-1 m=0,1,2,…,M−1
（2）蝶形单元
FFT计算结构由若干个蝶形运算单元组成，每个运算单元示意图如下：

蝶形单元的输入输入满足：
x m + 1 ( p ) = x m ( p ) + x m ( q ) ⋅ W N r x_{m+1} (p)=x_m (p)+x_m(q) \cdot W_N^r xm+1(p)=xm(p)+xm(q)⋅WNr
x m + 1 ( q ) = x m ( p ) − x m ( q ) ⋅ W N r x_{m+1} (q)=x_m (p)-x_m (q) \cdot W_N^r xm+1(q)=xm(p)−xm(q)⋅WNr

其中 q = p + 2 m q=p+2^m q=p+2m,每个蝶形运算单元运算时，进行一次乘法和两次加法，每一级中，均有N/2个蝶形单元。故完成一次FFT所需要的乘法次数为 N / 2 ⋅ l o g 2 ⁡ N N/2 \cdot log_2⁡N N/2⋅log2⁡N，所需要的加法次数为： N ⋅ l o g 2 ⁡ N N \cdot log_2⁡N N⋅log2⁡N

（3）组的概念
每一级 N/2 个蝶形单元可分为若干组，每一组有着相同的结构与 W N r W_N^r WNr因子分布。例如

m=0 时，可以分为 N/2=4 组。
m=1 时，可以分为 N/4=2 组。
m=M-1 时，此时只能分为 1 组。

（4） W N r W_N^r WNr因子分布
W 2 m + 1 r W_{2^{m+1}}^r W2m+1r因子存在于m级，其中 r = 0 , 1 , 2 , … , 2 m − 1 r=0,1,2,…,2^m-1 r=0,1,2,…,2m−1。举例说明：
在8点FFT第二级运算中，即m=1，蝶形因子运算为： W 4 0 W_4^0 W40和 W 4 1 W_4^1 W41，对应的8等分的旋转因子即为 W 8 0 W_8^0 W80和 W 8 2 W_8^2 W82，对应图2-1中的W(0)和W(2)；
在硬件实现时，旋转因子事先会计算好，所以会统一成最大等分数的旋转因子。更易理解的旋转因子分布式可以用如下关系表示：
假设N点FFT第m+1级的旋转因子为 W N k t W_N^{kt} WNkt，则其满足如下关系：
t = N / 2 m + 1 t=N/2^{m+1} t=N/2m+1
k = 0 , 1 , 2 , … , 2 m − 1 k=0,1,2,…,2^m-1 k=0,1,2,…,2m−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 ) X(000), X(001), X(010), X(011), X(100), X(101), X(110), X(111) 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 ) X(000), X(100), X(010), X(110), X(001), X(101), X(011), X(111) 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 ) X(0), X(4), X(2), X(6), X(1), X(5), X(3), X(7) X(0),X(4),X(2),X(6),X(1),X(5),X(3),X(7)
恰好对应FFT第一级输入数据的顺序。

2.2 FFT设计

（1）蝶形单元设计
端口定义如下：

module butterfly2(
  input                       clk        ,//系统时钟
  input                       rstn       ,//系统异步复位，低电平有效
  input                       en         ,//使能信号，表示输入数据有效
  input signed [15:0]         xp_real    ,//Xm(p)
  input signed [15:0]         xp_imag    ,
  input signed [15:0]         xq_real    ,//Xm(q)
  input signed [15:0]         xq_imag    ,
  input signed [15:0]         factor_real,//扩大8192 倍（左移 13 位）后的旋转因子
  input signed [15:0]         factor_imag,

  output                      valid      ,//输入数据有效执行信号
  output signed [15:0]        yp_real    ,//Xm+1(p)
  output signed [15:0]        yp_imag    ,
  output signed [15:0]        yq_real    ,//Xm+1(q)
  output signed [15:0]        yq_imag    
);

注意：蝶形单元为定点运算，输入factor_real和factor_imag为扩大 8192 倍（左移 13 位）后的旋转因子，且为了防止蝶形运算中的乘法和加法导致位宽逐级增大，每一级运算完成后，会对输出数据进行固定位宽的截位，也可去掉旋转因子倍数增大而带来的影响。
内部实现逻辑介绍:
蝶形单元内部使用流水线结构，将计算划分成了三级流水，每一级功能如下：

Xm(q) 与 factor进行复数乘法的乘法运算，Xm§ 左移13位（保证扩大相同的倍数）
Xm(q) 与 factor进行复数乘法的加法运算 a，Xm§进行寄存
进行蝶形运算的加法运算，并将得到的结果指定位宽截位，得到 X m + 1 ( p ) X_{m+1}(p) Xm+1(p)和 X m + 1 ( q ) X_{m+1}(q) Xm+1(q)

（2）输入缓存的设计
因为设计内容中要求输入为串行输入，所以需要对输入数据缓存。这里使用一个计数器来生成输入要寄存的地址信息，当valid_in有效时更新为counter_next的值，而counter_next的取值逻辑如下：

wire [5:0] counter_next = {6{sop_in }} & 6'b0 
                        | {6{~sop_in}} & (counter + 1);

注意，在写入操作时，为了保证数据与地址的一致，写输入缓冲的地址使用counter计数器的次态值，即counter_next的值。其输入地址代码如下：

// store input
reg signed [15:0] mem_xm_real[63:0];
reg signed [15:0] mem_xm_imag[63:0];
reg inv_reg ;
always @(posedge clk or negedge rst_n) begin
  if (~rst_n) begin
  end
  else if (valid_in) begin
    mem_xm_real[counter_next] <= x_re;
    mem_xm_imag[counter_next] <= {16{~inv}}&x_im | {16{inv}}&(-x_im); // inv==0 -> fft   inv==1 -> ifft
    inv_reg <= inv;
  end
end

（3）例化蝶形单元，实现fft/ifft
依据2.1 FFT算法原理可知，64点fft的实现，需要 l o g 2 ⁡ 64 = 6 log_2⁡64=6 log2⁡64=6级运算（对应for(m=0; m<=5; m=m+1)），且每一级需要64/2=32个蝶形单元（对应for (k=0; k<=31; k=k+1)）。所以例化的主体代码如下：

//butter instantiaiton
genvar m, k;
generate
   //6 stage
   for(m=0; m<=5; m=m+1) begin: stage
      for (k=0; k<=31; k=k+1) begin: unit
         butterfly2 u_butter(
            .clk         (clk                                      ) ,
            .rstn        (rst_n                                    ) ,
           这里做了省略
          );
      end
   end
endgenerate

另外fft算法中，要求输入数据索引二进制码位倒置，所以对xm_real[0]、xm_imag[0]初始化时，做如下处理：

genvar czw;
generate
  for (czw = 0; czw<=63; czw=czw+1) begin
    assign xm_real[0][czw] = mem_xm_real[index[czw]]; //输入地址二进制码翻转
    assign xm_imag[0][czw] = mem_xm_imag[index[czw]];
  end
endgenerate

其中index是一个一维向量，存放了0到63的二进制码位倒置值，这些值是事先用matlab计算出来的：

N=64;
index = bin2dec(fliplr(dec2bin(0:N-1)));

再者，第一级蝶形单元的计算要在输入缓冲满64个时才能开始，所以其使能的生成逻辑为counter==63，且需要第一辑32个蝶形单元同时有效，所以其代码如下：

wire   en_connect [223:0] ;
wire   en = counter==63;
assign {
  en_connect[ 0],en_connect[ 1]
           这里做了省略
  en_connect[30],en_connect[31] 
} = {
  en            ,en            ,
           这里做了省略
  en            ,en            
};

（4）输出缓存的设计

因为设计内容中要求输出为串行输出，所以需要对输入数据缓存。这里同样使用一个计数器来生成输出要寄存的地址信息，每个时钟将counter_o更新为counter_next_o的值，而counter_next_o的取值逻辑如下：

wire [5:0] counter_next_o = {6{sop_out_next }} & 6'b0 
                          | {6{~sop_out_next}} & (counter_o + 1);

注意，在写入操作时，为了保证数据与地址的一致，输出缓冲送到输出端口的地址使用counter_o计数器的次态值，即counter_next_o的值。其输出地址代码如下：

always @(posedge clk or negedge rst_n) begin
  if (~rst_n) begin
    counter_o <= 0;
    sop_out <= 0;
    valid_out <= 0;
    y_re <= 0;
    y_im <= 0;
  end
  else begin
    counter_o <= counter_next_o;
    y_re <= xm_real[6][counter_next_o];
    y_im <= {16{~inv_r[17]}}&xm_imag[6][counter_next_o] | {16{inv_r[17]}}&(-xm_imag[6][counter_next_o]); // inv==0 -> fft   inv==1 -> ifft
           这里做了省略
  end
end

说明

参考链接：https://www.runoob.com/w3cnote/verilog-fft.html

码农公寓

前言