Using Farrow Filter of Signal Digital Resampling on the Basis of Spline Interpolation

Contents
Introduction
In the previous section there was considered Farrow Filter of signal digital resampling on the basis of cubic Hermite splines (see Figure 1).
Figure 1. Functional Chart of the Optimized Farrow Filter
on the Basis of Cubic Hermite Splines.
The received filter requires only one multiplication by $\frac{1}{2}$ to calculate spline coefficients which can be considered trivial.
Earlier the functional chart of Farrow Filter of digital resampling on the basis of Lagrange polynomial interpolation and examples of its use were considered. We said that Farrow Filter on the basis of Lagrange polynomial interpolation requires one multiplication by $\frac{1}{6}$ and two trivial multiplications by $\frac{1}{2}.$
In this section we will review examples of using Farrow Filter on the basis of cubic Hermite splines and we will compare its characteristics with similar characteristics of the filter on the basis of Lagrange polynomial interpolation.
Recalculating Indexes of Input Signal Samples in case of Piecewise and Cubic Spline Interpolation
We considered this issue in the previous section when we were reviewing examples of using Farrow Filter on the basis of Lagrange polynomial interpolation. When using spline interpolation equations to recalculate indexes of input signal samples are constant, but to give a complete picture we will provide them once again.
Great efficiency of the digital resampler is reached due to using piecewise and cubic spline interpolation and recalculating indexes $n$ of input signal samples $s(n)$, $n =0,1,2, \ldots$ within the interval $[-2, \ldots 1]$ as it is shown in Figure 2.
Figure 2. Recalculating Indexes of Input Signal Samples within 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 everywhere in this section, as well as in the previous sections, the variable $n$ indexes samples of the input signal $s(n)$, and the variable $k$ indexes samples of the signal $y(k)$ at the resampler filter output.
As sampling frequencies of input and output signals are various, the same indexes of 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 shall specify according to what scale we index a sample. 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$ is within the interval $[-2 \ldots 1]$ to calculate coefficients of the cubic polynomial.
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 turn to the scale $t$ within the interval $[-2 \ldots 1]$ to calculate coefficients of the cubic polynomial (Figure 2b). At the same time the pickup time of the $k$ output sample according to the scale $n$ has to be recalculated into the value $0 \leq \Delta_k < 1$ according to the scale $t$.
The pickup time of the $k$ output sample $x_k$ 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 the cubic polynomial 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 $x_k$, and two other samples $s(n-2)$ and $s(n-3)$ have to be more to the left.
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 the next 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 there is shown the example of recalculating indexes of the input and output signals for the $5-$th of the output sample $y(5)$ under $P = 4$, $Q = 3$, $x_0 = 0.2$.
Figure 3. Example of Recalculating Indexes under $P = 4$, $Q = 3$, $x_0 = 0.2$
The moment $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 polynomial $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_1 = \frac{1}{2}\cdot \Big(s(n) - s(n-2)\Big), \\ a_3 = 2 \cdot \Big(s(n-2) - s(n-1)\Big) + a_1 + \frac{1}{2}\cdot \Big(s(n-1) - s(n-3)\Big), \\ a_2 = s(n-2)-s(n-1) + a_3 + a_1. \end{cases} \end{equation*}
(7)
then the value $y(k)$ can be received as the interpolation spline value $P(t)$ for $t=-\Delta_k$, i.e. $y(k) = P \left(-\Delta_k \right)$.
Using Farrow Filter on the Basis of Spline Interpolation to Compensate Fractional Signal Delay
Let the input signal $s(n)$ contain $N=8$ samples as it is shown in Figure 4. We have already used a similar signal to compensate the fractional signal delay when using Farrow Filter on the basis of Lagrange polynomial interpolation

Figure 4. Input Signal of Farrow Filter.
Compensate the fractional delay of this signal by the value $x_0 = 0.25$ of the sampling period with using the resampler on the basis of spline interpolation. In case of changing the fractional delay the sampling frequency is constant, 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 Farrow Filter input $N=8$.
The process of compensating the fractional signal delay $s(n),$ $n = 0 \ldots 7,$ shown in Figure 4 is given in Table 1.
Table 1. Compensating the Fractional Signal Delay $s(n)$.
Values $x_k,$ $n$ and $\Delta_k$ are calculated according to (1), (2) and (3) respectively. Further values $s(n-3) \ldots s(n)$ for each value $k$ are given. The zero values $s(n-3) \ldots s(n)$ for $n<3$ and $n>N-1$ being beyond the scope of indexing the signal $s(n),$ $n = 0 \ldots N-1$ are red. The values of coefficients of the cubic polynomial $a_0 \ldots a_3$ calculated according to (7) for the present moment $k$ are given in the next columns of Table 1. At last, in the column $y_S(k)$ there is given the signal value $y_S(k) = P(-\Delta_k)$ at the Farrow Filter output with using spline interpolation. In the last column $y_L(k)$ values are given at the Farrow Filter output with using Lagrange polynomial interpolation as ones to compare. As it is possible to see with the help of Table 1 we have received similar values at the Farrow Filter output when using Lagrange spline interpolation and Lagrange interpolation.
In Figure 5 there are shown the pickup signal $s(n),$ $n = 0 \ldots 7$ (black) and the signal $y_S(k),$ $n = 0 \ldots 7$ after compensating the fractional delay $x_0$ (red).
Figure 5. Input Signal $s(n)$ of Farrow Filter and the Signal after Compensating the Fractional Delay $y_S(k)$.
The program realizing calculation of Table 1 and data for creating Figure 5 is written using DSPL library. The source program code can be received in the DSPL library examples section.
Users of Matlab and GNU Octave can make calculation of Table 1 having executed the following script resample_spline_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_S(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);
a1 = 0.5 * (s(n) - s(n-2));
a3 = 2*(s(n-2) - s(n-1))+ a1 + 0.5 * (s(n-1) - s(n-3));
a2 =  s(n-2) - s(n-1) + a3 + a1;

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

fprintf('%d  %7.2f   %2.0d   %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;


It follows from Table 1 that when compensating the fractional delay parametric variables $\Delta_k = x_0$ are for all $k = 0 \ldots N-1$. In this case Farrow Filter can be interpreted as the third order IIR-filter with changing AFR $\left| H \left(\textrm{e}^{j\cdot \omega} \right) \right|^2$ and group delay $\tau(\omega)$ by setting the required value $x_0$.
Comparative characteristics of AFR $\left| H \left(\textrm{e}^{j\cdot \omega} \right) \right|^2$ and group delay $\tau(\omega)$ of Farrow Filter are given in Figure 6 when using spline interpolation (red) and Lagrange interpolation polynomial (blue) for the different parametric variable value $x_0$.
Figure 6. AFR and Group Delay of Farrow Filters for a Different Value of the Fractional Delay $x_0$ when Using Spline Interpolation (red) and Lagrange Interpolation Polynomial (blue)
Charts of AFR and group delay help draw a conclusion that Farrow Filter when using spline interpolation has smaller non-uniformity of AFR, but slightly bigger non-uniformity of group delay in comparison with Farrow Filter on the basis of Lagrange polynomial interpolation. At the same time changing the time delay is possible only up to $0.4\pi$ rad/s.
Calculating AFR and group delay of Farrow Filters with using spline interpolation for a different value of the fractional delay $x_0$ shown in Figure 6 was also made when using DSPL library. The source program code of calculating AFR and group delay of Farrow Filters can also be received in the corresponding DSPL library examples section.
Users of Matlab and GNU Octave can make calculation of AFR and group delay of Farrow Filters for a different value of the fractional delay $x_0$ having executed the following script resample_spline_fd_filter.m.

clear all;
close all;
clc;

% Check Matlab or GNU Octave
isOctave = exist('OCTAVE_VERSION', 'builtin');
if(isOctave)
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)
hl = resample_spline(s, 1, 1, x0(k));
hs = resample_lagrange(s, 1, 1, x0(k));

figure(1); subplot(2,5,k); hold on; stem(0:3, hl), stem(0:3, hs, 'sr');
axis([0, 3, -0.2, 1.2]); grid on;

HL = freqz(hl, 1, w);
HS = freqz(hs, 1, w);
tauL = grpdelay(hl,1,w);
tauS = grpdelay(hs,1,w);

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

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

end

To execute the script resample_spline_fd_filter.m and other scripts given below the following function resample_spline.m is required.

function y = resample_spline(s, p, q, x0)
% y = resample_spline(s, p, q, x0)
% Digital resampling by polynomial spline 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);
a1 = 0.5 * (s(n) - s(n-2));
a3 = 2*(s(n-2) - s(n-1))+ a1 + 0.5 * (s(n-1) - s(n-3));
a2 =  s(n-2) - s(n-1) + a3 + a1;

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

end

Using Farrow Filter on the Basis of Cubic Hermite Splines as a Digital Signal Interpolator
Consider an example of using the digital resampling filter on the basis of cubic Hermite splines for digital signal interpolation.
Let the input signal $s(n),$ $n = 0 \ldots N-1,$ have $N=8$ signal samples as it is shown in Figure 7a. We have already used this signal as an input signal of the digital filter for changing the fractional delay.
It is necessary to perform digital signal interpolation and to increase the sampling frequency of the signal $y(k),$ $k = 0\ldots K-1$ at the Farrow Filter output by $P=10$ times when using spline interpolation.
In case of digital interpolation the number of signal samples at the output will be equal to (see Figure 7b):
\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 made according to (1), (2) and (3) respectively.
The result of digital interpolation when using Farrow Filter is shown in Figure 7b. Interpolation nodes corresponding to the input signal are black.
Figure 7. Input Signal and the Result of Digital Interpolation when Using Farrow Filter.
The source program code realizing calculating the data for creating Figure 7 is also provided in the corresponding DSPL library reference section.
Users of Matlab and GNU Octave can make digital interpolation with using Farrow Filter having executed the following script resample_spline_interp.m.

clear all;
close all;
clc;

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

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

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

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

The digital interpolation filter can be interpreted as the filter of lower frequencies. The pulse response characteristic $h(n)$ and amplitude-frequency response $\left|H \left( \operatorname{e}^{j\cdot \omega}\right) \right|^2$ of Farrow Filter of digital signal interpolation are shown in Figure 8 for $P = 10$.

Figure 8. Pulse Response Characteristic and AFR of Farrow Filter on the Basis of Spline Interpolation (Red) and Lagrange Polynomials (Blue)
In Figure 8 it is possible to note that the pulse response characteristic has no continuous derivative in interpolation nodes (it has breakpoints) when using Lagrange polynomials in view of the fact that creating Lagrange interpolation polynomial is performed only according to signal samples without restrictions on derivative continuity. Therefore the rejection level in the scope of the filter AFR blocking is only 28 dB. On the other hand using Hermite splines provides the continuous derivative of the pulse response characteristic of Farrow Filter in interpolation nodes. As a result the rejection level in the scope of the Farrow filter AFR blocking on interpolation spline (the red chart) is 12 dB higher than the same filter, but it is on the basis of Lagrange polynomial interpolation (the blue chart).
The source program code realizing calculating the pulse response characteristic and AFR of Farrow Filter of digital signal interpolation (Figure 9) is also provided in the DSPL library reference section.
Users of Matlab and GNU Octave can make calculation of the digital interpolator filter characteristics having executed the following script resample_spline_interp_filter.m.


clear all;
close all;
clc;

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

s = [0 0 1 0 0];

hl = resample_lagrange(s, p, q, x0);
hs = resample_spline(s, p, q, x0);

w = 0:0.01:pi;
HL = freqz(hl,1,w);
HS = freqz(hs,1,w);

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

plot(0:length(hs)-1, hs, 'r'); hold on;
stem(0:length(hs)-1, hs, 'rs');

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

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

Using Farrow Filter on the Basis of Spline Interpolation for Fractional Frequency Change of Signal Sampling
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 have $N = 54$ samples taken with the sampling frequency $F_{\textrm{s}}=26.4$ kHz. We have already used this signal for an example of fractional resampling using Farrow Filter on the basis of Lagrange polynomial interpolation.
The sampling rate ratio $F_{\textrm{s}}$ of the sine signal $s(n)$ to this signal frequency $f_0$ is not the integer $\frac{F_{\textrm{s}}}{f_0} = \frac{26.4}{6}$. As a result the integral number of samples does not cover one sine signal period. The pickup signal is shown in the upper chart of Figure 9.
Change the sampling frequency of the input signal $s(n)$ by $\frac{P}{Q}$ times where $P = 20$, $Q = 11$ and get resampled signals $y_L(k)$ and $y_S(k)$ received on the basis of polynomial and spline interpolation respectively. Samples of the signals $y_L(k)$ and $y_S(k)$ 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 cover one signal period $y(k)$. Recalculating the values $x_k,$ $n$ and $\Delta_k$ is made according to (1), (2) and (3) respectively.
In the second chart of Figure 9 there is shown the signal $y_L(k)$ after digital resampling when using Farrow Filter on the basis of Lagrange polynomial interpolation. In the third chart there is given the signal $y_S(k)$ when using Farrow Filter on the basis of spline interpolation. In the last chart there are given signals of the resampling errors $e_L(k)$ and $e_S(k)$ at the output of Farrow Filters when using Lagrange polynomial interpolation and spline interpolation respectively.
In the lower chart there are given accessed values of normalized average square of the resampling errors $\textrm{NMSE} \big[ e_L(k)\big]$ and $\textrm{NMSE} \big[ e_S(k)\big]$.
Figure 9. Fractional Frequency Change of Signal Sampling.
Continuous lines in Figure 9 show the pickup continuous sine signal. It follows from Figure 9 that after resampling the eight signal samples cover one pickup signal period either when using Lagrange interpolation or when using spline interpolation. At the same time when using spline interpolation the normalized root-mean-square error $\textrm{NMSE} \big[ e_S(k)\big]$ is 0.001 (0.1%) higher.
You can also find the source program code of calculating the data of a digital resampling example in the corresponding DSPL library reference section.
Users of Matlab and GNU Octave can make digital sampling having executed the foolowing script resample_spline_pq.m.

clear all;
close all;
clc;

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

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

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

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

sr = sin(2*pi*tr*f0);

ys = resample_spline(s, p, q, x0);

figure(1);
subplot(3,1,1); hold on;
stem(t, s);
axis([0,2,-1.5,1.5]); grid on;

subplot(3,1,2);  hold on;
stem(tr, yl); hold on; stem( tr, ys,'rs');
axis([0,2,-1.5,1.5]);
grid on;

subplot(3,1,3);  hold on;
stem(tr, yl-sr); hold on; stem( tr, ys-sr,'rs');
axis([0,2,-0.1,0.1]);
grid on;

Conclusions
In this section we have reviewed examples of using Farrow Filter on the basis of spline interpolation and comparison of its characteristics with Farrow Filter on the basis of Lagrange polynomial interpolation. The main benefit of Farrow Filter on the basis of spline interpolation is in use only of trivial multiplications to calculate coefficients of an interpolation polynomial.
In this section we have considered a question of recalculating indexes of input signal samples in case of piecewise and cubic spline interpolation. There are given equations for recalculating sampling timepoints of the output signal according to the scale $n$ of the input signal, and also calculating the value $\Delta_k$ according to the scale $t$.
The given equations are used in examples of using Farrow Filter on the basis of spline interpolation to compensate the fractional delay, digital interpolation and fractional change of signal sampling frequency.
Comparisons of AFR and group delay of Farrow Filters families when using Lagrange polynomial interpolation and spline interpolation for a various fractional delay value are given.
Comparisons of Farrow Filter AFR of digital signal interpolation when using Lagrange polynomial interpolation and spline interpolation are also given. It is shown that the resampling filter on the basis of spline interpolation provides a continuous derivative of the interpolated signal at this expense rejection in the scope of the filter blocking is 12 dB better.
Comparison of fractional sine signal resampling when using Farrow Filters on the basis of Lagrange polynomial interpolation and spline interpolation is carried out. It is received that the rated RMSE at the Farrow Filter output on the basis of spline interpolation is 0.1% higher than the similar Farrow Filter on the basis of Lagrange polynomial interpolation has.
Links to the initial program codes realizing calculating all the shown examples with using DSPL library and packages of Matlab and GNU Octave are given.
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