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 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 and two trivial multiplications by
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 of input signal samples , within the interval as it is shown in Figure 2.
Figure 2. Recalculating Indexes of Input Signal Samples within 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 everywhere in this section, as well as in the previous sections, the variable indexes samples of the input signal , and the variable indexes samples of the signal 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 th sample of the output signal after resampling does not correspond to the th sample of the input signal . Therefore we shall specify according to what scale we index a sample. 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 is within the interval to calculate coefficients of the cubic polynomial.
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 turn to the scale within the interval to calculate coefficients of the cubic polynomial (Figure 2b). At the same time the pickup time of the output sample according to the scale has to be recalculated into the value according to the scale .
The pickup time of the output sample after resampling according to the scale of the input signal is equal to:
(1)
(1)
To calculate the cubic polynomial 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 , and two other samples and have to be more to the left.
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 the next floor.
The value according to the scale is equal to:
(3)
(3)
In Figure 3 there is shown the example of recalculating indexes of the input and output signals for the th of the output sample under , , .
Figure 3. Example of Recalculating Indexes under , ,
The moment 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 polynomial according to the following equations:
(7)
(7)
then the value can be received as the interpolation spline value for , i.e. .
Using Farrow Filter on the Basis of Spline Interpolation to Compensate Fractional Signal Delay
Let the input signal contain 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 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 and are equal to , and the number of samples at the Farrow Filter output is equal to the number of samples at the Farrow Filter input .
The process of compensating the fractional signal delay shown in Figure 4 is given in Table 1.
Table 1. Compensating the Fractional Signal Delay .
Values and are calculated according to (1), (2) and (3) respectively. Further values for each value are given. The zero values for and being beyond the scope of indexing the signal are red. The values of coefficients of the cubic polynomial calculated according to (7) for the present moment are given in the next columns of Table 1. At last, in the column there is given the signal value at the Farrow Filter output with using spline interpolation. In the last column 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 (black) and the signal after compensating the fractional delay (red).
Figure 5. Input Signal of Farrow Filter and the Signal after Compensating the Fractional Delay .
It follows from Table 1 that when compensating the fractional delay parametric variables are for all . In this case Farrow Filter can be interpreted as the third order IIR-filter with changing AFR and group delay by setting the required value .
Comparative characteristics of AFR and group delay of Farrow Filter are given in Figure 6 when using spline interpolation (red) and Lagrange interpolation polynomial (blue) for the different parametric variable value .
Figure 6. AFR and Group Delay of Farrow Filters for a Different Value of the Fractional Delay 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 rad/s.
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 have 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 at the Farrow Filter output by 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):
(8)
(8)
Recalculating values and 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 digital interpolation filter can be interpreted as the filter of lower frequencies. The pulse response characteristic and amplitude-frequency response of Farrow Filter of digital signal interpolation are shown in Figure 8 for .

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).
Using Farrow Filter on the Basis of Spline Interpolation for Fractional Frequency Change of Signal Sampling
Let the sine input signal
(9)
(9)
with the frequency kHz have samples taken with the sampling frequency 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 of the sine signal to this signal frequency is not the integer . 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 by times where , and get resampled signals and received on the basis of polynomial and spline interpolation respectively. Samples of the signals and are received with the sampling frequency kHz. As a result exactly 8 samples will cover one signal period . Recalculating the values and is made according to (1), (2) and (3) respectively.
In the second chart of Figure 9 there is shown the signal after digital resampling when using Farrow Filter on the basis of Lagrange polynomial interpolation. In the third chart there is given the signal when using Farrow Filter on the basis of spline interpolation. In the last chart there are given signals of the resampling errors and 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 and .
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 is 0.0001 (0.01%) higher.
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 of the input signal, and also calculating the value according to the scale .
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.01% higher than the similar Farrow Filter on the basis of Lagrange polynomial interpolation has.
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];
    a1 = 0.5*(s[n] - s[n-2]);
    a3 = 2.0*(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] = 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




def resampling_spline(s, p, q, x0):
    """
    """
    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];
        a1 = 0.5*(s[n] - s[n-2]);
        a3 = 2.0*(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] = 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_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")


    h = resampling_spline(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_spline_filter_fd_mag_%.1f.csv' % x0[k]
    np.savetxt(fname, np.transpose([w, H]), fmt="%+.9e")
    
    fname = 'dat/resample_spline_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



def resampling_spline(s, p, q, x0):
    """
    """
    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];
        a1 = 0.5*(s[n] - s[n-2]);
        a3 = 2.0*(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] = 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_spline_interp_s.txt", np.transpose([t, s]), fmt="%+.9e")

y = resampling_spline(s, p, q, x0)

t = np.linspace(0, len(y), len(y), endpoint = False, dtype = 'float64')
np.savetxt("dat/resample_spline_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")


h = resampling_spline(s, p, q, x0)

t = np.linspace(0, len(h), len(h), endpoint = False, dtype = 'float64')
np.savetxt("dat/resample_spline_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_spline_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
    
    

def resampling_spline(s, p, q, x0):
    """
    """
    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];
        a1 = 0.5*(s[n] - s[n-2]);
        a3 = 2.0*(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] = 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)

yL = resampling_lagrange(s, P, Q, 0.0)
yS = resampling_spline(s, P, Q, 0.0)
ty = float(Q) * np.linspace(0, len(yL), len(yL), endpoint = False, dtype = 'float64') / (Fs * float(P))

eL = yL - np.sin(np.pi*2.0*f0*ty)
eS = yS - np.sin(np.pi*2.0*f0*ty)

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[0:len(ty)-2], yL[0:len(ty)-2]]), fmt="%+.9e")
np.savetxt("dat/resample_lagrange_ex_fs_e.txt", np.transpose([ty[0:len(ty)-2],eL[0:len(ty)-2]]), fmt="%+.9e")
np.savetxt("dat/resample_spline_ex_fs_y.txt", np.transpose([ty[0:len(ty)-2], yS[0:len(ty)-2]]), fmt="%+.9e")
np.savetxt("dat/resample_spline_ex_fs_e.txt", np.transpose([ty[0:len(ty)-2], eS[0:len(ty)-2]]), fmt="%+.9e")


nmseL = np.sum(eL**2)/np.sum(np.sin(np.pi*2.0*f0*ty)**2)
nmseS = np.sum(eS**2)/np.sum(np.sin(np.pi*2.0*f0*ty)**2)

#print("nmseL = " ,nmseL) 
#print("nmseS = ", nmseS)