vivado FFT IP Core

前言

  由于计算资源受限，准备将上位机 FFT 的计算移动到 FPGA 侧，通过 FFT IP 进行并行计算加速。本文对 Xilinx FFT IP Core 进行介绍。

FFT IP 接口介绍

接口简介

  Fast Fourier Transform IP 核，接口如下

采用 AXI4-Stream 协议（与标准的 AXI4 协议略有不同），包含 S_AXI_CONFIG，S_AXI_DATA，M_AXI_DATA，M_AXI_STATUS 四个通道，每个通道均采用 valid/ready 握手机制。

• S_AXI_CONFIG

  配置通道，可以对工作模式（FFT/IFFT）、缩放因子、变换点数等进行配置。vaild、ready 信号不再赘述，下同；tdata 携带循环前缀长度 CP_LEN、工作模式 FWD/INV、变换点数 NFFT 和缩放因子 SCALE_SCH 等信息。

• S_AXI_DATA

  数据输入通道，tdata 携带待变换数据的实部 XN_RE 和虚部 XN_IM，tlast 指示是否是最后一个数据。

• M_AXI_DATA

  变换结果输出通道，tdata 输出结果的实部 XK_RE 和虚部 XK_IM，tlast 指示是否是最后一个数据，tuser 携带额外的样本信息，如输出数据索引 XK_INDEX、算术溢出指示器 OVFLO 和块指数 BLK_EXP。

• M_AXI_STATUS

  状态通道，tdata 传递块指数 BLK_EXP 或算术溢出指示器 OVFLO 两个状态。

• aclk, aresetn, aclken

  略；若 aclken 功能被启用，会降低 aclk 的最大运行频率；arestn 需保持至少 2 clk_cycle 以确保正确复位。

• event
  • event_frame_started，O，指示开始处理新帧的断言；
  • event_tlast_unexpected，O，当在非最后一个数据时看到 tlast 为 High 的断言；
  • event_tlast_missing，O，在随后一个数据看到 tlast 为 Low 的断言；
  • event_fft_overflow，O，数据输出通道出现数据溢出的断言；
  • event_data_in_channel_halt，O，数据输入通道需要数据而没有可用数据时的断言；
  • event_data_out_channel_halt，O，数据输出通道准备输出数据而无法进行时的断言；
  • event_status_channel_halt，O，状态输出通道准备写入状态信息而无法进行时的断言。

tdata 格式说明

• 配置通道 tdata

  s_axis_config_tdata，具有 FWD/INV、CP_LEN、NFFT、SCALE_SCH 等几个字段，8bit 对齐（所有的 tdata 和 tuser 都是 8bit 对齐的），组织格式如下（虚线字段为可选字段，下同；若在 IP 配置时没有选用这些可选字段，则在 tdata/tuser 中相应字段不存在）

  FWD/INV 配置 IP 核的工作模式，为 1 为正常的 FFT，为 0 表示进行 IFFT，每个 bit 对应一个 FFT 通道，LSB（bit0） 代表通道 0，可逐帧设置。

  NFFT，位宽 5 bits，设置运行时的 FFT 点数（需小于最大点数，最大点数在生成 IP 时设置，最大 65536（对应 5’h10）），值为 ${\log }_2(NFFT)$，比如要设置为 1024 点，则值应设置为 5‘d10。

  CP_LEN，位宽为 ${\log }_2(maxmumpointsize)$，循环前缀长度（Cyclic Prefix Length）：在输出整个变换结果之前，最初作为循环前缀输出的变换末端的样本数。可以设置为 0 到小于 NFFT 的任意数。

  SCALE_SCH，设置输出结果的缩放比例，可逐帧设置。对于 Pipelined Streaming I/O 和 Radix-4 burst I/O 架构，位宽为 $2\times ceil(\frac{NFFT}{2})$；对于 Radix-2 Burst I/O 和 Radix-2 Lite Burst I/O 架构，位宽为 $2\times NFFT$；这里的 NFFT 指 ${\log }_2(maxmumpointsize)$，即 FFT/IFFT 的蝶形算法阶段 (stage) 数，Radix-2 的每个蝶形单元拾取 2 个复数，而 Radix-4 的蝶形单元会拾取 4 个复数。每个 FFT 通道对应一个 SCALE_SCH，通道间的 SCALE_SCH 无 Padding，而整个字段需要向 8bit 对齐。对于 Burst I/O 架构，每个 stage 对应 2bit 的放缩系数，3/2/1/0，对应右移移位数；例如，1024 点 Radix-4，可设置为 [1 0 2 3 2]（LSB 对应第一个 stage），表示第一个 stage 移动 2 位，第二个 stage 移动 3 位，…，一共移动了 8 位，对应缩放系数 $1/256$；又如，128 点 Radix-2，可设置为 [1 1 1 1 0 1 2]，对应总缩放系数 $1/128$。 对于 Pipelined Streaming I/O 架构，将每两个相邻的 Radix-2 stage 看做一组，每组对应 2bit，3/2/1/0，对应右移移位数，当 FFT 点数不是 4 的幂时，最后一个 stage 将仅对应 1bit；例如 对于 256 点，可设置为 [2 2 2 3]，对应缩放系数 $1/512$；512 点的 [2 2 2 2 2] 将是非法值，而 [0 2 2 2 2] 和 [1 2 2 2 2] 则是合法的。一般而言，对于标准的 FFT，总放缩系数为 $1/N$，即对于 Radix-2 的每个 stage 应移位 1bit，而 Radix-4 每个 stage 应移位 2bit，这也正是 FFT IP 的 SCALE_SCH 默认配置；通过设置 SCALE_SCH 可对输出结果进行缩放。

• 数据输入通道 tdata

  s_axis_data_tdata，传递待变换数据的实部 XN_RE 和虚部 XN_IM，最大支持 12 路 FFT 通道。XN_RE/XN_IM 位宽 8~34 可选，可以是定点数，也支持单精度浮点（32bit，IEEE-754 Single-precision Floating Point，包括 1bit 符号位、8bit 指数以及 23bit 底数）；浮点的输入输出数据类型将导致 IP 消耗更多的资源，以进行浮点-定点类型的转换。s_axis_data_tdata 的组织格式如下

当输入输出数据为定点数格式时，为二进制补码形式。

• 数据输出通道 tdata

  m_axis_data_tdata，传递变换结果的实部 XK_RE 和虚部 XK_IM，格式同 s_axis_data_tdata，如下

• 数据输出通道 tuser

  m_axis_data_tuser，具有 XK_INDEX、BLK_EXP 和 OVFLO 三个字段，结构如下

  XK_INDEX，位宽 ${\log }_2(maxmumpointsize)$，对应输出数据的索引。

  BLK_EXP，位宽 8bit，缩放系数，即右移的位数，每个通道均对应一个。仅在 Block Floating Point 模式下存在。

  OVFLO，位宽 1bit，断言本数据帧中是否出现了数据溢出异常，在新的数据帧被重置，每个通道均对应一个。

• 状态输出通道 tdata

  m_axis_config_tdata，具有 BLK_EXP 和 OVFLO 两个字段，结构如下

BLK_EXP 和 OVFLO 同 m_axis_data_tuser 的，不再赘述。

其他细节

关于计算精度及缩放系数

  在 Radix-4 架构下，每经过一个蝶形 stage，值会出现 $1+3\sqrt{2}\approx 5.242$ 的增长（对应 3bit 的增长），而 Radix-2 则会出现 $1+\sqrt{2}\approx 2.414$ 的增长 （对应 2bit 的增长）。为了 处理 bit 的增长可能带来的问题（计算溢出 or 精度损失）：

1. 采用不缩放的全精度定点计算，以携带所有重要的整数位到计算结束；
2. 每个 stage 采用固定缩放因子的定点计算；
3. 采用浮点运算（Block Floating Point）。

  若采用全精度不缩放定点计算，则每个蝶形 stage 的乘法计算结果会截断（Truncation）或四舍五入（Convergent Rounding），最终输出宽度 $InputWidth+{\log }_{2}NFFT+1$，以适应比特增长的最坏情况。如果采用定点缩放，则需要合理设置 SCALE_SCH，以避免溢出和精度损失。采用 Block Floating Point 浮点运算模式时则将自动确定缩放倍数，缩放倍数可通过 BLK_EXP 字段进行追踪。

  蝶形单元输出数据的舍入模式可在 IP 界面配置。其中截断（Truncation）可能引入直流偏置，而收敛舍入（Convergent Rounding）可以避免这一问题，但消耗更多的资源，并略微增加转换延迟。

  是否启用缩放和是否采用浮点运算可以在例化生成 IP 核时设置，是否采用浮点计算与输入输出数据的定点/浮点格式可独立设置。在定点缩放和浮点运算的模式下，输入输出数据位宽相一致；全精度不缩放定点模式下输出数据位宽如前所述。

  在 FPGA 端实现完全的浮点运算是非常昂贵的，因此 FFT IP 核在 Block Floating Point 选项采用的实际是更高精度的定点运算，通过相位因子宽度 Phase Factor Width 的选择控制计算的精度。当然，浮点运算将消耗更多的资源，因此如果输入数据的状况是被充分了解的、可以通过固定的缩放因子保证数据不出现溢出问题或严重的精度损失时，建议采用固定缩放因子的定点运算，以节约 FPGA 资源。

  由于舍入误差的存在，实值 FFT （输入值均为实数，即虚部均置零）的输出结果不完全对称，同时考虑到低频部分的噪声水平更加明显，因此 Xilinx 建议在执行实值 FFT 时使用输出结果的上半部分（$N/2+1\sim N$）。

计算溢出

  当启用 OVFLO 功能时，若数据帧转换结果中出现溢出，则 OVFLO 会在数据卸载期间保持为高；如果有多个 FFT 通道，则每个通道各自对应一个。在发生溢出时，转换结果发生绕回（wrapped）而不是饱和（saturated），这会导致却大多数应用下数据不可用。

架构选择

  FFT IP 提供了四种架构，以在核的大小和转换速度方面进行权衡。Pipelined Streaming I/O 模式允许连续进行数据处理，规模最大，速度最快；Radix-4 Burst I/O 分别加载数据和处理数据，他的规模比 Pipelined 更小，但转换时间也更长；Radix-2 Burst I/O 使用比 Radix-4 Burst I/O 更小的蝶形单元，规模更小，转换时间更长；Radix-2 Lite Burst I/O 在 Radix-2 Burst I/O 的基础上进一步采用时分复用，规模更小，转换时间更长。在核心规模和数据吞吐方面，每一种架构几乎都是后一个架构的两倍，需根据项目需求决定采用何种架构。

  除 Radix-4 Burst I/O 仅支持 64 到 65536 点的变换外，其余架构均支持 8 到 65536 点变换。

  Pipeline 的每个 Radix-2 蝶形单元可以同时加载新帧的数据、处理当前帧的计算、卸载上一帧的计算结果，因此支持连续的数据输入，并支持一定延迟后的连续数据卸载。（不过数据流连续并不代表可以忽略 AXI 的流等待状态，实际上 FFT 核可能会插入等待状态来中断样本输入）。其架构如下

  Radix-4 Burst I/O 的数据处理与加卸载分开进行，首先进行数据加载，待数据帧完整后，开始进行转换，转换完成后再进行数据卸载。若卸载采用位翻转，则数据加卸载可以同时进行（即加载新帧的同时，卸载当前帧的转换结果）；若输出采取自然序，则无法同时加载、卸载数据。关于转换数据的组织顺序，详见下一小节。Radix-4 Burst I/O 架构如下

  Radix-2 Burst I/O 与 Radix-4 类似，数据加卸载与数据处理不能同时进行，而如果采用位翻转，则数据加卸载可以重叠进行。

  Radix-2 Lite Burst I/O 采用共享的加减法器，时分复用，以时间性能的损失换取更小的资源消耗。数据加卸载和转换不可同时进行，如果采用位翻转，则数据加载、卸载可以重叠进行。Radix-2 Lite 的架构如下

数据顺序

  待转换数据以自然序输入（就是正常的二进制数顺序），而数据输出可以选择自然序（Natrual Order）还是反序（Bit/Digit Reversed Order）。默认设置下，输出数据采用反序。采用自然序数据输出会导致性能损失，在 Burst I/O 架构下，会导致时间性能损失，而对于 Pipelined Streaming I/O 架构，需要额外的 RAM 资源进行数据重排，导致消耗更多的面积。

  反序输出时，对于 Radix-2 和 Pipelined 架构，输出结果存在位翻转，最高位实际是计算结果的 LSB，如 4’b0001 实际表示 4’d8（4’b1000）；对于 Radix-4 架构，则翻转是以 2bit 为单元（官方称之为 Digit）进行的，如 4’b0001 实际表示 4’d4（4’b0100），如果结果位宽是奇数，则最低 1bit 代表最高位，其余以 Digit 为单位翻转，如 5’b00101 实际表示 5’d24（5’b11000）。

  输出序号 XK_INDEX 值为 $0\sim N-1$，对应 N 个转换结果输出。在反序输出模式下，这个序号也进行了翻转。

  当输出顺序为自然序时，可使用循环前缀插入。循环前缀：接受 FFT 输出的一部分，前缀到开头，随后再输出转换结果。CP_LEN 可修改循环前缀的长度，取值 $0\sim N-1$，默认 0。

实时/非实时模式

  默认采用的是非实时模式（Non Real Time），AXI 接口会用到 tdata、tready、tvaild 三组信号，从而使 Master 和 Slave 端均具有随时暂停数据传输的能力。在实时模式下（Real Time），IP 会生成一个更小、更快的设计，以牺牲数据加载、卸载的灵活性为代价。

  在 Real Time 模式下，数据输出通道、状态输出通道的 tready 信号被删除，同时数据输入通道的 tvaild 信号在数据帧开始以后也将被忽略（也就是用户可以在数据帧尚未开始加载时阻止数据传输，但一旦开始加载数据，用户将丧失中止数据传输的能力）。在这种设计下，数据加载过程中若出现需要暂停数据传输的情况，将会导致 IP 加载到错误的数据；在卸载数据时，若收端不能及时响应，也将造成数据的丢失。

  Real Time 设计下，若在开始数据加载后，用户没有提供数据（s_axis_data_tvalid = L），则 IP 将重复上一个数据，并拉高 event_data_in_channel_halt 以指示出现数据加载异常，如下图所示

这种异常同时会进一步导致后续帧的数据错位（比如图中当前帧的最后一个数据 D7 将被 IP 认为是下一帧的第一个数据），因此发生这一问题后就必须重新进行数据同步。

数据输入输出时序

  从数据输入到数据输出的延迟由所选架构和循环前缀长度共同决定。

• Pipelined Stream I/O & CP_LEN=0

  在流水架构下，数据加载、处理、卸载同时进行，如上图所示。

• Burst I/O

  在 Burst I/O 架构下，在处理新帧前，必须先卸载当前帧的转换结果；若采用自然序，数据加卸载不能同时进行，如图 3-44 所示，而如果采用输出反序，则数据加卸载可以同时进行，如图 3-45 所示。注意到在加载完一帧数据后，IP 仍加载了一小段下一帧的数据，这是由于 AXI 包含了一个大小为 16 的输入缓冲区，缓冲区的数据在 IP 加载数据时才会被加载到处理核心的寄存器中进行转换计算。

关于配置信息的应用时间节点

  在核心开始转换新数据帧时，会检测是否有新的配置信息，若在加载新帧之前或在加载新帧时的同一上升沿看到新的配置信息，则会使用该配置信息处理该帧数据，否则使用最后看到的配置信息（如果从没有过配置，则采用默认配置）处理该帧。为确保新帧采用新的配置信息进行转换，应在数据输入通道的第一个数据前至少一个 clk_cycle 写入配置信息。

  在 Pipeline Stream 架构下，若核心非空闲，则难以确定新的配置信息被应用于了哪个帧，因为数据很可能已经进入了处理流程。官方推荐使用帧启动信号 event_frame_started 进行配置信息到帧的同步，当该信号为高时，是为下一帧更新配置的安全时间点，在这个时间点后的配置则可能无法确定被用于哪一帧。

FFT IP 例化介绍

  在 Vivado IP 界面可以找到 Fast Fourier Transform 的 IP 核

  在配置完成后，可以看左侧的 Implementation Details，里面给出了接口里各个字段的定义，这个很方便，不需要根据配置自己去推算各个字段是什么含义了：

控制代码实现 & 测试

  作为示例，这里例化一个 Radix-2 Burst I/O 架构、单 FFT 通道、1024 点变换、16bit 定点、定点缩放、截断舍入、自然序输出、Non Real Time 的 IP，同时启用 XK_INDEX 以方便我们观察数据结果的顺序。

  用 matlab 生成数据文件，该信号包含 3，37，200 三个频点：

clc,clear,close all

%% 生成 testbench $readmemh 需要读取的 FFT IP 测试数据
N = 1024;
DataWidth = 16;

t = 0:N-1;
data = 2*cos(2*pi/N*3*t) + sin(2*pi/N*37*t) + cos(2*pi/N*200*t+pi/3);

data = data./max(abs(data))*0.99;   %归一化
data = round(data*2^(DataWidth-1)); %量化

fp = fopen('./fft_test_data.txt','w');
for i=1:N
    if data(i)>=0
        s = dec2hex(data(i),ceil(DataWidth/4));
    else
        s = dec2hex(data(i)+2^DataWidth,ceil(DataWidth/4));
        %这个函数太蠢了，如果数比较小，则不会拓展高位的符号位，所以要自己计算一下
    end
    
    fprintf(fp,s);
    fprintf(fp,'\n');
end
clear i s

fclose(fp);

%% fft
data_fft = fft(data)./N;

figure('color','w')
subplot(3,1,1)
plot(data)
subplot(3,1,2)
plot(real(data_fft))
subplot(3,1,3)
plot(imag(data_fft))

  然后编写 Testbench 如下：

//测试FFT IP核
`timescale 1ns/100ps
`default_nettype none
module FFT_IP_tb();

reg		clk_100M	= 1'b1;
reg		rst_n		= 1'b1;

always #5 begin
	clk_100M	<= ~clk_100M;
end

//--------------------------FFT IP------------------------------------
//config
reg		[23:0]	s_axis_config_tdata		= {3'b0, 20'h55555, 1'b1};
reg				s_axis_config_tvalid	= 1'b0;
wire			s_axis_config_tready;

//data in
reg		[31:0]	s_axis_data_tdata		= 32'd0;
reg				s_axis_data_tvalid		= 1'b0;
wire			s_axis_data_tready;
reg				s_axis_data_tlast		= 1'b0;

//data out
wire	[31:0]	m_axis_data_tdata;
wire	[9:0]	m_axis_data_tuser;
wire			m_axis_data_tvalid;
reg				m_axis_data_tready		= 1'b0;
wire			m_axis_data_tlast;

//events
wire			event_frame_started;
wire			event_tlast_unexpected;
wire			event_tlast_missing;
wire			event_status_channel_halt;
wire			event_data_in_channel_halt;
wire			event_data_out_channel_halt;

//FFT IP
FFT_0 FFT_inst(
	.aclk							(clk_100M),
	.aresetn						(rst_n),

	//Config
	.s_axis_config_tdata			(s_axis_config_tdata),
	.s_axis_config_tvalid			(s_axis_config_tvalid),
	.s_axis_config_tready			(s_axis_config_tready),

	//DATA INPUT
	.s_axis_data_tdata				(s_axis_data_tdata),
	.s_axis_data_tvalid				(s_axis_data_tvalid),
	.s_axis_data_tready				(s_axis_data_tready),
	.s_axis_data_tlast				(s_axis_data_tlast),

	//DATA OUTPUT
	.m_axis_data_tdata				(m_axis_data_tdata),
	.m_axis_data_tuser				(m_axis_data_tuser),
	.m_axis_data_tvalid				(m_axis_data_tvalid),
	.m_axis_data_tready				(m_axis_data_tready),
	.m_axis_data_tlast				(m_axis_data_tlast),

	//event signal
	.event_frame_started			(event_frame_started),
	.event_tlast_unexpected			(event_tlast_unexpected),
	.event_tlast_missing			(event_tlast_missing),
	.event_status_channel_halt		(event_status_channel_halt),
	.event_data_in_channel_halt		(event_data_in_channel_halt),
	.event_data_out_channel_halt	(event_data_out_channel_halt)
);

//----------------------------Ctrl-------------------------------------
reg		signed	[15:0]	fft_test_data	[0:1023];

reg		[9:0]	cnt		= 10'd0;
always @(posedge clk_100M or negedge rst_n) begin
	if(~rst_n) begin
		cnt		<= 10'd0;
	end
	else begin
		if(s_axis_data_tvalid & s_axis_data_tready) begin
			cnt		<= cnt + 1'b1;
		end
		else begin
			cnt		<= cnt;
		end
	end
end

//KN_real/img
reg		signed	[15:0]	KN_Real;
reg		signed	[15:0]	KN_Img;
always @(*) begin
	KN_Real		<= fft_test_data[cnt];
	KN_Img		<= 16'sd0;
end

//s_axis_data_tdata
always @(*) begin
	s_axis_data_tdata	<= {KN_Img, KN_Real};
end

//s_axis_data_tvalid
always @(posedge clk_100M or negedge rst_n) begin
	if(~rst_n) begin
		s_axis_data_tvalid	<= 1'b0;
	end
	else begin
		s_axis_data_tvalid	<= 1'b1;
	end
end

//s_axis_data_tlast
always @(*) begin
	s_axis_data_tlast = (cnt==10'd1023)? 1'b1 : 1'b0;
end

//m_axis_data_tready
always @(posedge clk_100M or negedge rst_n) begin
	if(~rst_n) begin
		m_axis_data_tready	<= 1'b0;
	end
	else begin
		m_axis_data_tready	<= 1'b1;
	end
end

//XK_Real
reg		signed	[15:0]	XK_Real;
integer i;
always @(*) begin
    //如用反序输出，则需要调整XK_Real、XK_Img、XK_INDEX三个信号的顺序
    //自然序的结果与matlab仅相差+-3以内；但倒序时出的结果太奇怪了，序号和数据都对不上
	for(i=0; i<=15; i=i+1) begin
		XK_Real[i]		<= m_axis_data_tdata[i];//[15-i];//
	end
end

//XK_Img
reg		signed	[15:0]	XK_Img;
always @(*) begin
	for(i=0; i<=15; i=i+1) begin
		XK_Img[i]		<= m_axis_data_tdata[16+i];//[31-i];//
	end
end

//XK_INDEX
reg				[9:0]	XK_INDEX;
always @(*) begin
	for(i=0; i<=9; i=i+1) begin
        XK_INDEX[i]		<= m_axis_data_tuser[i];//[9-i];//
	end
end

initial begin
    $readmemh("./fft_test_data.txt", fft_test_data);
	rst_n	= 1'b0;
	#200;
	rst_n	= 1'b1;

	#100000;
	$stop;
end

endmodule

  仿真结果如下，对比 matlab 结果，误差在 $\pm 3$ 以内

  由于这个自然序不能同时加卸载数据，所以后面又尝试了反序输出，然而 XK_REAL 和 XK_IMG 的结果非常奇怪，不论是计算结果还是应当出现的频点，都没法被理解。有懂的兄弟可以说一下这个反序输出到底是怎么回事 o(╥﹏╥)o。

参考文献

• pg109-xfft.pdf