Farrow结构分数时延滤波器：系数构造与定时同步环路实现详解

简介：基于Farrow滤波器结构的时间同步算法MATLAB仿真，面向通信、声纳等领域需要处理分数时延和符号时间同步的工程师与研究人员，运行环境为MATLAB 2021a，适合算法验证与课程设计参考。资源包共6个文件，以4个.m脚本为主，涵盖主运行程序、Farrow滤波器系数构造、异步采样率转换等核心函数模块，另附1张结果截图和1段avi格式操作录像，压缩包约1.94MB，便于快速下载学习。Farrow滤波器是1988年提出的连续可变时延分数时延滤波器，可避免普通数字延时滤波器在延时参数快速变化时系数更新不及时、工程应用受限的问题。仿真代码完整展示时间同步算法与Farrow滤波器结构结合的实现流程，录像演示了MATLAB当前文件夹路径设置和运行步骤，即使初次接触也能快速复现实验并理解关键细节。已有726人学习下载。

1. 时间同步的分数时延瓶颈：为什么 Farrow 结构上了台面

解调链路里最先暴露工程味的往往不是载波同步，而是定时同步。上采样后的信号波形看着人畜无害，采样时刻偏偏落在符号间隔的 0.3 个码片上，星座点就被吹散成一团。整数倍的采样偏移靠找峰、搬移就能修掉，真正麻烦的是零点几个采样周期的分数时延。普通数字延时滤波器每改一次延时就要重新算一遍抽头系数，延时参数快速变化时系数更新跟不上，这在声纳、突发通信和雷达信号处理 matlab 仿真里都是老问题。1988 年 Farrow 提出的滤波器结构把时延参数从滤波器系数里拆出来，让分数时延变成一个多项式求值问题，后来成为符号同步和采样率转换的标准件。这篇博文把系数构造、环路实现、参数调优到工程化预计算完整拆一遍，读者对象是正在做时间同步算法仿真、或准备把 Farrow 结构搬进实时链路的工程师和研究生。

2. Farrow 滤波器系数构造：从 Lagrange 逼近到 construct_farrow.m

2.1 为什么不用一次一算的普通延时滤波器

数字延时滤波器的原理不复杂：要延迟 D 个采样周期，滤波器的理想冲激响应是h[n] = sinc(n - D)，再加窗截断。D 取 3.37 和取 3.71 是两套完全不同的系数，每更新一次 μ 就得重新生成 N 个抽头、重新做一次卷积。定时环路每符号要更新几千次，这个计算量在实时系统里是穿肠毒药。

Farrow 结构换了个思路：把冲激响应写成分数时延 μ 的 P 阶多项式，即 h_μ[n] = Σ_{p=0}^{P} c_p[n]·μ^p。系数 c_p[n] 与 μ 无关，仿真前算一次存好；运行时输入一个 μ，做 P 次乘加就得到新延时下的插值结果。μ 变，系数不动，这正是时间同步算法里最需要的性质。Farrow 最早就是为声纳学中的分数时延问题提出的，后来被通信接收机借用来做定时恢复，结构本身没有变化，变的是外部环路的组织方式。

2.2 用最小二乘拟合构造子滤波器系数

常见做法是取一组离散的归一化延时值，对每个延时值构造理想的分数延时响应，再用多项式拟合出 c_p[n] 矩阵。下面是construct_farrow.m的核心思路，参数化之后可以直接改：

function C = construct_farrow(N, P) % N: 子滤波器抽头数，P: 多项式阶数 d = linspace(-0.5, 0.5, 4 * (P + 1)); % 一组分数延时采样点 H = zeros(length(d), N); for k = 1:length(d) n = (0:N-1) - (N/2 - 1) - d(k); H(k, :) = sinc(n); % 理想分数延时响应 end V = ones(length(d), P + 1); for p = 1:P V(:, p + 1) = d(:) .^ p; end C = V \ H; % 最小二乘拟合 end

第一段循环构造的是不同分数延时下的理想响应，第二段用 Vandermonde 矩阵做多项式拟合，最后得到的 C 尺寸是(P+1) × N，第 p 行就是 μ 的 p 次幂对应的子滤波器系数。使用上有个细节：拟合点 d 的分布越均匀，拟合误差越平，但 V 矩阵在 P 较大时条件数会变差。实际里 P 取 4 或 5 就够用，继续加阶数对群延迟边缘的改善有限，反而容易在带外引入吉布斯抖动。

2.3 关键参数怎么选：N、P 与 μ 的边界

还有两个容易被忽略的坑。第一，Farrow 结构本身带固定群延迟 N/2 - 1 个采样点，串进定时环路后要把这部分从符号索引里减掉，否则锁定的符号位置永远偏一个固定量。第二，插值误差在信号带边缘最大，成型滤波器的滚降系数最好留 0.2 以上余量，别把 Farrow 当万能补偿器去扛陡峭的带外衰减。

3. 可变时延插值与时序环路：asrc_farrow 到 asrc_farrow_loop

3.1 工程文件的职责划分

这套仿真里三个 m 文件的边界很清楚：construct_farrow.m只生成系数矩阵 C，asrc_farrow.m负责对一段缓冲做单点分数时延插值，asrc_farrow_loop.m把插值器塞进定时误差检测环路里，逐符号输出同步后的样值。Runme.m是顶层脚本，把调制、成型、加噪声、以及上面三个模块串起来跑完整个仿真。

3.2 单点插值：霍纳法则求多项式

asrc_farrow.m对输入的 N 点缓冲做一次插值，输出目标时刻的样值。多项式求值用霍纳法则，避免直接算 μ 的幂次，数值稳定性好一些：

function y = asrc_farrow(buf, mu, C) % buf: 以目标采样时刻为中心的一段 N 点缓冲 % mu : 归一化分数时延，范围 [0, 1) P = size(C, 1) - 1; y = 0; for p = P:-1:0 y = y * mu + (C(p + 1, :) * buf(:)); % 内积是标量 end end

这段代码里C(p+1,:) * buf(:)计算的是第 p 阶子滤波器与缓冲的卷积结果（只取中心点），然后从最高次向下逐层乘 μ 累加。外面的定时环路每一拍都会重新给一个 μ，因此这里不能把 C 与 μ 的乘积缓存，否则就失去了 Farrow 结构随 μ 连续变化的意义。缓冲的取法需要注意：仿真里用round(k*sps + phase)定位整数采样点，再左右各取 N/2 个点；如果缓冲跨越序列边界，要么丢弃首尾若干符号，要么做边缘延拓，不要硬取越界索引。

3.3 环路实现：Gardner 定时误差检测加一阶相位更新

asrc_farrow_loop.m是整个时间同步算法的核心。它用 Gardner 算法做定时误差检测，公式是err = Re{x_mid · (x_k - x_{k-1})*}，其中 x_k 是当前符号时刻的插值样值，x_mid 是符号间隔中点处的插值样值。Gardner 的好处是对载波相位不敏感，可以在载波同步之前先做定时同步：

function sym_out = asrc_farrow_loop(rx, sps, C, Kp) % rx : 过采样接收序列 % sps : 每符号采样数，仿真里取 4 % Kp : 环路增益，常见取 0.005 ~ 0.05 Nf = size(C, 2); mu = 0; % 分数时延初值 phase = 0; % 定时相位，以采样周期为单位 nsym = floor(length(rx) / sps); sym_out = zeros(nsym, 1); x_prev = 0; for k = 1:nsym-1 n0 = round(k * sps + phase); % 当前符号采样参考点 nMid = n0 + round(sps / 2); % Gardner 中间点 if n0 - Nf/2 < 1 || nMid + Nf/2 > length(rx) break; end buf = rx(n0 - Nf/2 + 1 : n0 + Nf/2); bufM = rx(nMid - Nf/2 + 1 : nMid + Nf/2); xk = asrc_farrow(buf, mu, C); % 符号时刻插值 xmid = asrc_farrow(bufM, mu, C); % 中间时刻插值 err = real(xmid) * (real(xk) - real(x_prev)); % Gardner TED phase = phase + Kp * err; if phase > sps/2, phase = phase - sps; end % 相位回绕 if phase < -sps/2, phase = phase + sps; end x_prev = xk; sym_out(k) = xk; end end

相位累积和回绕是这段代码最容易出错的地方。phase 累加的是采样周期为单位的小数偏差，超出 ±sps/2 就要模回，否则整数采样点 n0 会一直朝一个方向漂移，最终跑出缓冲区。环路增益 Kp 不能贪大，大会让相位抖动直接反映到星座图上；也不能太小，太小收敛需要上千个符号。初值 μ 通常设 0，因为 Gardner 环路的捕获范围在 ±0.5 个符号周期内，仿真里插入的真实延时只要不超出这个范围，环路就能自己拉过去。

4. 时间同步仿真参数调优：环路增益与环路带宽怎么定

4.1 一阶环路还是二阶环路

上面的代码是一阶环路，只有比例支路 Kp。它能锁定恒定的定时相位偏差，但如果收发两端存在采样钟频偏，相位会随时间线性变化，一阶环路会留下一个固定的稳态跟踪误差。工程上通常给环路加一个积分支路，变成二阶环路：比例支路负责快速拉偏，积分支路负责消除稳态误差。简化实现是维护两个累加器：

alpha = 0.01; % 比例增益 beta = 0.0004; % 积分增益，取 alpha 的 1/25 上下 integr = 0; ... err = real(xmid) * (real(xk) - real(x_prev)); integr = integr + beta * err; phase = phase + alpha * err + integr;

alpha 和 beta 的比例直接决定环路的阻尼。beta 太大，相位会围绕锁定点震荡，收敛曲线出现明显过冲；beta 太小，遇到采样钟频偏时误差收敛得很慢。常见做法是先定 alpha，再取beta = alpha^2 / 4附近作为临界阻尼参考值，然后在仿真里做小范围扫描。

4.2 顶层仿真脚本 Runme.m 的参数组织

Runme.m里需要把调制、成型、定时偏移、噪声和 Farrow 环路串起来。一个可复现的结构如下：

clear; close all; rng(0); M = 4; % QPSK sps = 4; % 每符号采样数 nsym = 2048; data = randi([0 1], nsym*log2(M), 1); sym_tx = pskmod(bi2de(reshape(data, log2(M), []).', 'left-msb'), M, pi/4); up = upsample(sym_tx, sps); h = rcosdesign(0.35, 6, sps, 'sqrt'); % 平方根升余弦成型 tx = filter(h, 1, up); delay_sym = 0.37; % 插入真实定时偏差 rx = filter([zeros(1, round(delay_sym*sps)) 1], 1, tx); rx = awgn(rx, 25, 'measured'); C = construct_farrow(8, 4); sym_rx = asrc_farrow_loop(rx, sps, C, 0.01); scatterplot(sym_rx(200:end));

这个脚本里定时偏差用delay_sym = 0.37表示 0.37 个符号周期，先做整数部分搬移，再交给 Farrow 环路的分数部分。成型滤波器用rcosdesign(0.35, 6, sps, 'sqrt')，滚降 0.35，跨越 6 个符号，保证信号带宽不会顶到 Farrow 插值的带外误差区。信噪比设 25 dB 是为了让星座图散点直观可读，太低的话 Gardner 误差信号的信噪比也变差，环路收敛曲线会抖得不好看。

4.3 参数到现象的映射：调参时看什么

一个很有用的验证手段是画相位收敛曲线：把循环里的 phase 存下来，看它是否从 0 平滑趋向delay_sym * sps的余数附近。如果曲线发散先检查回绕逻辑，如果收敛但星座图带有规律性偏转，检查固定群延迟 N/2 - 1 有没有补偿。失败时不要直接怀疑 Farrow 系数，先把环路断开，用真实的 delay_sym 直接插值做对照，这样能区分是系数问题还是环路问题。

5. 跑通验证与 μ 查找表预计算：把 Farrow 仿真搬进实时链路

5.1 环境对齐与性能验证

这套仿真在 MATLAB 2021a 下运行，打开Runme.m后第一件事是检查 MATLAB 左侧当前文件夹路径，必须切到程序所在目录，否则Runme.m调用不到同目录下的construct_farrow.m。配套的操作录像 0042.avi 用 Windows Media Player 播放，对照录像里查看文件夹路径位置和运行顺序，能省掉大部分环境报错时间。运行结束后除了散点图，建议补一段 EVM 统计：

n_valid = 200:length(sym_rx)-50; evm = sqrt(mean(abs(sym_rx(n_valid) - sym_tx(n_valid)).^2) ... / mean(abs(sym_tx(n_valid)).^2)) * 100; fprintf('EVM = %.2f%%\n', evm);

25 dB 信噪比下，同步做得好的 QPSK 链路 EVM 通常能压在 12% 以内。如果 EVM 偏高，优先检查丢弃符号数，循环开头和结尾的瞬态样本不能参与统计。

5.2 μ 查找表预计算：去掉实时多项式求值

Farrow 结构虽然把系数和 μ 解耦了，但每符号仍要做 P+1 次内积。在 FPGA 或嵌入式实现里，更常见的做法是把 μ 量化成 M 个格点，提前把每个格点对应的 N 个插值系数算好存进查找表，运行时查表选系数，连霍纳循环都省掉。仿真阶段可以先用连续 μ 版本做黄金模型，再看量化损失：

Mq = 64; % μ 量化级数 LUT = zeros(Nf, Mq); for q = 1:Mq muq = (q - 0.5) / Mq; for n = 1:Nf LUT(n, q) = sum(C(:, n) .* muq.^(0:P).'); end end

LUT 的每一列对应一个量化后的 μ，查表时用q = round(mu * Mq)取列。量化级数取 64 时相邻 μ 间距只有 1/64，带来的 EVM 损失通常可以压到 0.1 dB 以内。这样做还有个附带好处：系数变成常数，定点化时可以把整数部分和小数部分分开存储，再用移位累加完成乘法，时序收敛压力小很多。验证时用同一段接收序列分别跑连续版本和查表版本，对比两者的星座图散点分布，只要量化损失在可接受范围，就可以放心把查表版本提交到硬件实现。

本文还有配套的精品资源，点击获取