**Interpolation in Digital
Signal Processing and Numerical Analysis **

*Saeed Babaeizadeh*

Spring 2003

In many practical applications of digital signal processing, one is faced with the problem of changing the sampling rate of a signal, either increasing it or decreasing it by some amount. An increase in the sampling rate by an integer factor L can be accomplished by interpolating L-1 new samples between successive values of the signal. The interpolation process can be accomplished in a variety of ways in both fields of digital signal processing or numerical analysis.

There is one very simple and straightforward approach to changing the sampling rate of a digital signal. This approach, called the analog approach, merely reconstructs the continuous-time signal from the original set of samples and then resamples the signal at the new rate. In practice this approach suffers from one major problem, namely, that the ideal operations required reconstructing the continuous-time signal from the original samples and to resample the signal at the new rate cannot be implemented exactly [1]. An alternative is the so-called direct digital approach that does the job without going through D/A or A/D conversion.

In DSP we use a process that preserves the spectral
shape of the signal sequence *x*(*n*) [2]. Here I describe a process
to increase the sample rate by an integer factor [18].

If the sampling rate is increased by the integer
factor *L*, then the new sampling rate will be samp./sec., thus the new sampling interval
will be .
Let the original signal be *x*(*n*) and the interpolated high rate
signal be *y*(*n*). If we define an intermediate signal with the
sample rate of

(1)

, it is known [1,2,3,18] that the Fourier transforms of them are related as:

(2)

Equation (2) says that the Fourier transform *V*
is just the same as that of *X* except that the frequency scale is
compressed by the factor *L*.

Now, the Fourier transform *Y* that we require,
should be the same as one period of *X* but compressed to an *L*’th
band of the full band .
There is also a scale factor of *L* which is necessary to make the inverse
Fourier transform come out right. Hence we need:

(3)

To obtain *Y* from *V*, we need to pull
out the first period of *V* and scale it by the factor *L*. This, can
be done by using the ideal low pass filter (LPF)

(4)

Then *y*(*n*)=*v*(*n*)**h*(*n*)
at rate .
Thus, the block diagram of the ideal interpolator is:

**
Figure 1)
**The block diagram of the ideal
interpolator

In practical designs we must replace the ideal LPF by a more practical finite order filter. Any low pass filter can be used, but typically, the filter is a FIR filter, because using a FIR one can do “random access” on the data with no extra computational cost. It is a useful property if N is high, because in such cases typically only a fraction of the samples in the oversampled signal are used.

When we replace *H* by an FIR filter, we must
allow for a transition band, thus, the filter cannot be flat up to

Define

The highest
frequency of the signal band (of signal *x*(*n*)) to be preserved.

Full signal
bandwidth of *x*(*n*), i.e. there is no energy in *x*(*n*)
above this frequency.

Now, the filter specifications can be designed to pass the band and to have stopband starting at . Hence, we have the following filter design problem:

(5)

There are so many different methods to design this FIR filter, but those methods are behind the scope of this paper.

We know that when we are using a finite number
sequence to design a low pass FIR filter, the result filter would have a
transition band for going from passband to stopband. Also the filter oscillates
in passband and stopband and it has some sidelobes in frequency domain which
result in leakage of energy to other frequencies. All these, affect the
accuracy of the interpolation and make some distortions in the result. Besides,
as it is seen in Fig(1), we have to filter *v*(*n*) and it means we
ignore the energy of signal components in frequencies more than the cut-off
frequency of the filter, and it is another source of error.

Hence, this method is good when the signal is band limited in frequency domain and we have to try to design a LPF with less sidelobes and smaller transition band. An slightly different interpretation of this method, which is also called the “FFT method”, has given by Fraser [17].

In the field of numerical analysis, the
interpolation is to create an approximation of the continuous signal, from the
information contained in the samples, and to sample that. A common
interpolation method is linear interpolation, where the continuous function is
approximated as piecewise-linear by drawing lines between the successive
samples. An even more crude form of interpolation is drop-sample interpolation,
drawing a horizontal line from each sample until the following sample. Both of
these methods fall into the category of piecewise polynomial interpolators. For
modeling or interpolating a function, we can use polynomial models. However, if
the order of the polynomial is too small, it will not result in a good
approximation for the rapid varying parts of the function. Conversely, if the
order is chosen too large, the estimated function may be too varying in the
other points. In other words, as shown in [4], with a polynomial approximation,
“*if the function to be approximated is badly behaved anywhere in the
interval of approximation, then the approximate is poor everywhere”*.
However, by using low order piecewise polynomials (For instance 2^{nd}
or 3^{rd} degree splines) we can well approximate rapid varying parts
of a function without affecting the other parts of the function.

A polynomial interpolator can and should be thought as a filter with a continuous-time impulse response [5]. A non-discrete impulse response yields a non-periodic frequency response that has an overall descending envelope. So, the spectral sidelobes attenuated by this continuous filter, making resampling a more sensible process. Ideally there would be no sidelobes, as the continuous signal that we are dealing with, is presumed to be bandlimited in range when is the sampling frequency. The goal is to have the sidelobes attenuated as much as possible. Polynomial interpolators don’t have a flat passband.

Piecewise polynomial interpolation in this context means that individual polynomials are created between successive sample points. These interpolators can be classified both by the number of sample points required from the neighborhood for calculating the value of the position, and by the order of the polynomials. Depending on the interpolator, the polynomial order is typically one less than the number of points, matching the number of coefficients in the polynomial to the number of samples, but there are many exceptions to this rule [5].

Here, I only consider the interpolators that operate on an odd number of points, which have impulse responses symmetrical around zero.

The only piecewise polynomials that are symmetrical and, thus, do not introduce phase distortions are those corresponding to the case where the order of the polynomials is odd and the number of constituent samples on either side of the interpolation interval is the same [9]. So, in the following section, I just will show some examples of odd-order piecewise polynomials.

B-splines are a family of interpolators that can be constructed by convolving a drop-sample interpolator by a drop-sample interpolator repeated times [5].

When contrasted with local or running polynomial
fits, the use of B-splines functions, which are piecewise polynomials as well,
seems to have a number of advantages [6]. First, higher order polynomials tend
to oscillate while spline functions are usually smooth and well-behaved.
Second, the juxtaposition of local polynomial approximations may produce strong
discontinuities in the connecting regions. B-spline surfaces, by contrast, are
continuous everywhere. The polynomial segments are patched together so that the
interpolating function and all derivatives up to order *n-1* are
continuous at all joining intersections. Third, there is exactly one B-spline
coefficient associated with each grid point and the range of these coefficients
is of the same order of magnitude as that of the initial gray level values [7].
Finally, either exact or approximate B-spline signal representations can be
evaluated quite efficiently [6].

Drop-sample is the 0^{th}-order B-spline
that operates only on one point, its impulse response is:

(6) *Drop-sample interpolator impulse response*

The first three odd-order B-spline impulse responses are

(7) *2-point, 1 ^{st}-order linear interpolator
impulse response.* (8)

(9) *6-point,
5 ^{th}-order B-spline impulse response*

These impulse and frequency responses are shown in Figure(2).

Please note that the higher-order B-splines don’t have zero-crossing at integer x, so the interpolated curve will not necessarily go thorough the points.

**
Figure 2)
**The
impulse and frequency responses of the first three odd-order B-splines,
including the linear interpolator.

It is seen that the signal components which have frequencies of multiples of 2π, will be attenuated heavily but those which have frequencies of multiples of π, will not have much attenuation. Besides, we observe that higher order means smaller sidelobes and as a result, better interpolating. So, when we use more neighbor points to estimate the value of the signal at the point we are looking for, i.e. we have a higher order polynomial, we will get a better result and it is seen that the associated filter would have smaller sidelobes.

Lagrange polynomials are forced to go thorough a
number of points. For instance, the 4-point Lagrange interpolator polynomial is
formed so that it goes through all of the four neighboring points, and the
middle section is used [5]. The 1^{st}-order (2-point) Lagrange
interpolator is the linear interpolator, which was already presented as a part
of the B-spline family. The order of the Lagrange polynomial is always one less
than the number of points. The third and 5^{th}-order Lagrange impulse
responses are:

(8)
*4-point, 3 ^{rd}-order Lagrange
impulse response*

**
Figure 3)
**The
impulse and frequency responses of the first three odd-order Lagrange, including
the linear interpolator.

It is seen that the sidelobes don’t get much lower when the order increases, but the passband behavior improves.

Hermite interpolation may be considered that type of osculatory interpolation where the derivatives up to some degree, are not only supposed to be continuous everywhere but are also required to assume prespecified values at the sample points. Hermite interpolators have a continuous first derivative.

Here are the impulse responses of some Hermite interpolators.

(10) *4-point,
3 ^{rd}-order Hermite impulse response.* (11)

*This
is also known as the Catmull-Rom spline.*

(12)* 6-point, 5 ^{th}-order Hermite impulse
response*

**
Figure 4)
**The impulse and frequency
responses of the 3^{rd}-(4-point and 6-point) and 5^{th}-order
Hermite.

An interesting feature in the Hermite frequency responses is that the largest sidelobes (compare for example to the frequency responses of Lagrange interpolators) have been “punctured”.

The definition of the 2^{nd}-order-osculating
interpolator is the same as of the Hermites, but in this case, the second
derivatives are also matched with the Lagrangians. The order of the
interpolation polynomial must be at least 5 [5].

Here are the impulse responses of some 2^{nd}-order-osculating
interpolators.

(10) *4-point,
5 ^{th}-order 2^{nd}-order-osculating impulse response*

**
Figure 5)
**The impulse and frequency
responses of the 2^{nd}-order-osculating(4-point and 6-point).

We know that if a signal has finite energy, the
minimum sampling rate is equal to two samples per period of the highest
frequency component of the signal. Specifically, if the highest frequency
component of the signal is *B* hertz, then the signal, *x*(*t*),
can be recovered from the samples by [8]:

(12)

The frequency *B* is also referred to the
signal’s bandwidth. The series of (12) is called the “*cardinal series*”[8]
but has many other names, including the *Whittaker-Shannon sampling theorem*
[Goodman (1968)], the *Whittaker-Shannon-Kotel’nikov-Kramer sampling theorem*
[Jerri (1977)].

Exact interpolation using the cardinal series assumes that (a) the values of the samples are known exactly, (b) the sample locations are known exactly and (c) an infinite number of terms are used in the series. Deviation these requirements results in interpolation error due to (a) data noise (b) jitter and (c) truncation respectively.

Suppose that the signal we sample is corrupted by real additive zero-mean wide-sense-stationary noise, . Then, instead of sampling the deterministic bandlimited signal, we would be sampling the signal . From these samples, we form the series

(13)

Therefore, is the stochastic process generated by the samples of alone and it is independent of the signal. Note that since is zero-mean, so is . Hence, expecting both sides of (13) leads us to the desirable conclusion that

(14)

A meaningful of the cardinal series’ noise
sensitivity is the *interpolation noise level,* which is the noise
variance, .
It can easily be shown[8] that .
That is, the cardinal series interpolation results in a noise level equal to
that of the original signal before sampling [Bracewell]. In many cases, one can
reduce the interpolation noise level by oversampling and filtering. If , for
example, we place *z*(*t*) in (13) thorough a filter that is unity
for and zero otherwise, then the output is

(15)

It is proved that [8]: .

Jitter occurs when samples are taken near to but not exactly at the desired sample locations. The cardinal series interpolation of jittered samples yields a biased estimate of the original signal. Although the bias can be corrected with an inverse filter, the resulting interpolation noise variance is increased [8].

The cardinal series requires an infinite number of
terms to exactly interpolate a bandlimited signal from its samples. In
practice, only a finite number of terms can be used. We can write the truncated
cardinal series approximation of *x*(*t*) as

(16)

The error resulting from using a finite number of terms is referred to as truncation error: . It can be shown that [Papoulis (1966)]:

(17)

Where *E* is the energy of and is the energy of .

The main difference between them is that in Shannon sampling theorem, we will assume that the given data is band limited to some band of frequencies and develop schemes which are optimal on this basis, while a numerical analyst typically assumes that his data is samples of polynomials (or very nearly so) and develops schemes to minimize the resulting error.

Please note that the sources of error discussed in 3.2.1.1 and 3.2.1.2 exist for piecewise polynomials interpolators as well. Besides, for piecewise polynomials interpolators, since to compute the value of each lost point, we just use the values of a finite number of its neighbor points, we will have another source of error.

Another approach to interpolation is to consider
data as realization of a random process at given spatial locations and to make
inferences on the observed values of the process by means of a statistically
optimal predictor. Here, “optimal” refers to minimizing the mean-squared
prediction error, which requires a model of the covariance between the data
points. So, it does not have an interpretation in frequency domain. This
approach to interpolation is generally known as *Kriging* [10]. In
Kriging, a crucial role is played by the *variogram*, a diagram of the
variance of the difference between the measurements at two input locations [11].
A unique feature of Kriging is that it provides an estimation of the error at
each interpolated point, providing a measure of confidence in the modeled
surface. The effectiveness of Kriging depends on the correct specification of
several parameters that describe the semivariogram and the model of the drift
(i.e., does the mean value change over distance). Because Kriging is a robust
interpolator, even a naive selection of parameters will provide an estimate
comparable to many other grid estimation procedures. The trade-off for
estimating the optimal solution for each point by Kriging is computation time.
Given the additional trial and error time necessary to select appropriate
parameters, Kriging should be applied where best estimators are required, data
quality is good, and error estimates are essential [16]. Besides, In this
method, samples need not be uniformly sampled.

A random process *Z*(.) can be described by where *D* is a fixed subset of and *Z*(*s*) is a random function at
location .
There are several types of Kriging, such as “Ordinary Kriging”, “Simple
Kriging”, and “Universal Kriging”, but I limit this subsection to *Ordinary
Kriging*, which makes the following two assumptions [11]:

(i)
The *model assumption* is that the
random process consists of a constant and an error term : (18)

(ii)
The *predictor assumption* is that
the predictor for the point - denoted by - is a weighted linear function of all the
observed output data:

(19)

To select the weights in (19),
the *criterion* is minimal mean-squared prediction error, ,
defined as

(20)

To minimize (20)
given (19), let *m* be the
Lagrangian multiplier ensuring ,
then we can write the prediction error as

(21)

Minimizing (21) gives us the optimal [11]:

(21) and

where denotes the vector and denotes the matrix whose element is .

The optimal weights (21) give the minimal mean-square prediction error. (20) becomes

(22)

However, in (21) and (22) is unknown. Usually it is estimated by

(23)

Where denotes the number of distinct pairs in . The estimator (23) is unbiased, if the process is indeed second-order stationary [12].

Irregular or non-uniform sampling constitutes another whole area of research. There are essentially two strategies for finding a solution [13].

Here, we consider the same kind of *shift-invariant*
spaces as in the uniform case and fit the model to the measurements [14].

Here, we define new basis functions (or spaces) that are better suited to the non-uniform structure of the problem, for instance, non-uniform splines [4,15].

I used MATLAB to simulate some of these methods. In
simulation, we have to note that since in MATLAB any signal is stored and
treated as a data vector, when we are dealing with continuous signals or
discrete-time signals with different sampling rate, we have to have the same
number of data points in each time unit for all the signals, so if they have
different sampling rates, we can fill the missing data in each unit time with
zeros. So, after this zero padding, all the signals will have the same number
of data in each time unit and we can work with them easily. For example, **Script 1** shows the MATLAB script I wrote to
do the simulation for piecewise polynomial interpolators. In this script, I
tried to recover the signal in the time interval .
First, I took some samples from the original signal, then I down-sampled the
result discrete-time signal and used the result data points with some different
piecewise polynomial interpolators to recover the original signal by convolving
the “interpolator” and “known data points”.
Since we are using the “convolution operation” and both our signals have
finite duration, the (NOP/2).Fs points in
the beginning and the end of the recovered signals are not correct and as you
can see in the MATLAB script, I got rid of them (NOP:
the number of points used in the interpolator, Fs: sampling frequency of the
original signal). Then, I computed the average error for the data points
I estimated, and the real data points that I had before the down-sampling. The
error measurement I used, is where is the real signal data point and is the estimated data point and *N* is
the number of data points we estimated. Please note that we are recovering an
analog signal, so we have to use “integral” instead of “sum” to compute the
average error, but as I said, all signals in MATLAB are discrete, so I used the
“sum”. **Figure 6** shows some of the
simulation results.

t=[0:0.35:10]; x = t.*(sin(t))+ 3.*t.^(1/2);

[w,n] = dnsample(x,2); [vv,n] = upsample(x,Fs/2); [v,n] = upsample(w,Fs);

y=conv(v,f); yy = y((NOP/2)*Fs+1:size(y,2)-(NOP/2)*Fs); % NOP: the number of points used in interpolator “f”

figure(1); hold on; tt = n/(max(n)/max(t)); plot(tt,yy,'b-');

stem(tt(find(vv~=0)),vv(find(vv~=0)),'r:');stem(tt(find(v~=0)),v(find(v~=0)),'g-.');

% vv: Real values, v: Known values, yy: Interpolated function

function [v,n] = upsample(x,L);

lenx = length(x); x = reshape(x,1,lenx); lenv = lenx*L-(L-1); n = [0:lenv-1]; v = zeros(1,lenv); v(1:L:end) = x;

function [w,n] = dnsample(x,M);

lenx = length(x); x = reshape(x,1,lenx); w = x(1:M:end); lenw = length(w); n = [0:lenw-1];

**Script1 ***-
The MATLAB script I used to simulate some of the piecewise polynomial
interpolators*

* *

* *

**
Figure 6)
**Some of the simulation
results.

In this paper, I talked about some of the methods being used to interpolation, but there are lot’s of other methods that I did not mention here. People with different backgrounds and different point of views have worked on these methods and found them useful in different areas. So, it is impossible to give a general comparison of them and say which one is the best. It depends on the case we are dealing with, to choose a method and I don’t think there will be any guarantee that the chosen method be the best possible one. Some people have compared some of these methods in their papers. For instance, Fraser says [17] “unlike small-kernel convolution methods, which have poor accuracy anywhere near the Nyquist limit, the accuracy of the FFT method is maintained to high frequencies, very close to Nyquist”. But, we shouldn’t forget that both those types of methods have been highly improved since then.

p.s.: I was one asked once about the accuracy of a naive approach using the shift property of DFT and a non-integer shift. This method mathematically is wrong and does not work, because in the proof of the shift property of DFT, it is an essential assumption that the shift is integer [2, 3]. I tried to test that in practice, and noticed it works (to some extent!) just in a case that the non-integer shift is very small in comparison with the DFT period (less than %10) and the original signal is smooth. We can say that in this case, the DFT coefficients which actually are the uniformly-spaced samples of the unit circle, have moved very little and so the IDFT gives us the values of the points that are very close to the original samples. But, it is not useful, because when we know that the signal is smooth and we want to find the values of the samples very close to the known samples, it is obvious that the result values will be very close the values of the original samples anyway!

**Bibliography**

[1] J.S. Lim and A.V.
Oppenheim, “*Advance Topics in Signal Processing,*” Prentice Hall, 1998.

[2] J.G. Proakis and
D.G. Manolakis, “*Digital signal processing: principles, algorithms and
applications,*” Third Edition, Prentice Hall, 1996.

[3] A. V. Oppenheim and R. W. Schafer, “*Discrete
Time Signal Processing,*” Prentice Hall, 1989.

[4] Carl Deboor, “*A practical guide to splines,*”
Springer-Verlag, 1978.

[5] Olli Niemitalo, “*Polynomial Interpolators for
High-Quality Resampling of Oversampled Audio,*” http://www.student.oulu.fi/~oniemita/DSP/INDEX.HTM

[6] M. Unser, A. Aldroubi and M. Eden, ”*B-Spline
Signal Processing: Part I-Theory,*” IEEE Transactions on Signal Processing,
vol. 41, no. 2, pp. 821-833, February 1993.

[7] C. de Boor, “*A practical Guide to Splines,*”

[8] Robert J. Marks II, “*Introduction to Shannon
Sampling Theory,*”

[9] R.W. Schafer and L.R. Rabiner, “*A Digital
Signal Processing approach to interpolation,*” Proc. IEEE, vol. 61, pp.
692-702, June 1973.

[10] Erik Meijering, “*A Chronology of
Interpolation: Form Ancient Astronomy to Modern Signal and Image Processing,*”
Proceedings of the IEEE, Vol. 90, pp. 319-342, March 2002.

[11] Win C.M. van Beers and Jack P.C. Kleijnen, “*Kriging
for Interpolation in Random Simulation,*” Discussion paper, ISSN 0924-7815.

[12] Cressie, “*Statistics for Spatial Data,*”

[13] Michael Unser, “*Sampling-50 Years After Shannon,*” Proceedings of the IEEE, Vol. 88, No. 4,
April 2000.

[14] J.J. Benedetto and W. Heller, “*Irregular
Sampling and the Theory of Frames,*” Math. Note, Vol. 10, pp. 103-125, 1990.

[15] L.L. Schumaker, “*Spline functions: Basic
Theory*”,

[16] “*Kriging Interpolation*”, http://lazarus.elte.hu/hun/digkonyv/havas/mellekl/vm25/vma07.pdf

[17] Donald Fraser, “*Interpolation by the FFT
Revisited An Experimental Investigation,*” IEEE
Transactions on Acoustics, Speech, and Signal Processing, Vol. 37, No. 5, pp.
665-675, 1989.

[18] R.W. Schafer and L.R. Rabiner, “*A digital
Signal Processing approach to interpolation*,” Proc. IEEE, Vol. 61, pp.
692-702, June 1973.