1. 项目缘起当“固定”的滤波器遇上“变化”的世界在数字信号处理的世界里我们常常面临一个经典的矛盾设计的灵活性与实现的效率。就拿一个最常见的需求来说——你需要一个截止频率可变的低通滤波器。最直接的想法是什么没错每次频率参数一变我就重新计算一遍滤波器的系数然后更新到我的滤波器结构中。这在离线处理或者对实时性要求不高的场合或许可行但一旦放到FPGA、DSP或者高性能嵌入式系统中频繁地重新计算并加载一组可能多达几十甚至上百个的滤波器系数比如一个高阶FIR滤波器带来的计算开销和延迟往往是不可接受的。尤其是在软件无线电、雷达信号处理、音频效果器这类对实时性要求极高的领域参数需要连续、平滑地调整这种“硬切换”系数的方式会引入信号的不连续产生可闻的咔嗒声或影响系统性能。那么有没有一种方法能让滤波器的某个关键参数比如截止频率、分数延迟像拧旋钮一样连续可调而无需在每次调整时都大动干戈地重新计算所有系数呢这就是Farrow结构滤波器要解决的核心问题。它本质上是一种实现可变分数延迟滤波器或可变参数滤波器的高效结构。我第一次在项目中接触它是在为一个多速率采样系统设计一个采样率转换模块时传统的多项式插值方法在精度和复杂度上难以平衡而Farrow结构提供了一种优雅的折中方案。今天我们就来彻底拆解这个在工程上极具魅力的结构从它要解决什么问题到它的核心原理再到如何一步步设计并实现它最后分享一些在硬件实现中容易踩的坑。2. Farrow结构的核心思想将参数变化“固化”到结构里要理解Farrow结构我们得先忘掉那些复杂的公式从一个更直观的视角来看。想象一下一个滤波器的输出是其系数与输入信号卷积的结果。如果这个滤波器的某个特性如群延迟需要连续变化传统做法是让系数本身成为这个变化参数的函数。也就是说每一个滤波器系数h[k]不再是固定值而是一个关于可变参数d通常代表分数延迟范围在0到1之间的函数h[k] f_k(d)。Farrow结构的巧妙之处在于它对这个函数f_k(d)做了一个关键的假设和简化每个系数关于可变参数d的变化可以用一个低阶多项式来近似。也就是说h[k] ≈ c_{k,0} c_{k,1} * d c_{k,2} * d^2 ... c_{k,L} * d^L这里L是多项式的阶数c_{k,l}就是我们需要预先计算并固定下来的“子滤波器”系数。看到了吗变化的部分d及其幂次被抽离出来了而需要存储和参与实时卷积运算的变成了固定不变的系数c_{k,l}。这个思想带来了巨大的优势实时性当需要改变延迟d时我们不再需要重新计算或加载任何滤波器系数c_{k,l}。我们只需要根据新的d值实时计算出一组“权重”(1, d, d^2, ..., d^L)然后用这组权重去组合那些固定的子滤波器的输出。结构化与模块化整个滤波器可以被分解为(L1)个并行的、系数固定的子滤波器每个对应多项式的一项以及一个多项式求值即权重组合模块。这种结构非常规整特别适合用硬件如FPGA进行并行流水线实现。设计灵活性我们可以通过设计多项式阶数L和子滤波器系数c_{k,l}来权衡逼近精度、计算复杂度和滤波器性能。那么这个结构具体长什么样呢一个典型的、用于实现分数延迟的Farrow结构框图如下所示我们以三次多项式即L3为例输入 x[n] | |---- 固定子滤波器 H0(z) (系数为 c_{k,0}) ---- 乘 1 (即 d^0) ---\ |---- 固定子滤波器 H1(z) (系数为 c_{k,1}) ---- 乘 d ----- 加法器 ---- 输出 y[n] |---- 固定子滤波器 H2(z) (系数为 c_{k,2}) ---- 乘 d^2 -----/ ---- 固定子滤波器 H3(z) (系数为 c_{k,3}) ---- 乘 d^3 -----/关键点H0(z),H1(z), ...,H3(z)都是普通的、系数固定的FIR滤波器。它们的输出分别乘以d^0,d^1,d^2,d^3然后求和就得到了最终具有分数延迟d的输出y[n]。d可以在每个采样时刻动态改变而所有子滤波器的系数是焊死在硬件或代码里的。3. 从零开始Farrow滤波器系数设计方法论理解了思想下一步就是如何得到那些固定的子滤波器系数c_{k,l}。这是Farrow滤波器设计的核心。设计目标通常是让整个滤波器在感兴趣的频带内其频率响应尽可能地逼近一个理想的、延迟为d的分数延迟器。理想分数延迟器的频率响应是H_ideal(e^{jω}) e^{-jωd}它具有线性相位群延迟恒为d。设计方法主要有两大类基于最小二乘误差准则和基于拉格朗日插值。这里我重点讲工程上最常用、也相对直观的最小二乘设计法。3.1 设计问题建模假设我们设计一个长度为N阶数为N-1的FIR滤波器用于近似分数延迟d。其传递函数为H(z, d) Σ_{k0}^{N-1} h[k] * z^{-k}而根据Farrow结构我们假设h[k]是d的L阶多项式h[k] Σ_{l0}^{L} c_{k,l} * d^l我们的目标是在d ∈ [0, 1]和ω ∈ [0, απ]α是带宽因子比如0.8表示利用80%的奈奎斯特带宽的范围内让实际响应H(e^{jω}, d)逼近理想响应e^{-jωd}。定义一个误差函数E(ω, d) H(e^{jω}, d) - e^{-jωd}设计目标就是找到一组系数c_{k,l}使得在定义的(ω, d)区域上误差E的某种范数如平方误差的积分最小。这是一个标准的优化问题。3.2 具体设计步骤与MATLAB示例虽然推导过程涉及积分和矩阵运算但幸运的是我们可以借助MATLAB等工具来高效完成。下面是一个手把手的步骤步骤一确定设计参数N: 主FIR滤波器长度抽头数。越大逼近精度越高计算量也越大。L: 多项式阶数。越高对参数d变化的拟合能力越强但结构也更复杂。通常L3三次是一个很好的折中。alpha: 归一化带宽0到1之间。我们只关心这个带宽内的性能。设计网格在d和ω的二维平面上需要划分密集的网格点来进行离散化优化。例如d从0到1步进0.01ω从0到alpha*pi步进pi/500。步骤二构建最小二乘问题对于每一个网格点(ω_i, d_j)我们可以写出H(e^{jω_i}, d_j) Σ_{k0}^{N-1} (Σ_{l0}^{L} c_{k,l} * d_j^l) * e^{-jω_i k}令C是一个将所有系数c_{k,l}按特定顺序例如先按k再按l排列成的列向量。那么对于所有网格点我们可以建立一个巨大的线性方程组A * C ≈ B其中A矩阵的每一行对应一个(ω_i, d_j)点其元素由d_j^l * e^{-jω_i k}构成B向量对应每个点的理想响应值e^{-jω_i d_j}。步骤三求解系数这是一个超定线性方程组我们用最小二乘法求解C (A^H * A) \ (A^H * B)这里^H表示共轭转置\是MATLAB中的左除运算符用于求解最小二乘解。步骤四验证与评估求解出系数C后将其重新排列成N x (L1)的矩阵C_mat其中第k行第l列就是c_{k,l}。 然后我们可以固定一个d值如0.5用C_mat生成对应的FIR系数h C_mat * [1; d; d^2; ... d^L]绘制其频率响应并与理想延迟e^{-jωd}比较。扫描d从0到1观察滤波器幅频响应应接近1和群延迟应接近d的变化是否平滑。实操心得在MATLAB中构建矩阵A时一定要注意索引和维度的对应关系这是最容易出错的地方。一个技巧是先用循环写一个清晰但低效的版本确保逻辑正确后再尝试用向量化操作meshgrid,kron等进行加速。另外带宽因子alpha不要设得太大如0.95以上因为逼近理想延迟器在频带边缘非常困难强求会导致通带内纹波增大。通常0.8~0.9是稳健的选择。4. 超越分数延迟Farrow结构的变体与应用扩展虽然Farrow结构最初是为可变分数延迟而生但其“用固定子滤波器组合实现参数可变”的思想可以被推广到更广泛的可变参数滤波器设计中。这才是它真正强大和有趣的地方。4.1 可变截止频率滤波器这是另一个非常实用的场景。假设我们需要一个截止频率fc可变的低通滤波器。理想情况下滤波器的系数应该是fc的函数。我们可以借鉴Farrow思想将每个系数h[k]用关于归一化截止频率ff fc / fsfs为采样率的多项式来近似h[k](f) ≈ Σ_{l0}^{L} c_{k,l} * f^l设计过程与分数延迟器类似但目标函数变了。此时我们需要让设计出的滤波器在通带[0, f]内响应接近1在阻带[fΔ, 0.5]内响应接近0Δ是过渡带。这同样可以转化为一个在(ω, f)二维区域上的最小二乘优化问题只是误差权重函数需要精心设计例如在通带和阻带赋予高权重在过渡带赋予低权重。实现结构和经典Farrow结构一模一样只是输入参数从延迟d变成了截止频率f子滤波器系数c_{k,l}是针对截止频率变化而优化的。4.2 应用于采样率转换多项式插值这是Farrow结构最早、也是最成功的应用之一。在异步采样率转换中我们需要计算输入序列在非整数采样点上的值这本质上就是一个分数延迟问题。例如常用的三次拉格朗日插值器其系数就是关于分数间隔μ相当于d的三次多项式。你可以验证这些多项式系数正好可以排列成一个4x4的C_mat矩阵从而完美地用Farrow结构实现一个高效的三次插值器。更一般地任何基于多项式的插值核如分段抛物线、样条插值都可以用Farrow结构来实现。这使得它在数字上下变频、软件无线电接收机中成为了标准模块。4.3 其他可变参数理论上只要滤波器性能是某个参数p的平滑函数并且可以用多项式较好地近似就可以尝试Farrow结构。例如可变带宽的带通滤波器中心频率固定带宽可变。可变Q值的谐振器谐振频率固定品质因数Q可变。可变滚降因子的升余弦滤波器。设计经验在将这些扩展时最大的挑战在于多项式近似的有效性。如果滤波器系数随参数p的变化非常剧烈或非线性那么可能需要很高的多项式阶数L才能较好近似这会急剧增加子滤波器的数量L1个和计算量。因此在决定采用Farrow结构前一定要先分析系数随参数变化的曲线是否相对平滑。通常在参数变化范围较小如d在0~1之间时Farrow结构表现最佳。5. 硬件实现考量与实战中的“坑”将Farrow滤波器从算法模型搬到FPGA或ASIC上是另一个充满细节的战场。这里分享几个我趟过的雷区。5.1 计算复杂度与资源优化一个N抽头、L阶多项式的Farrow滤波器需要(L1)个并行的N抽头FIR子滤波器。直接实现的乘法器数量是N*(L1)这看起来很大。但我们可以利用结构特点进行优化子滤波器合并计算观察结构每个输入采样x[n]需要与所有子滤波器的第一抽头系数c_{0,0}, c_{0,1}, ..., c_{0,L}相乘然后分别延迟、再与下一组系数相乘。我们可以将同一抽头位置k上的、属于不同子滤波器的系数c_{k,0} ... c_{k,L}视为一个向量。当计算该抽头的贡献时我们实际上是在计算这个系数向量与权重向量[1, d, d^2, ..., d^L]^T的点积。这可以转化为一个先乘累加、再乘以输入信号x[n-k]的过程有时能减少乘法器数量。多项式求值优化计算权重d^l可以使用霍纳法则将Σ c_{k,l} * d^l的计算转化为y_k c_{k,0} d*(c_{k,1} d*(c_{k,2} ... d*c_{k,L})...)这样只需要L次乘法和L次加法而不是L次幂运算和L次乘法。系数对称性利用如果目标频率响应是对称的如线性相位那么设计出的子滤波器系数c_{k,l}也可能呈现某种对称性。例如对于可变分数延迟滤波器其冲激响应关于中心点近似对称。利用这种对称性可以将乘法器数量几乎减半。5.2 动态参数更新的时序问题参数d或f是动态更新的。在硬件中这需要仔细处理更新速率d更新的速度不能超过数据处理流水线的“吞吐量”。如果d每个时钟周期都变那么权重d^l需要每个周期重新计算并且必须确保在用到新权重的时刻数据流水线中对应的是新的输入数据。通常d的更新速率远低于采样率。同步新的d值必须与输入数据流正确同步。一个常见的做法是使用一个“参数有效”信号当d更新时该信号拉高一个周期标志着从此之后进入滤波器的数据将使用新的d值。这需要在数据路径上插入相应的控制逻辑。5.3 有限字长效应与精度管理这是硬件实现中最容易出问题的地方尤其是在定点设计中。系数量化设计得到的c_{k,l}通常是高精度浮点数。我们需要将其量化为定点数如Q格式。量化会引入误差可能导致频率响应偏离设计目标甚至不稳定。必须进行充分的仿真扫描不同的量化位宽如12位、16位、18位观察通带纹波、阻带衰减等关键指标的变化找到满足性能要求的最小位宽。中间结果位宽扩展在计算d^l以及子滤波器乘累加的过程中中间结果的动态范围会扩大。例如计算d^2假设d是16位小数可能需要32位来保证精度不损失。在加法树中位宽扩展更明显。必须为每条数据路径仔细规划位宽防止溢出和精度过度损失。一个安全的方法是先做全精度仿真记录中间结果的最大最小值再据此确定定点位宽。权重计算的非线性d^l的计算尤其是高次幂对d的量化误差非常敏感。当d接近0或1时d^l的值可能非常小量化误差会占据主导。可以考虑采用查找表LUT来存储d^l的值或者使用分段线性近似等方法来计算高次幂以平衡精度和资源消耗。踩坑实录在一次通信接收机项目中我们使用Farrow结构做定时恢复。仿真时性能完美但上板后误码率总是差一点。用逻辑分析仪抓取中间数据发现当分数延迟d在0.1以下时输出信号有明显的失真。最终定位到问题计算d^3的模块为了节省乘法器我们采用了(d*d)*d的顺序计算。由于d很小d*d的结果在定点量化后直接变成了0导致三次项完全失效。教训是对于可能接近零的小数运算必须保证中间结果的精度或者改变计算顺序如先计算高精度浮点再量化或者采用查找表。后来我们改用一个小型的、针对d的LUT来存储d^2和d^3的值问题得以解决。6. 性能评估与设计实例一个可变延迟线让我们用一个具体的设计实例把前面所有的点串起来。目标设计一个用于音频处理的可变分数延迟线要求延迟d在0到1个采样周期内连续可调音频带宽为20kHz采样率fs48kHz因此alpha ≈ 0.833。设计参数选择主滤波器长度N32。这是一个折中能提供较好的阻带抑制。多项式阶数L3。三次多项式足以在0-1区间平滑近似。带宽因子alpha0.85略高于需求留有余量。采用最小二乘准则在d[0:0.02:1]和ω[0:π/1000:0.85π]的网格上设计。设计结果 使用MATLAB的farrow设计函数或按前述步骤自编代码得到32x4的系数矩阵C_mat。性能评估固定延迟响应取d0.5生成系数h C_mat * [1; 0.5; 0.25; 0.125]。绘制其频率响应。实测在0-20kHz通带内幅度起伏 0.01 dB群延迟波动 0.001个样本非常接近理想的0.5样本延迟。动态扫描让d从0线性变化到1观察滤波器幅频响应。在整个过程中通带内幅度保持平坦阻带20kHz衰减始终大于80dB。群延迟值紧密跟随d值变化误差在±0.005个样本以内。听感测试用一段正弦波扫频信号通过该延迟线同时用低频三角波调制d在0-1之间变化。输出信号听起来应该是纯净的、音高无变化因为只有延迟无频移但带有周期性相位调制的效果。如果设计不良可能会引入可闻的谐波失真或振幅调制。本例中听感干净。硬件实现概要FPGA系数存储32*4128个系数每个量化为18位有符号定点数存储在ROM中。数据处理路径输入音频数据x[n]为24位。设计4个并行的32抽头FIR滤波器H0-H3。利用系数的近似对称性每个FIR采用转置结构乘法器数量减半约需(32/2)*4 64个乘法器。参数d为16位无符号小数0x0000~0xFFFF 对应 0~1。权重计算模块使用一个时钟周期计算d2 d*d下一个周期计算d3 d2*d均为32位中间结果。同时使用一个三级流水线的霍纳法则单元计算每个抽头k的加权和sum_k c_k0 d*(c_k1 d*(c_k2 d*c_k3))。最终输出y[n]为24位由所有sum_k * x[n-k]的累加和得到。资源与性能在Xilinx Artix-7 FPGA上占用约70个DSP slice最高运行时钟频率可达150MHz远高于48kHz的音频采样率有充足资源进行多通道处理。通过这个实例可以看到一个精心设计的Farrow滤波器能够在参数连续变化时保持卓越且稳定的性能同时满足实时硬件处理的要求。它不再是教科书上一个抽象的数学结构而是一个能解决实际工程难题的强力工具。