Using Farrow filter on the basis of piecewise cubic polynomial interpolation for digital signal resampling

Contents
Introduction
In the previous section we considered Farrow filter for digital signal resampling. The structure of Farrow digital resampling filter on the basis of piecewise cubic polynomial interpolation (Figure 1) is got.

Figure 1. Functional chart of optimized Farrow filter
on the basis of piecewise cubic polynomial interpolation.
Calculating coefficients of an interpolation polynom requires one multiplication by $\frac{1}{6}$ and two trivial multiplication by $\frac{1}{2}$.
In this section we will consider examples of using the resampling filter to compensate fractional delay, digital interpolation and fractional resampling $\frac{P}{Q}$ times.
Recalculating sample indexes of an input signal in case of piecewise cubic interpolation
High performance of a digital resampler is reached due to using piecewise cubic interpolation and recalculating sample indexes of the input signal $s(n)$, $n =0,1,2, \ldots$ in the interval $[-2, \ldots 1]$, as it is shown in Figure 2.
Figure 2. Recalculating sample indexes of an input signal in the interval $[-2, \ldots 1].$
Let the sampling frequency of the input signal $s(n)$, $n =0,1,2, \ldots$ be equal to $F_\textrm{s}$. Then the signal $y(k)$, $k =0,1,2, \ldots$ at the resampler filter output will have the following sampling frequency $F_\textrm{y} = \frac{P}{Q} \cdot F_\textrm{s}$. Besides the first sample of the output signal $y(k)$ is delayed concerning the input one $s(n)$ by the value $0 \leq x_0 <1$.
Here and elsewhere in this section as well as in the previous section the variable $n$ indexes samples of the input signal $s(n)$, and the variable $k$ indexes samples of the output signal $y(k)$.
As sampling frequencies of the input and output signals are various, the same indexes of the input and output samples correspond to different timepoints. For example, the $5-$th sample of the output signal $y(5)$ after resampling does not correspond to the $5-$th sample of the input signal $s(5)$. Therefore we have to specify what scale we use to index samples. We will use three scales: the scale $n$ is for timepoints concerning samples of the input signal $s(n)$, the scale $k$ is for timepoints concerning the output signal and the scale $t$ in the interval $[-2 \ldots 1]$ is for calculating coefficients of a cubic polynom.
For the $k$ sample of the output signal $y(k)$ according to the scale $k$ it is necessary to calculate the corresponding indexes $n$ of the input signal $s(n)$ according to the scale $n$ (Figure 2a), and then to pass to the scale $t$ in the interval $[-2 \ldots 1]$ for calculating coefficients of a cubic polynom (Figure 2b). At that the output sample point $k$ according to the scale $n$ has to be recalculated under $0 \leq \Delta_k < 1$ according to the scale $t$.
The point $x_k$ of the $k$ output sample after resampling according to the scale $n$ of the input signal is equal to:
$$x_k = k \cdot \frac{Q}{P} - x_0.$$
(1)
To calculate a cubic polynom we need four samples of the input signal $s(n)$, $s(n-1)$,$s(n-2)$ and $s(n-3)$, as it is shown in Figure 2a. At that two samples $s(n)$ and $s(n-1)$ have to be more to the right of $x_k$ and two other samples $s(n-2)$ and $s(n-3)$ have to be more to the left of it.
Then the index $n$ of the sample $s(n)$ which corresponds to the timepoint $t = 1$ according to the scale $t$ from -2 to 1 is equal to:
$$n = \lfloor x_k \rfloor + 2,$$
(2)
where the functional $\lfloor x_k \rfloor$ means rounding to floor.
The value $\Delta_k$ according to the scale $t$ is equal to:
$$\Delta_k = \lfloor x_k \rfloor + 1 - x_k.$$
(3)
In Figure 3 the example of recalculating indexes of the input and output signals for the $5-$th output sample $y(5)$ under $P = 4$, $Q = 3$, $x_0 = 0.2$ is shown.
Figure 3. Example of recalculating indexes under $P = 4$, $Q = 3$, $x_0 = 0.2$
The point $x_5$ according to the scale $n$ in accordance with (1) is equal to:
$$x_5 = 5 \cdot \frac{3}{4} - 0.2 = 3.55.$$
(4)
Then the index $n$ in accordance with (2) is equal to:
$$n = \lfloor 3.55 \rfloor + 2 = 5,$$
(5)
and the parametric variable $\Delta_5$ according to the scale $t$ is equal to:
$$\Delta_5 = \lfloor 3.55 \rfloor + 1 - 3.55 = 0.45.$$
(6)
Thus, it is necessary to calculate coefficients of the cubic polynom $P(t) = a_0 + a_1 \cdot t +a_2 \cdot t^2 + a_3 \cdot t^3$ according to the following equations:
\begin{equation*} \begin{cases} a_0 = s(n-1),\\ a_3= \frac{1}{6} \cdot \left( s(n) - s(n-3)\right) + \frac{1}{2} \cdot \left( s(n-2) - s(n-1)\right), \\ a_1 = \frac{1}{2} \cdot \left( s(n) - s(n-2) \right)-a_3,\\ a_2 = s(n) - s(n-1) - a_1 - a_3,\\ \end{cases} \end{equation*}
(7)
then the value $y(k)$ can be got as the interpolation polynom value $P(t)$ for $t=-\Delta_k$, i.e. $y(k) = P \left(-\Delta_k \right)$.
Using Farrow filter on the basis of polynomial interpolation to compensate the signal fractional delay
Let the input signal $s(n)$ contain $N=8$ samples as it is shown in Figure 4.

Figure 4. Input signal of Farrow filter.
Compensate the fractional time delay of this signal by the value $x_0 = 0.25$ of the sampling period. Upon changing the fractional delay the sampling frequency is unchanged, the parametric variables $P$ and $Q$ are equal to $P = Q = 1$ and the number of samples at the Farrow filter output is equal to the number of samples at the input, i.e. $N=8$.
The calculation process of the signal fractional delay $s(n),$ $n = 0 \ldots 7,$ is given in Table 1. This process is shown in Figure 4.
Table 1. Equalization of the signal fractional delay $s(n)$.
The values $x_k,$ $n$ and $\Delta_k$ are calculated according to (1), (2) and (3) respectively. Then the values $s(n-3) \ldots s(n)$ for each value $k$ are given. Zero values $s(n-3) \ldots s(n)$ for $n<3$ and $n>N-1$ which fall outside the limits of the signal indexation $s(n),$ $n = 0 \ldots N-1$ are red. In the following columns of Table 1 the coefficient values of the cubic polynom $a_0 \ldots a_3$ calculated according to (7) for the current point $k$ are given. At last, the signal value $y(k) = P(-\Delta_k)$ at the Farrow filter output is given in the last column.
The pickup signal $s(n),$ $n = 0 \ldots 7$ (black) and the signal $y(k),$ $n = 0 \ldots 7$ after equalizing the fractional delay $x_0$ (red) are shown in Figure 5.

Figure 5. Input signal $s(n)$ of Farrow filter and the signal after equalizing the fractional delay $y(k)$.
Table 1 shows that upon equalizing the fractional delay the parametric variables $\Delta_k = x_0$ for all $k = 0 \ldots N-1$. In this case Farrow filter can be interpreted as the third order FIR filter with the variable group delay $\tau(\omega)$ by specifying the required value $x_0$.
The family of Farrow filter impulse response for a various fractional delay value $x_0$ is shown in Figure 6.

Figure 6. Farrow filter impulse responses for a various fractional delay value $x_0$
The magnitude $\left| H \left(\textrm{e}^{j\cdot \omega} \right) \right|^2$ and the group delay $\tau(\omega)$ of received filters are shown in Figure 7.

Figure 7. Magnitude and group delay of Farrow filters for a various fractional delay value $x_0$
The magnitude and the group delay charts allow to draw a conclusion that delay equalization is possible only within the scope of up to $0.4\pi$ rad/s. For a wider scope Farrow filter equalization on the basis of polynomial piecewise and cubic interpolation does not suit because of unacceptably high distortion of the filter characteristics. If it is required to use fractional delay equalization in a wider scope, it is necessary to use more difficult filters which use in their turn polynomial interpolation of a higher order or alternative FIR filters, or IIR filters of group delay equalization of signals.
The program implementing calculation of Table 1 and the data for plotting Figure 5 is written using DSPL library. The program source code can be received in the section of DSPL library examples.
GNU Octave and Matlab users can calculate the table 1 by script resample_lagrange_fd.m.

clear all;
close all;
clc;

p = 1;
q = 1;
x0 = 0.25;

s0 = [1 2 2 1 -0.5 -1 -2 -0.5];

tin = 0:length(s0)-1;

y = zeros(1, floor(length(s0)*q/p));

s = [0, 0, s0, 0, 0];

fprintf('---------------------------------------------------------');
fprintf('-------------------------------------------------------\n');
fprintf('k     x_k    n        d    s(n-3)    s(n-2)   s(n-1)    s(n)     ');
fprintf('a_0       a_1       a_2       a_3       y(k)\n');
fprintf('---------------------------------------------------------');
fprintf('-------------------------------------------------------\n');
for k = 0:length(y)-1
x =  k*q/p - x0       ;
n =  floor(x) + 2 + 3 ;
d =  floor(x) + 1 - x ;

tout(k+1) = x;

a0 = s(n-1);
a3 = 1/6 * (s(n) - s(n-3)) + 0.5*(s(n-2) - s(n-1));
a1 = 0.5 * (s(n) - s(n-2)) - a3;
a2 = s(n) - s(n-1) - a3 - a1;

y(k+1) = a0 - a1 * d + a2*d^2 - a3*d^3;

fprintf('%d  %7.2f   %d   %7.2f', k, x, n, d);
fprintf('%8.1f %8.1f %8.1f %8.1f ', s(n-3), s(n-2), s(n-1), s(n));
fprintf('%9.3f %9.3f %9.3f %9.3f ', a0, a1,a2,a3);
fprintf('%9.3f\n', y(k+1));

end

fprintf('---------------------------------------------------------');
fprintf('-------------------------------------------------------\n');

figure; stem(tin, s0);
hold on;
stem(tout, y, 'r');
axis([-0.5, 7.5, -2.5, 2.5]);
grid on;


Calculating family of impulse responses, magnitude and group delay of Farrow filters for a various fractional delay value $x_0$ which is shown in Figures 6 and 7 was also performed using DSPL library. The program source code of calculating magnitude and group delay of Farrow filters can also be received in the corresponding section of DSPL library examples.
GNU Octave and Matlab users can calculate family of impulse responses, magnitude and group delay of Farrow filters for a various fractional delay value $x_0$ which is shown in Figures 6 and 7 by script resample_lagrange_fd_filter.m.

clear all;
close all;
clc;

addpath('functions');

% Check Matlab or GNU Octave
isOctave = exist('OCTAVE_VERSION', 'builtin');
if(isOctave)
pkg load signal;
end

x0 = 0:0.1:0.9;

s = [0, 1, 0, 0];

figure(1); hold on;
figure(2); hold on;
figure(3); hold on;

w  = 0:0.001:pi;

for k = 1:length(x0)
h = resample_lagrange(s, 1, 1, x0(k));
figure(1); subplot(2,5,k); stem(0:3, h);
axis([0, 3, -0.2, 1.2]); grid on;

H = freqz(h, 1, w);
tau = grpdelay(h,1,w);

figure(2); plot(w/pi, 20*log10(abs(H)));
axis([0, 1, -20, 2]); grid on;
xlabel('w/pi');
ylabel('|H(e^{jw})|^2, dB');

figure(3); plot(w/pi, tau);
axis([0, 1, 0, 3]); grid on;
xlabel('w/pi');
ylabel('\tau(w)');
end

We need to use function resample_lagrange.m for the script resample_lagrange_fd_filter.m and other scripts bellow.

function y = resample_lagrange(s, p, q, x0)
% y = resample_lagrange(s, p, q, x0)
% Digital resampling by polynomial Lagrange interpolation.
% Function changes input signal s samplerate to p/q times and adds fractional
% delay.
%
% Input parameters
%  s   - input signal vector [N x 1];
%  p   - p paramter of samplarate conversion
%  q   - q paramter of samplarate conversion
%  x0  - fractional delay
%
% Ouptut parameters
%  y   - Resampled signal
%
% Author: Sergey Bakhurin (dsplib.org)
%

if(p>1 && q==1)
y = zeros(1, floor((length(s)-1)*p/q)+1);
else
y = zeros(1, floor(length(s)*p/q));
end

t = zeros(size(y));

s = [0, 0, s, 0, 0];

for k = 0:length(y)-1
x = k*q/p - x0;

t(k+1) = x;

n = floor(x) + 5;
d = floor(x)+1-x;

a0 = s(n-1);
a3 = 1/6 * (s(n) - s(n-3)) + 0.5*(s(n-2) - s(n-1));
a1 = 0.5 * (s(n) - s(n-2)) - a3;
a2 = s(n) - s(n-1) - a3 - a1;

y(k+1) = a0 - a1 * d + a2*d^2 - a3*d^3;

end

Using Farrow filter as a digital signal interpolator
Consider an example of using the digital resampling filter for digital signal interpolation.
Let the input signal $s(n),$ $n = 0 \ldots N-1,$ contain $N=8$ signal samples as it is shown in Figure 8a. We have already used this signal as an input signal of the digital fractional delay equalization filter.
It is necessary to perform digital signal interpolation and to increase the signal sampling frequency $y(k),$ $k = 0\ldots K-1$ at the Farrow filter output $P=10$ times.
In case of digital interpolation the number of signal samples at the output will be equal to (see Figure 8b):
\begin{equation*} K = (N-1) \cdot P +1 = (8-1)\cdot 10 +1 = 71. \end{equation*}
(8)
Recalculating values $x_k,$ $n$ and $\Delta_k$ is also performed according to (1), (2) and (3) respectively.
The result of digital interpolation when using Farrow filter is shown in Figure 8b. Interpolation nodes corresponding to the input signal are black.
Figure 8. Input signal and result of digital interpolation when using Farrow filter.
The program source code implementing data calculation to plot Figure 8 is also given in the corresponding DSPL library reference topic.
GNU Octave and Matlab users can use a script resample_lagrange_interp.m for digital interpolation by using Farrow filter.

clear all;
close all;
clc;

addpath('functions');

p = 10;
q = 1;
x0 = 0;

s = [1 2 2 1 -0.5 -1 -2 -0.5];

y = resample_lagrange(s, p, q, x0);

figure(1);
stem(0:length(y)-1, y,'r');
hold on;
stem((0:length(s)-1)*p, s, 'k');


The filter of digital interpolation can be interpreted as the low-pass filter. The Farrow filter impulse response $h(n)$ and the magnitude $\left|H \left( \operatorname{e}^{j\cdot \omega}\right) \right|^2$ of digital signal interpolation are shown in Figure 9 for $P = 10$.

Figure 9. Farrow filter impulse response and magnitude.
The impulse response under $P \rightarrow\infty$ is a dotted graph.
In Figure 9 it is possible to note that the pulse response has no continuous derivative in interpolation nodes (has salient points) in view of the fact that plotting Lagrange polynomial interpolation is carried out only according to signal samples without restrictions to derivative continuity. Therefore the filter suppression in stopband is 28 dB only.
The program source code of calculating the Farrow filter impulse response and magnitude of digital signal interpolation (Figure 9) is also given in the DSPL library reference topic.
GNU Octave and Matlab users can calculate interpolator frequency response by script resample_lagrange_interp_filter.m.

clear all;
close all;
clc;

addpath('functions');

p = 10;
q = 1;
x0 = 0;

s = [0 0 1 0 0];

h = resample_lagrange(s, p, q, x0);

w = 0:0.01:pi;
H = freqz(h,1,w);

figure(1); subplot(1,2,1);
plot(0:length(h)-1, h, 'r'); hold on;
stem(0:length(h)-1, h);

xlabel('n'), ylabel('h(n)');
grid on;

figure(1); subplot(1,2,2);
plot(w/pi, 20*log10(abs(H)));
axis([0, 1, -60, 30]);
xlabel('w/pi');
ylabel('|H(e^{jw})|^2, dB'); grid on;

Using Farrow filter for fractional sample rate conversion
Let the sine input signal
\begin{equation*} s(n) = \sin \left( 2\pi \cdot n \cdot \frac{f_0}{F_{\textrm{s}}} \right), ~ n =0 \ldots N-1 \end{equation*}
(9)
with the frequency $f_0 = 6$ kHz contain $N = 54$ samples taken with the sampling frequency $F_{\textrm{s}}=26.4$ kHz.
The ratio of the sampling frequency $F_{\textrm{s}}$ of the sine signal $s(n)$ to this signal frequency $f_0$ is not an integer $\frac{F_{\textrm{s}}}{f_0} = \frac{26.4}{6}$. As a result the integer of samples is not located in one sine signal period, and the pickup signal is shown in the upper chart of Figure 10.
Change the sampling frequency of the $s(n)$ input signal $\frac{P}{Q}$ times, where $P = 20$, $Q = 11$ and we will receive the resampled signal $y(k)$ samples of which are received with the sampling frequency $F_{\textrm{y}} = \frac{P}{Q} \cdot F_{\textrm{s}} = 48$ kHz. As a result exactly 8 samples will be located in one signal period $y(k)$. Fractional resampling will be performed when using Farrow filter on the basis of polynomial interpolation. Recalculating the values $x_k,$ $n$ and $\Delta_k$ is performed according to (1), (2) and (3) respectively.
The signal $y(k)$ after digital resampling when using Farrow filter is shown in the lower chart of Figure 10.
Figure 10. Fractional change of the signal frequency sampling.
The pickup continuous sine signal is shown in Figure 10 with the help of continuous lines. Figure 10 shows that after resampling exactly eight signal samples $y(k)$ are located in one pickup signal period.
You can also find the program source code of calculating the data of the digital resampling example in the corresponding DSPL library reference topic.
There is GNU Octave and Matlab script resample_lagrange_pq.m for Farrow filter digital resampling.

clear all;
close all;
clc;

addpath('functions');

Fs = 26.4;
f0 = 6;
N = 54;

p = 20;
q = 11;
x0 = 0;

t = (0:N-1)/Fs;
s  = sin(2*pi*t*f0);
t0 = 0:0.01:2;
s0 = sin(2*pi*t0*f0);

y = resample_lagrange(s, p, q, x0);
tr = (0:length(y)-1)/Fs/p*q;

figure(1);
subplot(2,1,1); hold on;
stem(t, s);
%plot(t0 ,s0, 'r');
axis([0,2,-1.5,1.5]); grid on;

subplot(2,1,2);  hold on;
stem(tr, y);
%plot(t0 ,s0, 'r');
axis([0,2,-1.5,1.5]);
grid on;

Conclusions
In this section we considered a question of recalculating sample indexes of the input signal in case of piecewise and cubic interpolation. The equations for recalculating output signal sampling points according to the scale $n$ of the input signal, and also value calculation $\Delta_k$ according to the scale $t$ were received.
The received equations were used in examples of using Farrow filter to compensate fractional delay, digital interpolation and fractional sample rate conversion.
Families of impulse responses, magnitudes and group delay of Farrow filters for a various fractional delay value were given.
The Farrow filter impulse response and magnitude of digital signal interpolation are also given.
Links to program source codes realizing calculation of all shown examples using DSPL library are provided.

You can provide your any feedbacks, questions and wishes in the message board or guestbook.
Reference
[1] Kahaner D., Moler C., Nash S. Numerical Methods and Software. Prentice Hall, 1988.

[2] Farrow C.W. A Continuously Variable Digital Delay Element. Circuits and Systems, IEEE International Symposium. 1988, p. 2641–2645. vol. 3

[3] Gardner Floyd M. Interpolation in Digital Modems-Part I: Fundamentals: IEEE Transactions on Communications, Vol. 41, No. 3, March 1993, P. 501-507.

[4] Erup L., Gardner Floyd M., Harris Robert A. Interpolation in Digital Modems-Part II: Implementation and Performance: IEEE Transactions on Communications, Vol. 41, No. 6, June 1993, p.998-1008.

[5] Franck A. Efficient Algorithms for Arbitrary Sample Rate Conversion with Application to Wave Field Synthesis. PhD thesis. Universitätsverlag Ilmenau, Ilmenau, 2012. [PDF]

[6] McConnell J. Analysis of Algorithms: An Active Learning Approach. Jones and Bartlett Publishers, 2001.

 Select spelling error with your mouse and press