基于FPGA的超长傅里叶变换实现（DDR4方案）

前言

       在本文中，我想谈谈超长快速傅里叶变换算法在 FPGA 上的实现。促使我写这篇文章是出于分享我个人实践经验的愿望，我不想失去它，只在我的脑海中留下信息。本文只讲FPGA及Vivado实现方法，不详细提及代码。
       本文展示了即使在最现代的 FPGA 晶体管上也不可能实现非常长的 FFT 的 “经典” 方案，并提出了一种允许您执行此作的算法。该算法的主要思想也逐步考虑：从数学组件到使用外部 DDR 存储器创建基于 FPGA 的完整解决方案。本文将涉及为此类任务设计多通道处理系统的复杂性，特别是描述了我的实践经验。

一、实验需求

       想象一下，您需要开发一个具有高频分辨率的宽带频谱分析仪。众所周知，FFT 长度越长，分辨率越高。也就是说，您需要设计一个系统，该系统接收来自 ADC 的信号，在 FPGA 内部进行一些处理，并通过高速接口（例如 PCIe、USB 等）将数据输出到数据处理的下一阶段。

让我们正式确定任务的要求：

• 算法：快速傅里叶变换
• FFT 长度从 256K 到 256M 点。
• 实时持续数据处理。
• 独立输入通道的数量从 1 到 N。
• ADC 信号的位深度为 16 位。
• 全面的输入信号。
• 单个 ADC 通道的数据速率为 200 MHz 或更高。
• 在 FFT 链路前面存在可切换窗口功能。

二、传统方法验证

经典方法行不通！

       在继续描述超长傅里叶变换算法之前，应该说在 FPGA 上实现一定长度的经典 FFT 方案从根本上是不可能的。正是因为这个原因，本文中将要描述的解决方案出现了。

       假设你面临的任务是 开发一个长度为256K到64M个采样点的FFT（快速傅里叶变换）节点。 目前，领先的供应商AMD（Xilinx）和Intel（Alter）并没有免费提供长度超过64K点的FFT核心。原因在于，随着FFT长度的增加，FPGA芯片中的块存储器（BRAM）消耗会变得非常巨大。即使Xilinx的芯片中出现了URAM（UltraRAM）块，也无法解决这个问题。即使你通过经典的FFT算法，通过蝶形节点和交叉交换节点的流水线连接，编写了一个超级优化的解决方案，你也无法在FPGA上实现它。

       为什么？

       当FFT的长度增加N倍时，占用的BRAM资源量至少也会按比例增加N倍。对于浮点格式的FFT，BRAM节点的数量会按比例增加，但对于定点格式的FFT，资源消耗的评估则更加悲观。这是因为在FFT计算的每个阶段（每次蝶形运算之后），中间数据的位数会增加1bit。这种位宽的增长导致定点格式的FFT需要更多的存储资源来容纳中间结果，从而进一步加剧了资源消耗的问题。
       FFT 单元输出处的最终位宽根据以下公式计算：

       例如，1M FFT 为输入信号增加 20 位。如果输入信号的位深度为 16 位，则输出信号将分别为 36 位。

       让我们以Xilinx为例，看看64K点FFT占用的资源表。粗略估计，如果你需要实现一个1M点的浮点FFT模块，那么你将需要大约400 * 16 = 6400个RAMB36K类型的块存储器单元，或者约220 Mbit的存储资源！只有顶级的Virtex UltraScale+（例如VU29P）才具备这样的资源，而这还依赖于URAM（UltraRAM）单元的支持。换句话说，实现如此大规模的FFT模块对硬件资源的要求极高，通常只有最高端的FPGA芯片才能满足这种需求。

vivado报错：

下表以 Kintex 系列为例，显示了现在Xilinx FPGA 有多少内存；
Kintex Ultrascale：

Kintex Ultrascale+

       正如所见，FPGA的块存储器资源非常有限，因此即使实现1M点长度的FFT，在传统架构下也几乎是不可能的。这种限制主要源于FFT算法在计算过程中对存储资源的高需求，尤其是随着点数增加，中间数据的存储和位宽扩展会显著增加资源消耗。因此，实现大规模FFT通常需要采用更优化的架构设计或依赖更高端的FPGA硬件资源（如URAM），甚至可能需要结合外部存储器或其他硬件加速方案来满足需求。

三、数学优化（超长 FFT）

       幸运的是，可以利用一个数学技巧，将一维FFT（快速傅里叶变换）基于二维FFT来实现。
       该算法的核心思想是将长度为N的信号向量分解为N1和N2个采样点（其中N1和N2是2的幂次方）。这个向量会被重新排列成一个维度为N1 x N2的矩阵，所有计算都在这个矩阵上进行。具体步骤如下：
FFT 的计算公式：

序列 N 的维度除以 N1 和 N2 的乘积，则：

然后公式采用以下形式：

转换转折乘数后：

结果：

       可以注意到，长度为N的一维FFT被转化为两个长度分别为N1和N2的FFT，并且需要将第一个FFT的结果乘以旋转因子（twiddle factors）。
超长 FFT 算法的结构图：

四、算法实现

具体步骤：

• 将一维序列转置为 K1 x K2 矩阵
• 按 K1 行 K2 乘以计算 FFT
• 转置结果矩阵
• 通过转动乘数乘以结果
• 按 K2 列 K1 乘以计算 FFT
• 将矩阵转置为原始一维向量

       因此，使用三个数据转置节点和两个短长度 FFT 单元，可以实现长 FFT 的计算。
       示例：采样长度为 1M 的 FFT 内核。可以将计算拆分为两个 1K 计数的 FFT：1024 x 1024 = 1048576。下图显示 1024 点 FFT 节点只需要 7 个 RAMB36K 个单元格。

       可以看出，小型 FFT 内核实际上并不占用 FPGA 资源，特别是它们几乎不使用 block memory。让我们来看看算法的所有元素，看看每个元素的主要功能。

分步实施：

1、FFT 模块

       FFT节点通常基于经典的流水线架构实现，即通过Radix-2或Radix-4蝶形运算单元和交叉交换节点（cross-communication nodes）的log₂(N)级联连接而成。开发者可以选择使用FPGA厂商提供的免费现成IP核，或者编写自己的优化内核。FFT节点可以根据任务需求，支持定点格式（fixed-point）、浮点格式（floating-point），甚至是某种自定义格式（custom format）进行计算。
FFT 单元：

• 流水线、流式 I/O
• 固定未缩放/浮点
• 基数-2 / 基数-4
• 位反转/自然顺序
• 复数乘法器：4-DSP 结构

       需要明确的是，FFT（快速傅里叶变换）是一种积分算法，这意味着在FFT的输出端，数据的位宽会增加（同时也会累积舍入误差）。这就引出了一个问题：所选数据格式的计算精度是否满足你的需求？根据我的经验，对于超长FFT（例如长度为1M点），浮点运算（IEEE 754）的舍入误差会由于尾数位数较少（25位）而破坏最终频谱。因此，我不得不在所有处理阶段使用定点格式（fixed-point）的FFT。在某些特殊情况下，可以在FPGA的DSP48模块上编写自定义的浮点数学运算，以扩展尾数位数（虽然实现复杂，但结果值得）。

解释一下IEEE 754舍入误差：
IEEE 754单精度浮点数的尾数（mantissa）只有23位（实际有效位数是24位，包括隐含的1位）。尾数位数决定了浮点数能够表示的精度，位数越少，精度越低。FFT是一种积分算法，计算过程中涉及大量的乘法和加法运算。每次浮点运算都会引入微小的舍入误差，这些误差在超长FFT中会逐渐累积。对于长度为1M点的FFT，计算量非常大，误差累积效应会更加明显。对频谱的影响：由于误差的累积，FFT计算结果的精度会下降。这种精度下降会体现在频谱中，可能导致频谱失真或出现虚假的频率成分，从而破坏频谱的准确性。
举例说明：
假设你在处理一个长度为1M点的信号，使用单精度浮点数进行FFT计算。由于尾数位数有限，每次蝶形运算都会引入微小的误差。随着计算步骤的增加，这些误差会逐渐累积，最终导致频谱中出现不应有的噪声或失真。

回到正文：
       FFT核可以被实现为结果以比特反转顺序输出，或者以自然（正常）顺序输出。在第一种情况下，你会更加节省块状内存资源，但在矩阵转置时会稍微复杂化数据排列算法。如果你刚开始实现超长FFT，那么从第二个选项开始会更好，因为它调试起来更简单。
       如果你正在实现一个多通道系统，并且FFT核心是你自己编写的，而不是使用FPGA制造商提供的现成解决方案，那么你可以额外节省BRAM存储用于存放FFT蝶形运算所需的旋转系数。对于N个独立的并行通道，只需要使用1个旋转因子存储模块。也就是说，系统的通道数越多——资源节约就越多。

2、旋转乘数

       根据算法流程，在将矩阵输入到第二个FFT模块之前，数据需要乘以旋转因子（twiddle factors）。这可以通过两种方式实现：使用DDS（直接数字合成）或采用CORDIC算法。 第一种方法需要存储大量数据，这会占用大量FPGA的块存储器（BRAM），而这正是我们从文章一开始就在努力优化的问题。理论上，可以使用泰勒级数近似来减少BRAM中存储的数据量。但在我的实践中，这种方法会由于旋转因子的阶梯状信号形式而导致最终频谱失真。

       第二种基于CORDIC的方法完全不需要BRAM，因为它采用了一种迭代的移位和加法/减法操作方案。CORDIC算法的缺点在于计算下一个值所需的时间较长（通常需要20到30个时钟周期，具体取决于数据位宽）。这一缺点导致需要额外的数据延迟线，这会占用一定的FPGA逻辑资源。例如，对于一个位宽为512（2个复数采样，每个32位，8个通道）的多通道系统，额外需要512 * 30 = 15,000个触发器。在Xilinx FPGA中，可以使用SLICEM单元来构建基于移位寄存器的延迟线。或者，也可以使用几个BRAM块来实现延迟线。
CORDIC算法原理详解：
CORDIC(Coordinate Rotation Digital Computer, 坐标旋转数字计算方法)算法就是一种化繁为简的算法，通过基本的加减和移位运算代替乘法运算，逐渐逼近目标值，得出函数的数值解。

回到正文：
对 CORDIC 核心提出了以下要求：

• 并行处理方案
• 输入数据的位宽不小于 FFT 最终长度的 log2（N），以避免信号频谱失真。

对于多通道方案，也可以使用一个旋转乘法器节点来处理整个数据流，从而节省 FPGA 资源。

3、转置节点

       内存控制器用于存储中间数据向量，并在算法的各个阶段对矩阵进行转置。这可以是任何可用的内存类型：QDR SRAM、DDR3 或 DDR4 SDRAM。在我的项目中，我使用了后两种。但它们的工作原理是相同的：内存控制器对数据进行转置——按行获取一批数据，然后按列输出一批数据。

       从算法流程图中可以看出，此任务需要三个外部内存节点，这些节点需要满足两个主要要求。

1、首先

       内存必须存储至少 2 个 FFT 长度的最小数据块。这是为了在写入一批数据（按矩阵行）的同时，能够成功读取另一批数据（按矩阵列）。然而，最佳的数据存储方案不是通过多路复用器，而是当外部内存以 FIFO（先进先出）的方式实现时。在这种情况下，外部内存可以存储多个长度为 N 的数据块，并高效利用其资源。

       此外，在实践中，这种方案可以有效应对 FFT 节点输出接口的短暂中断。特别是当 PCIe 总线传输出现短暂中断时，FIFO 模式下的内存溢出概率显著低于在两种数据块之间切换的方案。在我实现的项目中，当 PCIe 总线出现中断时，DDR 内存在多路复用器模式下几乎总是会溢出，而在 FIFO 模式下则从未发生过溢出。

       我们来计算外部存储器中存储数据所需的容量。假设输入数据的位宽为 32 位（单精度浮点数，single floating-point），信号为复数形式（I/Q），FFT 的长度为 1M 个采样点，且采用“乒乓”模式作为最低要求。那么，存储数据将需要 2×32/8=8 MB2×32/8=8MB 的内存。如果采用深度为 32 的 FIFO 模式存储数据，则需要 256 MB 的内存（即1Gb，很多开发板上面的颗粒都是2Gb，可以满足）。

2、第二

       内存必须能够实时读写数据。根据算法，控制器输入和输出处的数据传输方案不同。因此，有必要以最快的速度正确组织数据传输过程，以免数据出现间隙和失真。为此，即使在设计或购买可以使用超长 FFT 的现成电路板的阶段，您也需要计算整个系统的最大吞吐量。也就是说，根据数据表找到流向外部存储器的最大 I/O 流量，并确定来自多通道 ADC 的数据馈送速率。FFT的长度并不影响带宽（吞吐量）。

       此外，对FPGA芯片提出了与三个内存控制器交互的要求。例如，每个SO-DIMM DDR4 SDRAM x64内存控制器需要三个FPGA Bank（相当于约150个芯片的物理I/O引脚）。总共至少需要450个I/O端口或9个HP（高性能）FPGA Bank，这不包括多千兆位收发器的Bank和配置Bank0。

DDR4 SDRAM IP 内核设置示例：

3、内存控制器交互算法：

• input stream 以处理速率进入一个小的 FIFO ，并以 memory controller 频率（例如 333 MHz）转换为 stream；
• 数据以直接形式写入内存，DDR 内存地址计数器呈线性增加。
• 在写入大小为 N = K1 * K2 的突发后，它以相反形式读取。
• DDR 读存储器地址计数器根据特定算法以指定步骤进行更改，以允许转置 K1 x K2 矩阵。
• 在读取的同时，对另一个空闲区域的外部存储器的写入仍在继续。
• 读取数据流从 controller 频率转换为 processing frequency （使用较小的 FIFO）。
• 如果由于某种原因 read 端冻结，则写入不会停止并继续，直到 controller 设置 FULL FIFO 标志。

在 FIFO 模式下与外部存储器交互的最终方案如下：

       需要注意的是，为了实现最大性能，在读取时需要遍历FPGA中内存控制器的所有FSM组（参见Xilinx DDR MIG的数据手册），也就是说需要访问所有bank和bank组的物理内存。这带来了额外的限制，并导致需要在FPGA中的内存控制器之后设置数据缓冲区（有一个 data buffer）。其作用是组织和打包从每个内存bank处理的事务，以便进一步传输到电路的下一个节点。对于数据缓冲区的实现，Xilinx现代FPGA中的URAM模块（深度为4K、位宽为72位的块）是理想的选择。

我还注意到，最后一个转置节点可以用作累加器，实现一种累加频谱分量以平均结果频谱的方案。

4、窗口函数及实践

       输入信号上的窗口叠加是可选的，但如果您需要向长 FFT 添加窗口函数，则可能会遇到困难。

       由于我们使用的是超长 FFT，因此没有地方存储窗口函数的数组。N 长度的 FFT 将需要大量的 block memory 资源。例如，考虑到其对称性，长度为 32M 的 16 位信号形式的窗口函数将需要读数：32 * 4 / 2 = 256 Mbit。即使对于高端 FPGA ，这也是很多。如果您希望能够连续切换功能（至少两个独立的数据缓冲区）怎么办？

       通过使用 Blackman-Harris 窗口的顺序和计算窗口系数的标准公式，这个问题的解决方案非常简单。通过使用众所周知的 CORDIC 生成所需频率的谐波信号，可以在不使用 FPGA 块存储器的情况下实现任何阶次的 Blackman-Harris 窗口函数！

窗函数的阶数越高，频谱的旁瓣水平越低。实际上，必须使用最高 5 阶的窗口。

       信号流经超长FFT算法的节点如下所示：选择某一值的峰值信号和线性调频信号作为输入信号。
Delta Function：


       Ultra-Long FFT 的每个阶段生成数据的脚本。它是用 Python 编写的，需要 numpy、scipy 和 matplotlib 库才能工作：

通过网盘分享的文件：Ultra-Long FFT脚本.py
链接: https://pan.baidu.com/s/1-6KFIpK3wz3A9hy2VaAIQQ?pwd=ymmp 提取码: ymmp
–来自百度网盘超级会员v7的分享

FFT 的一个阶段之前和之后的转置线性调频信号：

线性调频信号通过 FFT 单元：

5、个案研究

假设你有一个包含以下输入的任务：

• 输入信号是复数，16 位。
• 通道数为 1。400 MHz 采样率。
• IEEE-754 格式的中间数据。
• FFT 长度：1600 万点。无窗口函数。

计算：
1、输入数据流：（2 * 16） * 400e6 / 2^30 = ~1.5 GB/s。
2、中间信息流：~3 GB/s。
3、存储一个 FFT 包的内存容量：（2 * 32 / 8） * 16M = 128 MB。
4、读写内存带宽至少为 ~6 GB/s。
5、DDR4-2400 SDRAM x32 根据以下公式满足此要求：频率 * 双倍速率 *（宽度/字节）= 1.2e9 * 2 * （32 / 8） / 2^30 = 8.94 GB/s。

       通过IP Catalog，我们将创建DDR、FFT和CORDIC核心。完整的16M点FFT被分解为4K x 4K的FFT矩阵。假设在DDR控制器前后需要对数据进行重新打包，这需要占用4个RAMB36K块存储单元的FIFO。
       对实现算法所需的FPGA资源进行粗略估算。我们需要3个内存控制器、2个4K点FFT节点、1个CORDIC模块、6个FIFO块（位于内存控制器前后）以及3个内存bank缓冲区。为了简化计算，暂不考虑用于流同步的延迟线、复数乘法器以及其他逻辑。

综合表：

       如你所见，该项目不占用太多资源，并且适合相对便宜的 FPGA。但是，不要忘记 3 个内存控制器需要一定数量的 FPGA I/O 端口（7系列稍微贵些的基本能满足），并非所有芯片都能正常工作。

       本文展示了在实时高精度信号频谱分析任务中，如何在FPGA上实现超长FFT节点的方法。设计此类算法需要特定的硬件配置——FPGA必须连接三个独立的内存控制器。然而，这使得能够实现点数从256K到256M的超长FFT电路。由于算法在硬件实现中有许多细节，因此需要提前计算电路中所有节点的参数，并确保所选配置能够实现该核心。