1.算法描述
目前的通用做法就是采用多项式插值滤波器去实现一些分数倍比较大的采样率转换。同时采用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核心程序
N = 23; % filter order, odd better
L = N+1; % filter length;
Npt = 256; % no. of frequency points for plots
w = (0:1:Npt-1)/Npt; % frequenc scan (0,1)
delay = [0 0.1 0.2 0.3 0.4 0.5]; % delay range x=0..0.5
Nfil = length(delay); % number of filters
h = zeros(1,L); % impulse response vector
hvec=zeros(Nfil,L); % impulse response coefficient matrix
magresp = zeros(Nfil,Npt);
phasdel = zeros(Nfil,Npt-1);
xvec=zeros(Nfil,1); % fractional delay vector
P = 2; % polynomial order for FARROW structure (ca. 1-5)
C=zeros(P+1,N+1); % polynomial coeff. matrix
wp = 0.8; % normalized bandwidth (0-1.0)
for i=1:Nfil
...............................................................
end
for k=1:N+1
cc=polyfit(xvec,hvec(:,k),P); % fit P:th-order polynomial to each coeff set
C(:,k)=cc';
end
for j=1:Nfil
d=delay(j);
if d==0
d=d+0.0000001; % add 0.001 to avoid sin(0)/0;
end
h = C(P+1,:); % coeffs. via pol. approximation
for n=1:P
h=h+d^n*C(P+1-n,:);
end
h=h/sum(h); % scale response at zero freq. to unity
H = freqz(h,1,w*pi);
magresp(j,:) = abs(H);
uwphase=-unwrap(angle(H));
phasdel(j,:) = uwphase(2:Npt)./(w(2:Npt).*pi); % avoid divide by zero
end