Electronics-Related.com
Blogs

Interpolation Basics

Neil RobertsonAugust 20, 201917 comments

This article covers interpolation basics, and provides a numerical example of interpolation of a time signal.  Figure 1 illustrates what we mean by interpolation.  The top plot shows a continuous time signal, and the middle plot shows a sampled version with sample time Ts.  The goal of interpolation is to increase the sample rate such that the new (interpolated) sample values are close to the values of the continuous signal at the sample times [1].  For example, if we increase the sample rate by the integer factor of four, the interpolated signal is as shown in the bottom plot.  The time between samples has been decreased from Ts to Ts/4.

The simplest technique for interpolation is linear interpolation, in which you draw a straight line between sample points, and compute the new samples that fall on the line.  However, in many cases, linear interpolation is not accurate enough.  Other techniques involve fitting polynomials to the time function.  Here, we’ll show the power of approaching the problem in the frequency domain.

Figure 1.  Top:  Continuous signal       

Middle:  Signal with sample time Ts  

Bottom:  Interpolated signal with sample time Ts/4  


Before we describe how interpolation works, let’s look at two sampled signals and their spectra.  In the top row of Figure 2, the time signal u is sampled at 400 Hz, and its magnitude spectrum is shown on the right.  In the bottom row, the time signal x consists of every fourth sample of u, and thus has sample rate of 100 Hz.  The spectrum is again shown on the right, with the image spectrum above 100 Hz shown by a dashed line.  Note the amplitude of the spectrum is ¼ that that of u.

Now let’s look at an example of interpolation by an integer factor of four.  For the signal x just described, interpolation by four should result in the signal u.  Matlab code for the example interpolator is provided at the end of the article.

If we simply increase the sample rate of x from 100 Hz to 400 Hz, we get the signal shown in Figure 3.  Here we’ve just inserted three zeros between each of the original samples of x, a process called upsampling.  Upsampling includes scaling the amplitude of x by the upsampling ratio of 4, which gives a maximum spectrum magnitude of 1, instead of 1/4.  The spectrum of the upsampled signal xup is shown on the right.  Given that xup sample values have the same time spacing as x, the spectrum is the same as that as x, except its amplitude is a factor of 4 larger.

Now, if we compare the spectrum of xup to that of u, they match from 0 to 50 Hz, but xup has additional image spectra centered at 100 and 200 Hz.  So, it seems that we should be able to approximate u by low-pass filtering xup to attenuate the images.  If this works, we will have interpolated x by four.  A block diagram of the interpolator is shown in Figure 4.  In our case L = 4.

Figure 2.  Top:  sampled signal u with fs = 400 Hz and its spectrum (linear amplitude scale)

Bottom:   sampled signal x with fs = 100 Hz and its spectrum

Figure 3.  Upsampled signal xup with fs = 400 Hz and its spectrum (linear amplitude scale)

Figure 4.  Interpolator block diagram (for our example, L = 4)             


To design the low-pass filter, we’ll look at the dB spectrum of xup, as shown in the top of Figure 5.  We need a filter that passes the desired spectrum and attenuates the undesired images.  Thus, we have the following passband and stopband:

Passband:  0 to 36 Hz

Stopband:  66 to 135 Hz and 166 to 200 Hz

We leave the ranges from 50 to 66 Hz and 135 to 166 Hz unspecified, since the amplitude of xup’s spectrum is low in these ranges.  The Matlab code at the end of this article includes synthesis of a 41-tap FIR interpolation filter.  The filter response is shown in the middle of Figure 5.  The spectrum of xinterp is shown in the bottom plot.

The time signal xinterp is shown on the left of Figure 6, with the linear-amplitude spectrum on the right.  As hoped, xinterp resembles the plot of u in Figure 2.   We can compute the error between xinterp and u as:

% error = 100*(xinterp – u)/max(u)

The error is plotted in Figure 7.   Reflecting on how we achieved this result, it is worth noting that we performed interpolation in time by applying knowledge of the input signal’s frequency spectrum.

Figure 5.  Top:  Spectrum of xup            

Middle:  Interpolation low-pass filter response          

Bottom:  Spectrum of xinterp          



Figure 6.  Left:  Interpolator output xinterp (compare to u in Figure 2)

Right:  spectrum of xinterp (linear amplitude scale)

Figure 7.  Interpolation percent error                 


Matlab code demonstrating interpolation by four

The following code uses the same signal names as the text:

u                        signal sampled at 400 Hz

x                        interpolator input signal (formed from every 4th sample of u)

x_up                  upsampled version of x, fs = 400 Hz

x_interp             interpolation filter output, fs = 400 Hz

The signal u is a Chebyshev window function of length 32 with -70 dB sidelobes [2,3].  Other functions could be used, for example a Blackman window.  As previously discussed, the calculation of xup includes scaling by 4.  The interpolator output x_interp has length 72.  Samples 21:52 are used as an approximation of the signal u.

The interpolation filter has fs = 400 Hz and is synthesized using the Parks-McClellan algorithm (Matlab function firpm).  The coefficients are plotted in Figure 8, and the filter’s frequency response is shown in the center plot of Figure 5.  Note that while this design approach is straightforward, it does not result in the most efficient interpolator.  One alternative approach would cascade two interpolate-by-two sections, as shown in Figure 9.  This allows the use of halfband interpolation filters [4].  For a discussion of efficient approaches to interpolator design, see [5].


%interp_demo.m  8/17/2019 Neil Robertson
% Demonstrate interpolation by 4
fs= 400;                % Hz sample rate
p= chebwin(32,70)';     % signal = Chebyshev window
u= p/sum(p);            % normalize for sum = 1
x= [u(1:4:end)];        % interpolator input signal, fs= 100
% upsampling
x_up= zeros(1,32);
x_up(1:4:32)= 4*x(1:8); % upsampled signal
%
% interpolation filter using Parks-McClellan algorithm
fn= fs/2;
f= [0 36 66 135 166 199]/fn;    % frequency vector
a= [1 1 0 0 0 0];               % amplitude goal vector
Ntaps= 41;
b= firpm(Ntaps-1,f,a);          % synthesize filter coeffs
b= round(b*2^13)/2^13;          % fixed point coeffs
%
x_interp= conv(x_up,b);         % filter x_up
interp_error= 100*(x_interp(21:52) - u)/max(u);
%
[Xup,f]= freqz(x_up,1,256,fs);      % spectrum of x_up
Xinterp= freqz(x_interp,1,256,fs);  % spectrum of x_interp
%
%
% plotting
subplot(211),plot(u,'.','markersize',10),grid
axis([1 32 0 .1]),xticklabels({}),text(23,.06,'u')
subplot(212),plot(x,'.','markersize',10),grid
axis([1 9 0 .1]),xticklabels({}),text(6.5,.06,'x'),figure
stem(b),grid
axis([1 41 -.05 .3]),title('Interpolation Filter Coefficients b')
figure
subplot(211),stem(x_up),grid
axis([1 32 0 .35]),xticklabels({}),text(23,.24,'x_{up}')
subplot(212),plot(x_interp(21:52),'.','markersize',10),grid
axis([1 32 0 .1]),xticklabels({}),text(23,.06,'x_{interp}'),figure
plot(interp_error,'.','markersize',10),grid
axis([1 32 -1 1]),xlabel('sample'),ylabel('percent error')
title('Interpolation Percent Error'),figure
subplot(211),plot(f,abs(Xup)),grid
xlabel('Hz')
axis([0 fs/2 0 1]),title('Spectrum of x_{up} (linear amplitude scale)')
subplot(212),plot(f,abs(Xinterp)),grid
xlabel('Hz')
axis([0 fs/2 0 1]),title('Spectrum of x_{interp} (linear amplitude scale)')


Figure 8.  Interpolate-by-four filter coefficients.     




Figure 9.  Interpolation by four using two Interpolate-by-two stages.



References

  1. Lyons, Richard G., Understanding Digital Signal Processing, 3rd Ed., Prentice-Hall, 2011, section 10.4.
  2. The Mathworks website, https://www.mathworks.com/help/signal/ref/chebwin.html
  3. Wikipedia, “Window Function”, https://en.wikipedia.org/wiki/Window_function
  4. Robertson, Neil, “Simplest Calculation of Halfband Filter Coefficients”, https://www.dsprelated.com/showarticle/1113.php
  5. Harris, Fredric J., Multirate Signal Processing for Communication Systems, Prentice Hall, 2004, Ch 7.



Neil Robertson               August, 2019





[ - ]
Comment by jimi75September 9, 2019

Hi,

Thank you for this interesting post. I tried to implement it in python. If anyone was interested just look at this link: https://colab.research.google.com/github/Strato75/.... Please, let me know for bugs or any other issue.

Giovanni

[ - ]
Comment by neiroberSeptember 9, 2019

Giovanni,

Thanks for the Python code.

Neil


[ - ]
Comment by omersayliSeptember 9, 2019

Off the topic a little bit but this Python & math and dsp packages combined with Jupyter notebooks / Google co-lab type running environments is really a big power and facility.. 

Thanks for the code

[ - ]
Comment by woodpeckerSeptember 12, 2019

Hi Neil,

Thanks for the interesting article. I have succesfully run the Matlab code in Octave ver 5.1, by simply replacing the firpm function with the remez function.

remez is contained within the 'signal' add-on package for Octave.

[ - ]
Comment by neiroberSeptember 12, 2019

Thanks.  It certainly confused matters when Matlab changed the name of their Parks-McClellan function from remez to firpm.


[ - ]
Comment by NewhufSeptember 14, 2019

Hi, 

I want to simulate a bidirectional full-duplex communication network and I don’t know how to go about it in MATLAB. Please do you have material which could be of help. Thank you

[ - ]
Comment by DanielFZLMay 4, 2021

Thank you very much for the article.

I was trying to implement a Farrow structure using Lagrange interpolation. The output signal's (the interpolated signal) spectrum was compressed as if the fundamental tone was close to zero. Also, there were some spurs that could be seen from fundamental frequency till fs/2. 

Is it true that after interpolation, the spectrum will be compressed in the frequency domain?

[ - ]
Comment by neiroberMay 4, 2021

Hi Daniel,

No, the frequency at the output should be the same as that at the input.  Note, however, that for a sine of frequency f0, if you interpolate by two, f0/fs_out is 1/2 of f0/fs_in, where fs_in is input sample rate and fs_out is output sample rate.  For example, if f0 = 5 Hz, fs_in= 20 Hz, and fs_out= 40 Hz, then f0/fs_in = 1/4 and f0/fs_out= 1/8.

There is a discussion of Farrow interpolators in Michael Rice's book "Digital Communications, A Discrete-Time Approach".

regards,

Neil

[ - ]
Comment by DanielFZLMay 6, 2021

Thank you for the response and for the book.

The problem is that when representing the PSD, the fundamental frequency is very close to zero. and I suspect I am doing something wrong about that or it is the nature of the interpolation. 

Also, I cannot detect the spurs of the signal in low frequencies(maybe I am trying to find out the SFDR). since my PSD representation is not perfect the spurs should really show up to be able to detect them using the findpeak function of Matlab.

with regards,

[ - ]
Comment by vfreyMay 19, 2025

Thanks for the article, I definitely had the wrong definition of interpolation (or at least we did not have the same in our team). I have a question regarding my process, here is the context first:

I want to filter out a cadence (like pedal cadence on a bike), from a sensor proving information every 125ms seconds (8Hz). There are some message drops. My approach was to upsample (generate sample at 16Hz). I would linear interpolate between each new samples - as opposed to a DSP interpolation where it is zeros. And then I would apply an IIR or FIR (still not clear which one would apply best). 

My question would be: is the linear interpolation a good approach here? Can I achieve the same with zero?

[ - ]
Comment by neiroberMay 19, 2025

Hi,

I don't really understand your application, but I can make some general comments.  (If you want to provide more information, I may be able to give a more detailed answer).

Linear interpolation can be performed by using upsampling and lowpass filtering as shown in Figure 4.  In your stated case, you would upsample by 2 and then filter with a linear interpolate-by-2 FIR filter with coefficients:

b = [1/4 1/2 1/4]

The frequency response magnitude of this linear interpolate by 2 filter is shown in the first plot.  

In the time domain, consider the signal x= [1 2 3].  If you upsample by 2, you get u = [1 0 2 0 3 0];  If you multiply this sequence by two and convolve it with b:

y = conv(2*u,b),

you get the interpolator output shown in the 2nd plot.

regards,

Neil

freq_response_25731.png

time_response_91434.png


[ - ]
Comment by vfreyMay 21, 2025

Thanks for the time you spent on answering my message Neil. 

My application is a rider cadence sensor on a bicycle, which transmit data over a wireless protocol. The messages are supposed to come at a fixed interval (4Hz), but we observed some jitter and lost messages.

The receiving device is filtering out the incoming data from the messages. The issue I have, and if I want to apply DSP, I need discrete time, which are meant to be sampled at a constant rate. However, due to the jitter and the low frequency, I had to come up with a - debatable - solution:

  • I have defined a fix sample interval, let's say 25ms (40 HZ). The incoming messages are supposed to arrive every 250ms, but due to the time jitter, I artificially "re-align" the sample to a discrete-time, and this is where I used linear interpolation, to evaluate the cadence for the closest discrete time.

This is essentially how I would build my samples in the time-domain, from 4 samples a second (in the ideal case, there might be some dropped messages), I create 40 samples, with their value interpolated. This is really rudimentary, I was wondering if there is a better approach.


Victor

[ - ]
Comment by neiroberMay 21, 2025

Victor,

Basic questions:  why can you not use the cadence sensor in the manner intended by the manufacturer?  What is the transmission system?  (Bluetooth, etc).  What is the data format?  Is the message interval fixed at 0.25 seconds?

Are the transmitter and receiver part of a system made by a single vendor?

-- Neil

[ - ]
Comment by neiroberMay 22, 2025

Update:  It looks like the answer to this problem was to run the jittery data through a FIFO.


To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.

Please login (on the right) if you already have an account on this platform.

Otherwise, please use this form to register (free) an join one of the largest online community for Electrical/Embedded/DSP/FPGA/ML engineers: