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 and two trivial multiplication by .

In this section we will consider examples of using the resampling filter to compensate
fractional delay, digital interpolation and fractional resampling 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 ,
in the interval , as it is shown in Figure 2.

Figure 2. Recalculating sample indexes of an input signal in the interval

Let the sampling frequency of the input signal ,
be equal to .
Then the signal , at the resampler filter output
will have the following sampling frequency .
Besides the first sample of the output signal is delayed
concerning the input one by the value .

Here and elsewhere in this section as well as in the previous section
the variable indexes samples of the input signal ,
and the variable indexes samples of the output signal .

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 th sample of the output signal after resampling
does not correspond to the th sample of the input signal .
Therefore we have to specify what scale we use to index samples.
We will use three scales: the scale is for timepoints
concerning samples of the input signal ,
the scale is for timepoints concerning the output signal and the scale
in the interval is for calculating coefficients of a cubic polynom.

For the sample of the output signal according to the scale
it is necessary to calculate the corresponding indexes of the input signal
according to the scale (Figure 2a), and then to pass to the scale
in the interval for calculating coefficients of a cubic polynom (Figure 2b).
At that the output sample point according to the scale has to be recalculated under
according to the scale .

The point of the output sample after resampling according
to the scale of the input signal is equal to:

(1)

(1)

To calculate a cubic polynom we need four samples of the input signal ,
, and , as it is shown in Figure 2a.
At that two samples and have to be more to the right of
and two other samples and have to be more to the left of it.

Then the index of the sample which corresponds to the timepoint
according to the scale from -2 to 1 is equal to:

(2)

(2)

where the functional means rounding to floor.

The value according to the scale is equal to:

(3)

>

In Figure 3 the example of recalculating indexes of the input and output
signals for the th output sample under , , is shown.

Figure 3. Example of recalculating indexes under , ,

The point according to the scale in accordance with (1) is equal to:

(4)

(4)

Then the index in accordance with (2) is equal to:

(5)

(5)

and the parametric variable according to the scale is equal to:

(6)

(6)

Thus, it is necessary to calculate coefficients of the cubic polynom
according
to the following equations:

(7)

(7)

then the value can be got as the interpolation polynom value
for , i.e. .

Using Farrow filter on the basis of polynomial interpolation to compensate the signal fractional delay

Let the input signal contain 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 of the sampling period.
Upon changing the fractional delay the sampling frequency is unchanged, the parametric
variables and are equal to and the number of samples at the Farrow
filter output is equal to the number of samples at the input, i.e. .

The calculation process of the signal fractional delay
is given in Table 1. This process is shown in Figure 4.

Table 1. Equalization of the signal fractional delay .

The values and are calculated according to (1), (2) and (3) respectively.
Then the values for each value are given.
Zero values for and which fall outside the
limits of the signal indexation are red.
In the following columns of Table 1 the coefficient values of the cubic polynom
calculated according to (7) for the current point are given.
At last, the signal value at the Farrow filter output is
given in the last column.

The pickup signal (black) and the signal
after equalizing the fractional delay (red) are shown in Figure 5.

Figure 5. Input signal of Farrow filter and the signal after equalizing the fractional delay .

Table 1 shows that upon equalizing the fractional delay the parametric variables
for all .
In this case Farrow filter can be interpreted as the third order FIR filter
with the variable group delay by specifying the required value .

The family of Farrow filter impulse response for a various fractional delay value
is shown in Figure 6.

Figure 6. Farrow filter impulse responses for a various fractional delay value

The magnitude
and the group delay of received filters are shown in Figure 7.

Figure 7. Magnitude and group delay of Farrow filters for a various fractional delay value

The magnitude and the group delay charts allow to draw a conclusion that delay
equalization is possible only within the scope of up to 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.

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 contain 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 at the Farrow filter output times.

In case of digital interpolation the number of signal samples at the output will be equal to (see Figure 8b):

(8)

(8)

Recalculating values and 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 filter of digital interpolation can be interpreted as the low-pass filter.
The Farrow filter impulse response and the magnitude of digital signal interpolation are shown in Figure 9 for .

Figure 9. Farrow filter impulse response and magnitude.

The impulse response under 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.

Using Farrow filter for fractional sample rate conversion

Let the sine input signal

(9)

(9)

with the frequency kHz contain samples taken with the
sampling frequency kHz.

The ratio of the sampling frequency of the sine signal
to this signal frequency is not an integer
.
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 input signal times,
where , and we will receive the resampled signal samples
of which are received with the sampling frequency
kHz.
As a result exactly 8 samples will be located in one signal period .
Fractional resampling will be performed when using Farrow filter on
the basis of polynomial interpolation.
Recalculating the values and is performed according
to (1), (2) and (3) respectively.

The signal 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 are
located in one pickup signal period.

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 of the input signal, and also value calculation
according to the scale 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.

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.

Appendix

Python simulation scripts:

# -*- coding: utf-8 -*-
"""
Created on Sun Aug 27 23:40:00 2017
@author: sergey
"""
import numpy as np
p = 1
q = 1
x0 = 0.25
s0 = np.array([1., 2., 2., 1., -0.5, -1., -2., -0.5])
tin = np.linspace(0.0, len(s0), len(s0), endpoint=False)
y = np.zeros(int(np.floor( float(len(s0)*q)/float(p))))
s = np.concatenate((np.array([0., 0.]), s0, np.array([0., 0.])))
print('---------------------------------------------------------', end='');
print('-------------------------------------------------------\n', end='');
print('k x_k n d s(n-3) s(n-2) s(n-1) s(n) ', end='');
print('a_0 a_1 a_2 a_3 y(k)\n', end='');
print('---------------------------------------------------------', end='');
print('-------------------------------------------------------\n', end='');
tout = np.zeros(int(np.floor( float(len(s0)*q)/float(p))))
for k in range(len(y)):
x = k*q/p - x0
n = int(x) + 1 + 3
d = int(x) + 1 - x
tout[k] = 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] = a0 - a1 * d + a2*d**2 - a3*d**3
print('%(d0)d %(f0)7.2f %(d1)2d %(f1)7.2f' %
{'d0': k, 'f0': x, 'd1': n, 'f1': d}, end='')
print('%(s0)8.1f %(s1)8.1f %(s2)8.1f %(s3)8.1f ' %
{'s0': s[n-3],'s1': s[n-2], 's2': s[n-1], 's3':s[n]}, end='')
print('%(a0)9.3f %(a1)9.3f %(a2)9.3f %(a3)9.3f '%
{'a0': a0, 'a1': a1, 'a2': a2, 'a3': a3}, end='')
print('%(y)9.4f\n' % {'y': y[k]}, end='')
print('---------------------------------------------------------', end='');
print('-------------------------------------------------------\n', end='');
"""
figure; stem(tin, s0);
hold on;
stem(tout, y, 'r');
axis([-0.5, 7.5, -2.5, 2.5]);
grid on;
"""

# -*- coding: utf-8 -*-
"""
Created on Mon Aug 28 23:34:59 2017
@author: sergey
"""
import numpy as np
import scipy.signal as signal
def resampling_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):
if(q==1):
y = np.zeros(int(float((len(s)-1) * p) / float(q)) + 1)
else:
y = np.zeros(int( float(len(s) * p ) / float(q)))
else:
y = np.zeros(int( float(len(s) * p ) / float(q)))
t = np.zeros(len(y))
s = np.concatenate((np.array([0., 0.]), s, np.array([0., 0.])))
for k in range(len(y)):
x = k*q/p - x0
t[k] = x
n = int(np.floor(x)) + 4
d = np.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] = a0 - a1 * d + a2*d**2 - a3*d**3
return y
x0 = np.linspace(0.0, 0.9, 10)
s = np.array([0., 1., 0., 0.])
t = np.array([0., 1., 2., 3.])
for k in range(len(x0)):
h = resampling_lagrange(s, 1, 1, x0[k])
w, H = signal.freqz(h)
H = 20.0 * np.log10(np.abs(H))
w, gd = signal.group_delay((h, 1))
fname = 'dat/resample_lagrange_filter_fd_time_%.1f.csv' % x0[k]
np.savetxt(fname, np.transpose([t, h]), fmt="%+.9e")
fname = 'dat/resample_lagrange_filter_fd_mag_%.1f.csv' % x0[k]
np.savetxt(fname, np.transpose([w, H]), fmt="%+.9e")
fname = 'dat/resample_lagrange_filter_fd_gd_%.1f.csv' % x0[k]
np.savetxt(fname, np.transpose([w, gd]), fmt="%+.9e")

# -*- coding: utf-8 -*-
"""
Created on Mon Aug 28 23:34:59 2017
@author: sergey
"""
import numpy as np
import scipy.signal as signal
def resampling_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):
if(q==1):
y = np.zeros(int(float((len(s)-1) * p) / float(q)) + 1)
else:
y = np.zeros(int( float(len(s) * p ) / float(q)))
else:
y = np.zeros(int( float(len(s) * p ) / float(q)))
t = np.zeros(len(y))
s = np.concatenate((np.array([0., 0.]), s, np.array([0., 0.])))
for k in range(len(y)):
x = k*q/p - x0
t[k] = x
n = int(np.floor(x)) + 4
d = np.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] = a0 - a1 * d + a2*d**2 - a3*d**3
return y
p = 10
q = 1
x0 = 0
"""
Пример интерполяции сигнала с использованием фильтра Фарроу
"""
s = np.array([1., 2., 2., 1., -0.5, -1., -2., -0.5])
t = np.linspace(0, len(s)*p, len(s), endpoint = False, dtype = 'float64')
np.savetxt("dat/resample_lagrange_interp_s.txt", np.transpose([t, s]), fmt="%+.9e")
y = resampling_lagrange(s, p, q, x0)
t = np.linspace(0, len(y), len(y), endpoint = False, dtype = 'float64')
np.savetxt("dat/resample_lagrange_interp_y.txt", np.transpose([t, y]), fmt="%+.9e")
"""
Расчет импульсной характеристики h[k] и частотной характеристики H(w)
"""
s = np.array([0., 0., 1., 0., 0.])
h = resampling_lagrange(s, p, q, x0)
t = np.linspace(0, len(h), len(h), endpoint = False, dtype = 'float64')
np.savetxt("dat/resample_lagrange_interp_filter_time.txt", np.transpose([t, h]), fmt="%+.9e")
w, H = signal.freqz(h)
H = 20.0 * np.log10(np.abs(H))
np.savetxt("dat/resample_lagrange_interp_filter_freq.txt", np.transpose([w, H]), fmt="%+.9e")

# -*- coding: utf-8 -*-
"""
Created on Sat Sep 2 12:18:52 2017
@author: sergey
"""
# -*- coding: utf-8 -*-
"""
Created on Mon Aug 28 23:34:59 2017
@author: sergey
"""
import numpy as np
import scipy.signal as signal
def resampling_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):
if(q==1):
y = np.zeros(int(float((len(s)-1) * p) / float(q)) + 1)
else:
y = np.zeros(int( float(len(s) * p ) / float(q)))
else:
y = np.zeros(int( float(len(s) * p ) / float(q)))
t = np.zeros(len(y))
s = np.concatenate((np.array([0., 0.]), s, np.array([0., 0.])))
for k in range(len(y)):
x = k*q/p - x0
t[k] = x
n = int(np.floor(x)) + 4
d = np.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] = a0 - a1 * d + a2*d**2 - a3*d**3
return y
P = 20
Q = 11
N = 54
L = 10
Fs = 26.4
f0 = 6.0
tc = np.linspace(0, N*L, N*L, endpoint = False, dtype = 'float64')/(Fs * float(L))
c = np.sin(np.pi*2.0*f0*tc)
ts = np.linspace(0, N, N, endpoint = False, dtype = 'float64')/Fs
s = np.sin(np.pi*2.0*f0*ts)
y = resampling_lagrange(s, P, Q, 0.0)
ty = float(Q) * np.linspace(0, len(y), len(y), endpoint = False, dtype = 'float64') / (Fs * float(P))
np.savetxt("dat/resample_lagrange_ex_fs_c.txt", np.transpose([tc[0:(N-1)*L], c[0:(N-1)*L] ]), fmt="%+.9e")
np.savetxt("dat/resample_lagrange_ex_fs_s.txt", np.transpose([ts, s]), fmt="%+.9e")
np.savetxt("dat/resample_lagrange_ex_fs_y.txt", np.transpose([ty, y]), fmt="%+.9e")