1.算法描述
采样速率转换（SRC）在通信中非常普遍。一般有两种方法：一种是通过D/A重构信号，再采样，从而实现采样速率的转换；另一种是利用数字滤波器直接进行采样转换。数字滤波器有CIC，多相，FARROW。
在一个采样间隔T内，FARROW滤波器结构的系数不变，可变的是分数延迟，在一个采样间隔内，可任意改变分数延迟，提高采样率。一般拉格朗日内插采用直接的FIR滤波器实现，虽然简单，但系数随迟延变化而变化，一旦延迟D变化，需重新计算或保存滤波器系数，且滤波器系数与内插因子成比例，实现不灵活。因此提出一种更为灵活高效的实现结构，FARROW滤波器结构。FARROW滤波器结构不仅可应对变化的分数延迟而且大大减少运算，提高运行效率。理想的分数延迟滤波器是since滤波器，在时域里，它是一个无限长冲击响应滤波器，不可实现。
目前的通用做法就是采用多项式插值滤波器去实现一些分数倍比较大的采样率转换。同时采用farrow结构实现更简便，即采用farrow类型滤波器来实现更为简便，farrow类型滤波器也称抽取滤波器。一般的数学模型为重采样模型为，采样信号x(mts)经过内插器h(t)，输出信号：在时刻t＝kti对信号进行重采样，输出信号：假设h(t)是特定的脉冲响应，这里的目的是计算t＝kti时刻y(kti)的采样值，因此首先需要定义x(mts)的采样基准时刻mkts，这个时刻刚好在t＝kti时刻之前，因此其中int[z]表示不大于z的最大整数，mk为插值基点，决定输入序列中参与运算的采样点，由插值时刻t＝kti决定。因此插值时间t＝kti可表示为mkts加上一个正小数部分的求和形式：代入得到最后的结果：其中参数uk为误差间隔，决定内插滤波器冲击响应系数，其范围为uk∈[0,1)。插值基点mk和误差间隔uk表示了ts和t之间的关系，如图2所示。而对于多项式插值滤波器来讲，得到通用的插值公式如下：其中lk(uk)的如果采用4阶拉格朗日插值得到的结果为：和得到转化后的结果：式中的a1、a2、a3、a4分别为：和a4＝x(mk+1)。最终转换成farrow类型滤波器实现方式：y(n)＝((a1uk+a2)uk+a3)uk)+a4。其中a1、a2、a3、a4可以分别写成滤波器组的形式。滤波器系数coeff分别为{-1/6、1/2、-1/2、1/6}，{1/2、-1、1/2、0}，{-1/3、-1/2、1、-1/6}，{0、1、0、0}。其中插值基点mk和误差间隔uk的计算过程在前面已经有描述过。而在fpga实现的过程当中主要以定点类型的数据进行运算，需要将正常的浮点型数据量化成整形数据。在上述公式中，需要量化的系数主要是滤波器系数coeff和误差间隔uk。对于farrow类型滤波器来说，需要选择合适的量化位宽，保证对于频谱的幅频响应满足，同时也要防止滤波器不要溢出。
farrow滤波器是一种连续可变时延的分数时延滤波器，这种滤波器的结构是由farrowcw于1988年提出，起初是用来解决声纳学中的分数时延问题。普通数字延时滤波器虽然结构简单，但系数计算过程复杂，在延时参数快速变化时，系数更新速度无法满足实时性要求，在工程应用上受限制。采用farrow结构数字延时滤波器能够更加灵活高效地进行分数延时滤波，延时参数改变时，无需重新计算滤波器系数，更容易在现场可编程门阵(fpga)上实现。信号处理的fpga实现过程中，往往需要大量消耗的乘法资源，从而导致fpga的乘法器资源成为系统瓶颈，本设计介绍了一种基于fpga的farrow滤波器设计方法，该方法采用对称结构的滤波器求解方法，充分利用乘法资源，高效实现farrow滤波器功能。————————————————
farrow滤波器的基本结构如下：
Farrow滤波器的原理：
其中,μm表示当前输出抽样点与其前面输入抽样点之间的距离,且有0≤μm<1。由(5)式即可得出一种实现SRC滤波的多项式滤波器,一般称为Farrow结构,该结构的联络图如图3所示。对于距离原来抽样位置为μm的任何输出抽样值,若用t=μm代入所在位置的分段多项式就可以计算出来,而不需要存储这些抽样值。
2.仿真效果预览
matlab2022a仿真结果如下：
3.MATLAB核心程序
x = zeros(1,L);
x(1) = 1;
% set the initial conditions
x0 = zeros(1,size(P,1)-1);
t0 = 0;
if 1,
[h t N t0 x0] = asrc_farrow(x,1/M,t0,x0,P);
else,
[h t N t0 x0] = asrc_farrow_loop(x,1/M,t0,x0,P);
end;
% plot the approximated impulse response and overlay the filter output
subplot(5,1,3);
plot(hh);
hold on;
plot(h,'k');
hold off;
% plot the magnitude response
[H,f] = freqz(h/M,1,M*np);
subplot(5,1,4);
% first drop a marker at the original nyquist
plot([pi/M pi/M],[FLOOR 0],'k');
hold on;
% plot dB magnitude with a -60dB floor
plot(f,max(20*log10(abs(H)),FLOOR));
% zoom in on the passband
%plot(f(1:np)*M,max(20*log10(abs(H(1:np))),-60/10)*10,'r');
hold off;
% plot the phase response
subplot(5,1,5);
plot(f,angle(H));
hold on;
% zoom in on the passband
plot(f(1:np)*M,angle(H(1:np)),'r');
hold off;
end;