🏠

Farrow分数延迟滤波器算法深度解析（基于Vitis-Tutorials）

1. 问题陈述

1.1 形式化定义

给定离散时间信号 \( x[n] \sim X(e^{j\omega}) \)，任意分数延迟 \( \tau = k + \mu \)（其中 \( k \in \mathbb{Z} \) 为整数延迟，\( \mu \in [0,1) \) 为分数延迟部分），要求生成输出信号 \( y[n] = x[n - \tau] \)，满足以下约束：

  1. 实时流式处理：输入/输出为连续样本流，允许有限因果延迟
  2. 动态延迟调整：分数延迟 \( \mu \) 可随时间变化，无需重新计算滤波器系数
  3. 硬件效率：在Xilinx Versal AIE架构上实现，最大化向量MAC利用率
  4. 数值稳定性：对于16位复数输入 \( x[n] \in \mathbb{C}^{16} \)，输出误差不超过1 LSB

1.2 经典问题边界

传统多速率滤波器（如多相FIR）仅支持固定有理因子 \( \frac{P}{Q} \) 的重采样，对于任意 \( \mu \) 需预计算海量系数表；而直接多项式插值（如拉格朗日）存在串行依赖，无法适配AIE的SIMD并行架构。

2. 直觉

2.1 朴素方法的失败

  1. 直接时域插值：对每个 \( \mu \) 实时计算 \( x[n - \mu] = \sum_{k} x[k] \cdot \text{sinc}(n - \mu - k) \)，sinc函数无限长，截断引入混叠
  2. 固定系数多相FIR：为每个可能的 \( \mu \) 预存系数表，当 \( \mu \) 分辨率为16位时，表大小达 \( 2^{16} \times N \)，完全超出AIE片上存储
  3. 串行霍纳法则：虽将多项式求值从 \( O(K^2) \) 降至 \( O(K) \)，但存在数据依赖（\( \mu \) 需逐层传递），无法并行化

2.2 核心洞察

转置Farrow结构：将分数延迟 \( \mu \) 的依赖从信号路径转移到系数路径，重构为 \( K \) 阶多项式形式的线性时变滤波器：

\[ y[n] = \sum_{m=0}^{K-1} \mu^m \cdot \sum_{k=0}^{N-1} c_{m,k} \cdot x[n-k] \]
其中：

  • \( K \) 为多项式阶数（通常取3或4）
  • \( c_{m,k} \) 为固定系数，与 \( \mu \) 无关
  • 内层和为 \( K \) 个并行FIR，外层为对 \( \mu \) 的多项式求值

3. 形式化定义

3.1 拉格朗日插值形式的Farrow滤波器

对于 \( K \) 阶拉格朗日插值，滤波器输出为：

\[ y[n + \mu] = \sum_{k=0}^{K} x[n - k + 1] \cdot \prod_{\substack{0 \leq m \leq K \\ m \neq k}} \frac{\mu + m - 1}{k - m} \]
将其展开为 \( \mu \) 的多项式：
\[ y[n + \mu] = \sum_{m=0}^{K} \mu^m \cdot f_m[n] \]
其中 \( f_m[n] = \sum_{k=0}^{K} c_{m,k} \cdot x[n - k + 1] \)，系数 \( c_{m,k} \) 由拉格朗日基函数的多项式展开得到。

3.2 转置结构的矩阵形式

令 \( \boldsymbol{\mu} = [1, \mu, \mu^2, \dots, \mu^{K-1}]^T \)，\( \boldsymbol{f}[n] = [f_0[n], f_1[n], \dots, f_{K-1}[n]]^T \)，则输出可表示为向量内积：

\[ y[n] = \boldsymbol{\mu}^T \cdot \boldsymbol{f}[n] \]
其中每个 \( f_m[n] \) 由独立的FIR滤波器生成，状态方程为：
\[ f_m[n] = c_{m,0}x[n] + c_{m,1}x[n-1] + \dots + c_{m,K-1}x[n-K+1] \]

3.3 约束条件

  1. 因果性：所有FIR滤波器 \( f_m[n] \) 仅依赖当前和历史输入 \( x[n], x[n-1], \dots \)
  2. 数值范围：\( \mu \in [0,1) \)，确保多项式求值无溢出
  3. 存储约束：固定系数 \( c_{m,k} \) 存储在AIE的系数内存（256位对齐）中，状态缓冲区大小为 \( K \times N \)（\( N \) 为FIR长度）

4. 算法

4.1 伪代码

Algorithm: Farrow Fractional Delay Filter (Transposed Structure)
Input: 
    x[n] ∈ ℂ¹⁶: Complex input sample stream
    μ[n] ∈ [0,1): Fractional delay stream
    c[m,k] ∈ ℝ¹⁶: Fixed Farrow coefficients (0 ≤ m < K, 0 ≤ k < N)
    f_state[m][0..N-1]: State buffers for K FIR filters
Output:
    y[n] ∈ ℂ¹⁶: Fractionally delayed output stream

Parameters:
    K: Polynomial order (e.g., 4)
    N: FIR filter length (e.g., 4)

Initialization:
    for all m in 0..K-1:
        for all k in 0..N-1:
            f_state[m][k] ← 0
    Set symmetric_inf rounding, saturate mode (AIE-specific)

Main Loop (per sample):
    1. Read inputs: x ← x[n], μ ← μ[n]
    2. Update FIR states (shift registers):
        for all m in 0..K-1:
            for k in N-1 downto 1:
                f_state[m][k] ← f_state[m][k-1]
            f_state[m][0] ← x
    3. Compute parallel FIR outputs f[m]:
        for all m in 0..K-1:
            f[m] ← MAC(c[m,0], f_state[m][0], ..., c[m,N-1], f_state[m][N-1])
    4. Compute μ powers (Horner form):
        μ_pows[0] ← 1
        for m in 1..K-1:
            μ_pows[m] ← μ_pows[m-1] * μ
    5. Compute final output (inner product):
        y ← 0
        for m in 0..K-1:
            y ← y + μ_pows[m] * f[m]
    6. Write output: y[n] ← y

4.2 执行流程图

flowchart TD Start([Start Per Sample]) --> Read[Read x[n] and mu[n]] Read --> UpdateState[Update K FIR State Buffers] UpdateState --> ParallelFIR[Compute K Parallel FIR Outputs f[m]] ParallelFIR --> MuPowers[Compute mu Powers via Horner] MuPowers --> InnerProd[Compute Inner Product mu^T * f] InnerProd --> WriteOut[Write y[n] to Output Stream] WriteOut --> CheckDone{More Samples?} CheckDone -->|Yes| Read CheckDone -->|No| End([End]) subgraph AIE_Optimizations[AIE-Specific Optimizations] VecMAC[Vector MAC for Parallel FIR] VecLUT[Optional mu Quantization + LUT] CyclicBuf[Cyclic State Buffers] end ParallelFIR -.- VecMAC MuPowers -.- VecLUT UpdateState -.- CyclicBuf

5. 复杂度分析

5.1 时间复杂度

阶段 操作 复杂度（每样本）
状态更新 移位寄存器 \( O(K \cdot N) \)（AIE中通过循环缓冲区优化为 \( O(1) \) 指针更新）
并行FIR 向量MAC \( O(K) \)（单周期内完成 \( K \) 个并行MAC，\( N \leq 8 \) 时）
\( \mu \) 幂次 标量乘法 \( O(K) \)
内积 向量MAC \( O(1) \)（单周期向量内积）

整体时间复杂度：

  • 理论（串行）：\( O(K \cdot N) \)
  • AIE优化（并行）：\( O(1) \) 每样本（当 \( K \leq 8, N \leq 8 \) 时，适配AIE向量宽度）

5.2 空间复杂度

存储类型 大小
固定系数 \( K \times N \times 16 \) 位
状态缓冲区 \( K \times N \times 32 \) 位（复数16位实部+虚部）
临时变量 \( O(K) \)（存储 \( \mu \) 幂次和FIR输出）

整体空间复杂度：\( O(K \cdot N) \)，完全适配AIE片上存储（\( K=4, N=4 \) 时仅需128字节系数+256字节状态）。

6. 实现说明（基于Vitis-Tutorials代码）

6.1 代码与理论的差异

  1. 循环缓冲区替代移位寄存器：
    • 理论中状态更新为显式移位（\( O(K \cdot N) \)）
    • 实现中使用AIE cyclic_add 硬件指令，仅更新指针（\( O(1) \)），旧数据自然被覆盖
  2. 量化 \( \mu \) 与查找表（Final版本）：
    • 理论中实时计算 \( \mu^m \)
    • 实现中对 \( \mu \) 进行8位量化（256个条目），预计算所有 \( \mu^m \) 并存入LUT，单周期向量查表
  3. 向量FIR并行化：
    • 理论中 \( K \) 个FIR独立计算
    • 实现中使用AIE 256位向量寄存器，将 \( K \) 个FIR的系数和状态打包为向量，单周期完成所有MAC

6.2 工程折衷

  1. 多项式阶数选择：
    • 教程中选择 \( K=4 \)（4阶拉格朗日插值），在精度和硬件开销间取得平衡
    • \( K=3 \) 时节省25%存储，但通带波纹增加约3dB
  2. 舍入与饱和模式：
    • 教程中显式设置 aie::rounding_mode::symmetric_inf 和 aie::saturation_mode::saturate
    • 对称无穷舍入比默认截断减少约1 LSB的误差
  3. 状态初始化：
    • 所有状态缓冲区初始化为零，避免冷启动时的瞬态噪声
    • 初始版本（farrow_initial）为每个FIR独立初始化状态，最终版本合并为单一循环缓冲区数组

6.3 关键代码片段分析

// 来自 farrow_initial/farrow_kernel.cpp 的状态初始化
farrow_kernel::farrow_kernel( void )
{
  aie::set_rounding(aie::rounding_mode::symmetric_inf);
  aie::set_saturation(aie::saturation_mode::saturate);
  // 为4个并行FIR（K=4）分别初始化状态缓冲区
  for(int i=0;i<STATE_LEN;i++)
  {
    f3_state[i] = (cint16){0,0};
    f2_state[i] = (cint16){0,0};
    f1_state[i] = (cint16){0,0};
    f0_state[i] = (cint16){0,0};
  }
}

7. 对比

7.1 与经典多相FIR的对比

特性 Farrow分数延迟滤波器 经典多相FIR
延迟灵活性 任意 \( \mu \in [0,1) \)，动态调整 固定有理因子 \( \frac{P}{Q} \)
系数存储 \( K \times N \) 固定系数 \( P \times N \) 相系数（需预存所有相）
计算复杂度 \( O(K) \) 每样本 \( O(N) \) 每输出样本
AIE适配性 完美适配SIMD并行 需多相分解后并行，控制复杂
数值精度 取决于多项式阶数 \( K \) 取决于滤波器长度 \( N \)

7.2 与直接拉格朗日插值的对比

特性 转置Farrow结构 直接拉格朗日插值
并行性 \( K \) 个FIR完全并行 串行依赖（需先计算基函数）
\( \mu \) 更新 仅需更新外层多项式 需重新计算所有基函数
硬件实现 向量MAC友好 标量乘法为主，难以并行
AIE性能 1样本/周期（\( K=4 \)） ~4样本/周期（串行霍纳法则）

7.3 与Vitis-Tutorials中其他版本的对比

版本 优化技术 性能（样本/周期） 存储开销
farrow_initial 直接多项式求值 ~0.25 低
farrow_opt1 霍纳法则 ~0.5 低
farrow_opt2 转置结构 ~1 中
farrow_final \( \mu \) 量化+LUT+向量MAC ~1（峰值） 中高

参考文献

  1. Farrow, C. W. (1988). A continuously variable digital delay element. IEEE International Symposium on Circuits and Systems (ISCAS).
  2. Vitis-Tutorials: AI Engine Development — Farrow Filter Design Tutorial (Xilinx, 2024).
  3. Versal AI Engine Architecture Manual (AM001, Xilinx, 2024).
On this page