# What is an algorithm to re-sample from a variable rate to a fixed rate?

I have a sensor that reports its readings with a time stamp and a value. However, it does not generate readings at a fixed rate.

I find the variable rate data difficult to deal with. Most filters expect a fixed sample rate. Drawing graphs is easier with a fixed sample rate as well.

Is there an algorithm to resample from a variable sample rate to a fixed sample rate?

• This is a cross post from programmers. I was told this is a better place to ask. programmers.stackexchange.com/questions/193795/… – FigBug Apr 3 '13 at 20:20
• It tries to sample every 5ms or 10ms. But it is a low priority task, so it may get missed or delayed. I have the timing accurate to 1 ms. The processing is done on the PC, not in real time, so a slow algorithm is ok if it's easier to implement. – FigBug Apr 3 '13 at 21:08
• Have you took a look at a fourier reconstruction? There is a fourier transform based on unevenly sampled data. The usual aoproach is to transform a fourier image back to evenly sampled time domain. – mbaitoff Apr 4 '13 at 1:42
• Do you know any characteristics of the underlying signal that you are sampling? If the irregularly-spaced data is still at a reasonably high sample rate compared to the bandwidth of the signal being measured, then something simple like polynomial interpolation to an evenly-spaced time grid might work fine. – Jason R Apr 4 '13 at 13:11

The simplest approach is to do some kind of spline interpolation like Jim Clay suggests (linear or otherwise). However, if you have the luxury of batch processing, and especially if you have an overdetermined set of nonuniform samples, there's a "perfect reconstruction" algorithm that's extremely elegant. For numerical reasons, it may not be practical in all cases, but it's at least worth knowing about conceptually. I first read about it in this paper.

The trick is to consider your set of nonuniform samples as having already been reconstructed from uniform samples through sinc interpolation. Following the notation in the paper:

$$y(t) = \sum_{k=1}^{N}{y(kT)\frac{\sin(\pi(t - kT)/T)}{\pi(t - kT)/T}} = \sum_{k=1}^{N}{y(kT)\mathrm{sinc}(\frac{t - kT}{T})}.$$

Note that this provides a set of linear equations, one for each nonuniform sample $y(t)$, where the unknowns are the equally-spaced samples $y(kT)$, like so:

$$\begin{bmatrix} y(t_0) \\ y(t_1) \\ \cdots \\ y(t_m) \end{bmatrix} = \begin{bmatrix} \mathrm{sinc}(\frac{t_0 - T}{T}) & \mathrm{sinc}(\frac{t_0 - 2T}{T}) & \cdots & \mathrm{sinc}(\frac{t_0 - nT}{T}) \\ \mathrm{sinc}(\frac{t_1 - T}{T}) & \mathrm{sinc}(\frac{t_1 - 2T}{T}) & \cdots & \mathrm{sinc}(\frac{t_1 - nT}{T}) \\ \cdots & \cdots & \cdots &\cdots \\ \mathrm{sinc}(\frac{t_m - T}{T}) & \mathrm{sinc}(\frac{t_m - 2T}{T}) & \cdots & \mathrm{sinc}(\frac{t_m - nT}{T}) \end{bmatrix} \begin{bmatrix} y(T) \\ y(2T) \\ \cdots \\ y(nT) \end{bmatrix}.$$

In the above equation, $n$ is the number of unknown uniform samples, $T$ is the inverse of the uniform sample rate, and $m$ is the number of nonuniform samples (which may be greater than $n$). By computing the least squares solution of that system, the uniform samples can be reconstructed. Technically, only $n$ nonuniform samples are necessary, but depending on how "scattered" they are in time, the interpolation matrix may be horribly ill-conditioned. When that's the case, using more nonuniform samples usually helps.

As a toy example, here's a comparison (using numpy) between the above method and cubic spline interpolation on a mildly jittered grid: (Code to reproduce the above plot is included at the end of this answer)

All that being said, for high-quality, robust methods, starting with something in one of the following papers would probably be more appropriate:

A. Aldroubi and Karlheinz Grochenig, Nonuniform sampling and reconstruction in shift-invariant spaces, SIAM Rev., 2001, no. 4, 585-620. (pdf link).

K. Grochenig and H. Schwab, Fast local reconstruction methods for nonuniform sampling in shift-invariant spaces, SIAM J. Matrix Anal. Appl., 24(2003), 899- 913.

--

import numpy as np
import pylab as py

import scipy.interpolate as spi
import numpy.random as npr
import numpy.linalg as npl

npr.seed(0)

class Signal(object):

def __init__(self, x, y):
self.x = x
self.y = y

def plot(self, title):
self._plot(title)
py.plot(self.x, self.y ,'bo-')
py.ylim([-1.8,1.8])
py.plot(hires.x,hires.y, 'k-', alpha=.5)

def _plot(self, title):
py.grid()
py.title(title)
py.xlim([0.0,1.0])

def sinc_resample(self, xnew):
m,n = (len(self.x), len(xnew))
T = 1./n
A = np.zeros((m,n))

for i in range(0,m):
A[i,:] = np.sinc((self.x[i] - xnew)/T)

return Signal(xnew, npl.lstsq(A,self.y))

def spline_resample(self, xnew):
s = spi.splrep(self.x, self.y)
return Signal(xnew, spi.splev(xnew, s))

class Error(Signal):

def __init__(self, a, b):
self.x = a.x
self.y = np.abs(a.y - b.y)

def plot(self, title):
self._plot(title)
py.plot(self.x, self.y, 'bo-')
py.ylim([0.0,.5])

def grid(n): return np.linspace(0.0,1.0,n)
def sample(f, x): return Signal(x, f(x))

def random_offsets(n, amt=.5):
return (amt/n) * (npr.random(n) - .5)

def jittered_grid(n, amt=.5):
return np.sort(grid(n) + random_offsets(n,amt))

def f(x):
t = np.pi * 2.0 * x
return np.sin(t) + .5 * np.sin(14.0*t)

n = 30
m = n + 1

# Signals
even   = sample(f, np.r_[1:n+1] / float(n))
uneven = sample(f, jittered_grid(m))
hires  = sample(f, grid(10*n))

sinc   = uneven.sinc_resample(even.x)
spline = uneven.spline_resample(even.x)

sinc_err   = Error(sinc, even)
spline_err = Error(spline, even)

# Plot Labels
sn = lambda x,n: "%sly Sampled (%s points)" % (x,n)
r  = lambda x: "%s Reconstruction" % x
re = lambda x: "%s Error" % r(x)

plots = [
[even,       sn("Even", n)],
[uneven,     sn("Uneven", m)],
[sinc,       r("Sinc")],
[sinc_err,   re("Sinc")],
[spline,     r("Cubic Spline")],
[spline_err, re("Cubic Spline")]
]

for i in range(0,len(plots)):
py.subplot(3, 2, i+1)
p = plots[i]
p.plot(p)

py.show()

• Nice method and code. But for $t_j = jT$ with a few dropouts (e.g. [0 1 - 3 4 5 - 7 8] T), which I think is the OPs question, aren't the sinc s in the matrix all 0 ? Sure there are ways to fix that, but. – denis May 17 '13 at 14:56
• As I understand it, the OP's question is about samples that get dropped and/or delayed. This method is basically just an overdetermined system of equations, so dropped samples only show up as unknowns (not as data points with a value of 0). Or maybe that's not what you're asking? – datageist May 17 '13 at 20:09
• What happens if the $t_j$ are all integers (T=1) ? Say we have data points [$j, y_j$] for $j \in J$, a set of non-zero integers e.g. {-1 1} or {-2 -1 1 2}; isn't the interpolated $y_0 = 0$, regardless of the $y_j$ -- or have I missed something ? – denis May 18 '13 at 16:20
• If the sample rates are exactly identical (w/ missing points), then the interpolation matrix will be sparse (because each output depends on only one input). In general, the average sample rate of the nonuniform samples needs to be greater than the uniform reconstruction rate. In other words, you'd need to reconstruct at a lower rate to "fill in the gaps" (T>1 for your example). I see your point though. – datageist May 19 '13 at 7:55
• Answers like this are pure gold. – Ahmed Fasih Jun 7 '15 at 3:37

This sounds like a problem of asynchronous sample rate conversion. To convert from one sample rate to another, we can compute the continuous time representation of the signal by performing sinc interpolation, then resample at our new sample rate. What you are doing is not much different. You need to resample your signal to have sample times that are fixed.

The continuous time signal can be calculated by convolving each sample with a sinc function. Since the sinc function continues indefinately, we use something more practical like a windowed sinc of a practical finite length. The tricky part is that because your samples move about in time, a sinc with a different phase offset may need to be used for each sample when resampling.

Continuous time signal from sampled signal:

$x(t) = \sum\limits_{n=-\infty}^\infty x[n]sinc({t-nT_s\over T_s})$

where $T_s$ is your sample time. In your case, however, your sample time is not fixed. So I think you need to replace it with the sample time at that sample.

$x(t) = \sum\limits_{n=-\infty}^\infty x[n]sinc({t-nT_s[n]\over T_s[n]})$

From this you can resample the signal:

$y[n] = x(nT_{ns}$)

where $T_{ns}$ is the desired sample time.

Putting it all together you get:

$y[m] = \sum\limits_{n=-\infty}^\infty x[n]sinc({mT_{ns}-nT_s[n]\over T_s[n]})$

Since this is not causal or tractable, the sinc function can be replaced with a function of finite support and the limits of summation adjusted accordingly.

Let kernel(t) be a windowed sinc or other similar function of length 2k then:

$y[m] = \sum\limits_{n=-k}^k x[n]kernel({mT_{ns}-nT_s[n]\over T_s[n]})$

I hope this helps..., but I may have made a mistake along the way and it might be a little math intensive. I would recommend researching sample rate conversion for more information. Maybe somebody else here could give a better explanation or solution as well.

• Using a sinc function to reconstruct a continuous version of a signal requires that the samples be equally spaced, so the sinc function will have to adapt to arbitray sample spacing. Could be rather difficult to implement. – user2718 Apr 5 '13 at 11:43
• yeah, this would not be very efficient to do exactly as seen here. It would require calculating new kernel coefficients for each different sample time. A collection of several kernels could be computed, however, and the time quantized to one of these. There would be a performance hit relative to the quantization error. – Jacob Apr 5 '13 at 19:32
• You could also pre-compute a single sinc lookup table and interpolate between the points of that lookup table. – jms Jun 6 '18 at 17:46

I think Jacob's answer is very workable.

An easier method that's probably not quite as good in terms of introducing distortion is to do polynomial interpolation. I would use either linear interpolation (easy, not as good signal performance-wise) or cubic splines (still not too hard, better signal performance) to produce samples at any time you want from your arbitrary time samples.

• You answer seems a lot easier to implement than Jacob's, so I've gone with it first. It seems to be working, but I haven't run a full set of tests yet. – FigBug Apr 5 '13 at 0:20
• @FigBug -If you have time, add a comment with you're final solution. – user2718 Apr 6 '13 at 0:32

(A month later) there are two main choices for any interpolation method:
1) the number of data points nearest the missing point to use, $Nnear$ 2 4 6 ...
2) the class of basis functions to use: linear, polynomial, sine-cosine (Fourier), piecewise cubic (B-spline or interpolating spline), sinc-like ...
(Choice 0 is whether to use somebody else's method and code, or do-it-yourself.)

Fitting a straight line to $Nnear$ points is easy:
2 points [-1, $y_{-1}$], [1, $y_1$ ]:
$\qquad$ estimate $y_0 \sim (y_{-1} + y_1) / 2$
points $[x_i, y_i]$ with average $x_i = 0$:
$\qquad y_0 \sim$ average $y_i$
general $[x_i, y_i]$:
$\qquad$ see e.g. Numerical Recipes p. 781: fit a line $a + b x$ and estimate $y_0 \sim a$.
$\qquad$ One can fit quadratics, cubics, sine-cosines ... in the same way.

I understand that you have uniformly-spaced data with a few points missing, is that right ?
How well does linear interpolation work for this case ?
Well, let's try cos $2 \pi f t$ with $f$ = 0.25: 1 0 -1 0 1 0 -1 0 ...
2 neighbors of any point average to 0, terrible.
4 neighbors: average of [ 1 0 (missing -1) 0 1] = 1/2, terrible.
(Try the 4-neighbor filter [-1 3 3 -1] / 4 on this.)

Linear interplation with 4 or 6 or 8 neighbors might work well enough for your data.
I'd suggest starting with a method that you understand thoroughly before diving into splines, sinc-like ... although those can be fun too.

Another, quite different method is Inverse distance weighting . It's easy to implement (see idw-interpolation-with-python on SO), works in 2d 3d and up too, but is afaik hard to analyze theoretically.

(Obviously, NO single interpolation method can possibly fit the zillions of combinations of
[signal, noise, error metric, test function] that occur in reality.
There are more methods in the world, with more knobs, than test functions.
Nonetheless a gallery of methods and test functions might be useful.)

If you work with matlab you can do it by working with timeseries.

time  % is your starting vector of time

data % vector of data you want to resample

data_TS = timeseries(data,time); % define the data as a timeseries

new_time = time(0):dt:time(end); % new vector of time with fixed dt=1/fs

data_res = resample(data_TS,new_time); % data resampled at constant fs


Before you get into doing some exotic processing you could try something simple like this (pseudo code - no interpolation, but that could be added)

TimeStamp[]  //Array of Sensor TimeStamps -NULL terminated – TimeStamp[i] corresponds to Reading[i]

AlgorithmOut   //Delimited file of of readings in fixed sample time (5ms)

SampleTime=TimeStamp //ms virtual sample time, 5ms fixed samples

i = 0 // loop index
While(TimeStamp[i] != NULL)
{
FileWrite (CurrentSavedReading, AlgorithmOut)//write value to file
SampleTime = SampleTime + 5//ms
if(SampleTime > TimeStamp[i])
{
i++