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 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")