radix16_dft_stage 模块深度解析

一句话概括

radix16_dft_stage 是素因子 FFT（Prime Factor Algorithm, PFA）流水线中的基-16 DFT 计算核心，它通过两个级联的 AIE 内核将 16 点 DFT 分解为高效的矩阵乘法运算，专门处理长度为 1008（=7×9×16）的复合长度 FFT 中最后的 16 点变换阶段。可以把整个 PFA-1008 FFT 想象成一条装配线：数据先经过 dft7 和 dft9 的粗加工，再通过转置重排，最后到达 radix16_dft_stage 进行精密的基-16 变换——就像把一个大工程拆分成多个小车间协作完成。

问题空间与设计动机

为什么需要这个模块？

在数字信号处理中，FFT 的长度往往是制约性能的关键。传统的 Cooley-Tukey 算法要求长度为 2 的幂次，但许多实际应用（如通信系统中的 OFDM）需要非 2 的幂次长度。素因子算法（PFA）允许我们将长度为 $N = N_1 \times N_2 \times ... \times N_k$ 的 DFT 分解为多个较短长度的 DFT，只要这些因子互质。

对于 PFA-1008（ $1008 = 7 \times 9 \times 16$ ），我们需要实现三个互质的短长度 DFT：7 点、9 点和 16 点。其中 16 点 DFT 是最复杂的，因为：

计算密度高：16 点意味着 256 个复数乘法（朴素实现）
精度要求高：AIE-ML 使用定点运算，需要仔细的舍入和饱和处理
吞吐量需求：整个流水线需要在严格的时序约束下运行

为什么不直接用库函数？

AIE 架构的特殊性决定了我们必须手写优化内核：

向量指令集：AIE 提供 512-bit 向量寄存器和专门的 MAC（乘累加）指令
内存层次：本地内存（LMEM）和流式 DMA 访问需要显式管理
确定性时序：实时 DSP 要求可预测的延迟，通用库无法满足

核心抽象与心智模型

类比：精密钟表匠的分工协作

想象 radix16_dft_stage 是一个制作精密钟表的工坊，有两个工匠（dft16_0 和 dft16_1）协同工作：

dft16_0（第一工匠）：负责读取原材料（输入样本），进行前半部分的齿轮加工（矩阵乘法的前半部分），然后将半成品通过传送带（cascade 连接）传递给下一个工匠
dft16_1（第二工匠）：接收半成品（累加器状态），完成后半部分的加工（矩阵乘法的后半部分），进行最终的抛光（右移归一化），产出成品（输出样本）

这种分工的核心洞察是：16 点 DFT 可以分解为两次 8 点的部分结果计算，中间通过 cascade 通道传递高精度的累加器状态（64-bit），避免在边界处丢失精度。

关键抽象层

┌─────────────────────────────────────────────────────────────┐
│                    dft16_graph (Graph 层)                    │
│              ADF 框架的图定义，描述拓扑结构                   │
│         ┌──────────┐      cascade       ┌──────────┐        │
│  sig_i ─┤ dft16_0  ├────────────────────→│ dft16_1  ├── sig_o │
│         └──────────┘                     └──────────┘        │
└─────────────────────────────────────────────────────────────┘
                              ↓
┌─────────────────────────────────────────────────────────────┐
│                  dft16_0 / dft16_1 (Kernel 层)               │
│              AIE 内核实现，使用 aie_api 向量指令             │
│                                                             │
│   MMUL<1,4,8> compute;  // 1×4 × 4×8 矩阵乘法单元           │
│   aie::vector<cint16,32> vc0-vc3;  // 旋转因子向量          │
└─────────────────────────────────────────────────────────────┘

架构与数据流

整体数据流图

详细数据流分析

1. 输入分发（Fan-out）

输入端口 sig_i 同时连接到两个内核的不同输入端口：

kk0.in[0]：接收偶数索引的样本向量（由迭代器 itr0 访问）
kk1.in[1]：接收奇数索引的样本向量（由 itr0 += 2 偏移后访问）

这是一种空间并行策略：同一个输入缓冲区被两个内核以不同的偏移量读取，没有数据竞争，因为每个内核只读取自己负责的子集。

2. Cascade 通道（精度保持）

kk0 的输出通过 connect<cascade> 连接到 kk1 的输入。Cascade 是 AIE 架构特有的高精度累加器传递机制：

数据类型：cacc64（64-bit 复数累加器）vs cint16（16-bit 复数整数）
目的：在第一级计算的部分结果传递到第二级时，保持完整的精度，避免过早量化
物理实现：直接连接相邻 AIE 核的专用 cascade 总线，低延迟、高带宽

3. 输出汇聚

kk1 完成最终计算后，通过 to_vector<TT_DATA>(DNSHIFT) 将 64-bit 累加器结果右移 15 位（DNSHIFT = 15），饱和截断为 16-bit 输出。

组件深度剖析

1. `dut_graph` —— 测试封装图

文件: dft16_app.cpp

这是 radix16_dft_stage 的顶层测试封装，继承自 ADF graph 类：

template<int X, int Y>
class dut_graph : public graph {
  dft16_graph dft16;  // 实例化核心图
public:
  input_plio  sig_i;  // PL 侧输入接口
  output_plio sig_o;  // PL 侧输出接口
  
  dut_graph(void) {
    // 创建 PLIO 端口，64-bit 位宽
    sig_i = input_plio::create("PLIO_i", plio_64_bits, "data/sig_i.txt");
    sig_o = output_plio::create("PLIO_o", plio_64_bits, "data/sig_o.txt");
    
    // 连接 PLIO 到 DFT16 图
    connect<stream>(sig_i.out[0], dft16.sig_i);
    connect<stream>(dft16.sig_o, sig_o.in[0]);
    
    // 物理位置约束（X,Y 为基准 tile 坐标）
    location<kernel>(dft16.kk0) = tile(X+3, Y+0);  // 相对偏移 +3,+0
    location<kernel>(dft16.kk1) = tile(X+4, Y+0);  // 相对偏移 +4,+0
  }
};

设计要点：

模板参数 <X,Y>：允许在更大的 AIE 阵列中灵活放置，通过编译时参数指定基准坐标
位置约束：kk0 和 kk1 必须放在相邻的 tile（X+3 和 X+4），以满足 cascade 连接的物理限制
运行时比例：runtime<ratio>(kk0) = 0.9 表示该内核占用 tile 90% 的计算资源，预留 10% 给数据传输开销

2. `dft16_graph` —— 核心计算图

文件: dft16_graph.h

class dft16_graph : public graph {
public:
  kernel kk0, kk1;           // 两个 AIE 内核实例
  port<input> sig_i;         // 图级输入端口
  port<output> sig_o;        // 图级输出端口
  
  dft16_graph(void) {
    // 创建内核对象，传入旋转因子系数表
    kk0 = kernel::create_object<dft16_0>(std::vector<typename dft16_0::TT_TWID>{DFT16_TWID});
    kk1 = kernel::create_object<dft16_1>(std::vector<typename dft16_1::TT_TWID>{DFT16_TWID});
    
    // 指定源码文件（用于编译链接）
    source(kk0) = "dft16_mmul0.cpp";
    source(kk1) = "dft16_mmul1.cpp";
    
    // 运行时资源分配
    runtime<ratio>(kk0) = 0.9;
    runtime<ratio>(kk1) = 0.9;
    
    // 连接定义
    connect<>(sig_i, kk0.in[0]);        // 输入 → kk0
    connect<>(sig_i, kk1.in[1]);        // 输入 → kk1（奇数样本）
    connect<cascade>(kk0.out[0], kk1.in[0]);  // kk0 → kk1（累加器）
    connect<>(kk1.out[0], sig_o);       // kk1 → 输出
    
    // FIFO 深度配置（缓冲流式数据）
    fifo_depth(n0) = 32;
    fifo_depth(n1) = 32;
  }
};

连接拓扑解读：

sig_i ──┬──→ kk0.in[0] ──[cascade]──→ kk1.in[0] ──┐
        │                                          ├──→ sig_o
        └──→ kk1.in[1] ────────────────────────────┘

注意到 kk1 有两个输入：in[0] 来自 kk0 的 cascade（累加器状态），in[1] 来自 sig_i（原始输入的奇数样本）。这是双源输入模式，允许 kk1 合并前级部分结果和新的输入数据。

3. `dft16_0` —— 第一级计算内核

文件: dft16_mmul0.h, dft16_mmul0.cpp

类型定义（关键常量）：

typedef cint16 TT_DATA;      // 输入/输出数据类型：16-bit 复数
typedef cint16 TT_TWID;      // 旋转因子类型：16-bit 复数
typedef cacc64 TT_ACC;       // 累加器类型：64-bit 复数

static constexpr unsigned NUM_FFT = 4*7*9;     // 252 次 DFT-16 变换
static constexpr unsigned NSAMP_I = 16*NUM_FFT; // 4032 个输入样本
static constexpr unsigned DNSHIFT = 15;         // 输出右移位数
static constexpr unsigned COEFF_DEPTH = 16*16;  // 256 个旋转因子

Run 方法执行流程：

void dft16_0::run(input_buffer<TT_DATA, extents<NSAMP_I>>& __restrict sig_i,
                  output_cascade<TT_ACC>* __restrict acc_o)
{
  // 1. 创建向量迭代器，每次读取 4 个 cint16（128 bits）
  auto itr0 = aie::begin_vector<4>(sig_i);
  
  // 2. 定义矩阵乘法单元：1×4 × 4×8 = 1×8
  using MMUL = aie::mmul<1, 4, 8, TT_DATA, TT_TWID>;
  MMUL compute;
  
  // 3. 预加载旋转因子到向量寄存器（避免重复内存访问）
  aie::vector<TT_TWID, 32> vc0 = aie::load_v<32>(&coeff[0]);
  aie::vector<TT_TWID, 32> vc1 = aie::load_v<32>(&coeff[32]);
  aie::vector<TT_TWID, 32> vc2 = aie::load_v<32>(&coeff[128]);
  aie::vector<TT_TWID, 32> vc3 = aie::load_v<32>(&coeff[160]);
  
  // 4. 主循环：处理所有 DFT-16 变换
  for (unsigned rr = 0; rr < NUM_FFT; rr++)
    chess_loop_range(NUM_FFT/8,)     // 向编译器提示循环次数
    chess_prepare_for_pipelining     // 启用软件流水线
  {
    // 第一次 MAC：计算输出的第 0-7 个频点（前半部分）
    compute.mul(*(itr0 + 0), vc0);   // 初始乘法
    compute.mac(*(itr0 + 1), vc1);   // 乘累加
    writeincr(acc_o, compute.to_accum());  // 写入 cascade
    
    // 第二次 MAC：计算输出的第 8-15 个频点（后半部分）
    compute.mul(*(itr0 + 0), vc2);
    compute.mac(*(itr0 + 1), vc3);
    writeincr(acc_o, compute.to_accum());
    
    itr0 += 4;  // 跳转到下一组 16 个样本（4 个向量 × 4 元素）
  }
}

关键设计决策：

MMUL<1,4,8> 配置：
- 第一个维度 1：输出向量的行数
- 第二个维度 4：输入向量和系数向量的共享维度（内积长度）
- 第三个维度 8：输出向量的列数（频点数的一半）
- 一次 MMUL 运算产生 8 个频点，两次覆盖完整的 16 点 DFT

系数布局（关键！）：

coeff[0:31]   → vc0 → 用于计算频点 0-7 的第一部分
coeff[32:63]  → vc1 → 用于计算频点 0-7 的第二部分
coeff[128:159]→ vc2 → 用于计算频点 8-15 的第一部分
coeff[160:191]→ vc3 → 用于计算频点 8-15 的第二部分

注意 coeff[64:127] 被跳过——这是留给 dft16_1 使用的另一半系数。

__restrict 关键字：向编译器保证指针不重叠，启用更激进的优化。

4. `dft16_1` —— 第二级计算内核

文件: dft16_mmul1.h, dft16_mmul1.cpp

Run 方法签名差异：

void dft16_1::run(
    input_cascade<TT_ACC>* __restrict acc_i,     // 新增：来自 kk0 的累加器输入
    input_buffer<TT_DATA, extents<NSAMP_I>>& __restrict sig_i,
    output_buffer<TT_DATA, extents<NSAMP_O>>& __restrict sig_o  // 改为输出缓冲区
)

执行流程对比：

特性	`dft16_0`	`dft16_1`
输入来源	`sig_i`（原始样本）	`acc_i`（cascade）+ `sig_i`（原始样本）
输出目标	`acc_o`（cascade）	`sig_o`（最终输出缓冲区）
样本偏移	`itr0`（偶数向量）	`itr0 += 2`（奇数向量）
系数范围	`[0:63]`, `[128:191]`	`[64:127]`, `[192:255]`
结果处理	`to_accum()`（保持精度）	`to_vector<TT_DATA>(DNSHIFT)`（量化输出）

核心逻辑：

// 迭代器初始化：从奇数向量开始
auto itr0 = aie::begin_vector<4>(sig_i);
auto itw0 = aie::begin_vector<8>(sig_o);  // 写入偶数输出向量
auto itw1 = aie::begin_vector<8>(sig_o);  // 写入奇数输出向量
itr0 += 2;   // 跳过前两个向量（偶数样本由 kk0 的逻辑对应）
itw1 += 1;   // 奇数输出向量偏移

for (...) {
  // 处理偶数输出位置（0, 2, 4...）
  compute = readincr_v<8>(acc_i);  // 读取 kk0 的部分结果
  compute.mac(*(itr0 + 0), vc0);   // 继续累加
  compute.mac(*(itr0 + 1), vc1);
  *itw0 = compute.template to_vector<TT_DATA>(DNSHIFT);  // 量化并写入
  itw0 += 2;  // 跳到下一个偶数位置
  
  // 处理奇数输出位置（1, 3, 5...）
  compute = readincr_v<8>(acc_i);
  compute.mac(*(itr0 + 0), vc2);
  compute.mac(*(itr0 + 1), vc3);
  *itw1 = compute.template to_vector<TT_DATA>(DNSHIFT);
  itw1 += 2;  // 跳到下一个奇数位置
  
  itr0 += 4;  // 下一组输入
}

5. `dft16_twiddle.h` —— 旋转因子表

文件: dft16_twiddle.h

这是一个纯头文件，包含预计算的 256 个旋转因子（ $W_{16}^{nk} = e^{-j2\pi nk/16}$ ），以 Q15 定点格式存储：

#define DFT16_TWID {\
{32767,0},      // W^0 = 1 + 0j (Q15 格式的 1.0)
{32767,0},      // 重复用于向量加载对齐
...
{30274,-12540}, // W^1 ≈ cos(π/8) - j*sin(π/8)
...
{-32768,0},     // W^8 = -1 + 0j
...
}

设计细节：

Q15 格式：16-bit 有符号整数表示 [-1, 1) 范围的浮点数，32767 ≈ 0.99997
冗余存储：某些值重复出现，确保 32 元素向量加载时的对齐要求
静态数组引用：通过构造函数传入引用，避免拷贝，内核间共享同一份只读数据

依赖关系与交互契约

上游依赖（谁调用我）

调用者	关系类型	契约说明
pfa1008_graph	组合	作为 PFA-1008 流水线的最后一级，接收 `transpose1` 的输出
`dut_graph`（本模块测试）	实例化	独立测试时使用，通过 PLIO 连接仿真数据

在完整 PFA-1008 图中：

// pfa1008_graph.h
connect<>(transpose1.sig_o, dft16.sig_i);  // 转置后数据进入 DFT16
connect<>(dft16.sig_o, sig_o.in[0]);       // DFT16 输出到 PL

下游依赖（我调用谁）

本模块是纯计算叶节点，不调用其他用户模块，仅依赖：

ADF 框架：<adf.h> 提供 graph, kernel, port 等基础设施
AIE API：<aie_api/aie.hpp> 提供向量指令封装（aie::mmul, aie::vector 等）

同级模块（协作关系）

模块	角色	与本模块的关系
dft7	PFA 第一阶段	处理 7 点 DFT，输出到 transpose0
dft9	PFA 第二阶段	处理 9 点 DFT，位于 transpose0/1 之间
transpose0/transpose1	数据重排	在 DFT 阶段之间进行矩阵转置，实现维度映射

设计权衡与决策记录

1. 为什么拆分为两个内核？

选择：将 16 点 DFT 拆分为 dft16_0 和 dft16_1 两个级联内核

替代方案：单内核完成全部计算

权衡分析：

维度	单内核方案	双内核级联方案（当前）
寄存器压力	高：需要同时保存 16 个频点的累加器	低：每个内核只处理 8 个频点
代码大小	大：循环展开程度高	适中：可维护性好
cascade 开销	无	有：但 cascade 是专用硬件，延迟极低
并行度	低：顺序执行	高：两核可流水线并行
精度	需中途量化	中间结果保持 64-bit 精度

结论：双内核方案在 AIE 架构上是更优选择，cascade 通道的硬件支持使得拆分几乎没有额外代价，而获得了更好的资源利用和精度保持。

2. 为什么选择 MMUL<1,4,8>？

选择：使用 1×4 × 4×8 的矩阵乘法粒度

替代方案：

MMUL<1,2,8>：更细粒度，更多循环迭代
MMUL<1,8,8>：更粗粒度，可能超出寄存器容量

决策依据：

AIE-ML 的向量寄存器宽度为 512 bits = 32 × 16-bit 元素
MMUL<1,4,8> 的输入：1×4（4 个 cint16 = 128 bits）+ 4×8（32 个 cint16 = 1024 bits，分两次加载）
输出：1×8（8 个 cacc64 = 1024 bits）
这个配置完美匹配 AIE 的 SIMD 能力，单次 MAC 指令完成 32 次乘加运算

3. 定点舍入策略

代码体现：

// 构造函数中设置全局舍入模式
dft16_0::dft16_0(...) {
  aie::set_rounding(aie::rounding_mode::symmetric_inf);
  aie::set_saturation(aie::saturation_mode::saturate);
}

选择：对称向无穷大舍入（symmetric_inf）+ 饱和模式

原因：

对称舍入：避免直流分量（全 1 输入）的系统性偏差
向无穷大：对正负数对称处理，符合 IEEE 754 风格
饱和：防止溢出时回绕，在信号处理中通常比模运算更安全

4. 位置约束的硬编码

代码体现：

location<kernel>(dft16.kk0) = tile(X+3, Y+0);
location<kernel>(dft16.kk1) = tile(X+4, Y+0);

风险：X+3 和 X+4 必须是物理上相邻的 tile，否则 cascade 连接会失败

缓解：在更大的系统（如 pfa1008_graph）中，通过统一的基准坐标 (X,Y) 管理，确保所有 tile 分配的一致性。

使用指南与示例

独立测试

// 实例化测试图，放置在 tile(18,0) 附近
dut_graph<18, 0> aie_dut;

int main(void) {
  aie_dut.init();     // 初始化图结构
  aie_dut.run(8);     // 运行 8 次迭代（每次处理 252 个 DFT-16）
  aie_dut.end();      // 结束运行，清理资源
  return 0;
}

集成到 PFA-1008

#include "dft16_graph.h"

class pfa1008_graph : public graph {
  dft7_graph dft7;           // 第一阶段
  transpose0_graph trans0;   // 第一次重排
  dft9_graph dft9;           // 第二阶段
  transpose1_graph trans1;   // 第二次重排
  dft16_graph dft16;         // 第三阶段（本模块）
  
public:
  pfa1008_graph(void) {
    // 数据流连接
    connect<>(dft7.sig_o, trans0.sig_i);
    connect<>(trans0.sig_o, dft9.sig_i);
    connect<>(dft9.sig_o, trans1.sig_i);
    connect<>(trans1.sig_o, dft16.sig_i);  // 进入 DFT16
    
    // 统一的位置约束（确保 tile 不冲突）
    location<kernel>(dft16.kk0) = tile(X+3, Y+0);
    location<kernel>(dft16.kk1) = tile(X+4, Y+0);
  }
};

边缘情况与陷阱

1. 输入数据格式要求

陷阱：输入文件 sig_i.txt 必须是特定的文本格式，每行一个样本：

(real, imag)
(12345, -6789)
...

验证：使用 gen_vectors.m MATLAB 脚本生成标准测试向量。

2. FIFO 深度配置

当前设置：fifo_depth(n0) = 32;

风险：如果上游数据源突发传输速率超过 DFT 处理速率，32 深度的 FIFO 可能溢出。

调试：在仿真中监控 fifo_status 信号，若出现 overflow 标志，需增加深度或优化调度。

3. 旋转因子表的对齐要求

陷阱：coeff 数组通过 alignas(16) 强制 128-bit 对齐，这是 AIE 向量加载指令的要求。

后果：若移除对齐属性，可能导致加载指令异常或性能下降。

4. 循环次数的编译期假设

代码：chess_loop_range(NUM_FFT/8,) 向编译器提示循环至少执行 31 次（252/8）。

风险：若运行时实际迭代次数少于提示值，可能导致错误优化；若远大于提示值，优化效果降低。

5. 多图实例化的 tile 冲突

场景：在同一 AIE 阵列中实例化多个 dft16_graph

风险：默认位置约束会导致 tile 冲突（多个内核争夺同一 tile）。

解决：为每个实例指定不同的基准坐标 (X,Y)，确保 tile(X+3,Y) 和 tile(X+4,Y) 不重叠。

性能特征

指标	数值	说明
单次调用处理量	252 × 16 = 4032 个复数样本	`NUM_FFT` 次 DFT-16 变换
内核周期估算	~252 × 8 = ~2016 周期	每次 DFT 约 8 次迭代，每次 2 次 MMUL
吞吐量	~2 样本/周期	考虑流水线填充和排空
资源占用	2 个 AIE tile，各 90% 利用率	`runtime<ratio> = 0.9`
内存占用	256 × 4 bytes = 1KB 旋转因子	只读，可被多个内核共享

总结

radix16_dft_stage 是 AIE-ML 架构上高效实现非 2 的幂次 FFT 的典型范例。其核心设计思想包括：

任务拆分：利用 cascade 通道将复杂 DFT 分解为多个简单内核，平衡寄存器压力和计算效率
向量化：通过 aie::mmul 充分挖掘 SIMD 并行度，单次指令完成 32 次乘加
精度管理：中间结果保持 64-bit，仅在最终输出时量化，最大化动态范围
显式控制：通过模板参数和位置约束，在编译期确定资源映射，满足实时系统的确定性需求

对于新加入团队的开发者，理解本模块的关键在于把握数据如何在两个内核间流动以及旋转因子如何被切分和复用。建议从修改 NUM_FFT 参数开始实验，观察循环行为和资源占用的变化。