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

Question

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?

Answer 1

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)[0])

    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[0].plot(p[1])

py.show()

Answer 2

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.

Answer 3

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.

Answer 4

(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.)

Answer 5

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]
Reading[]      //Array of Sensor Readings       -NULL terminated

AlgorithmOut   //Delimited file of of readings in fixed sample time (5ms) 
CurrentSavedReading = Reading[0]

SampleTime=TimeStamp[0] //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++
      CurrentSavedReading = Reading[i]
   }
}

Answer 6

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

Answer 7

IMHO Datageist's answer is correct, Jacob's answer is not. An easy way to verify this is that datageist's suggested algorithm is guaranteed to interpolate through the original samples (assuming infinite numerical precision), whereas Jacob's answer does not.

• For the uniform sampling case, the set of sinc functions is orthogonal: if each shifted sinc function is discretized over the input samples, they form an infinite identity matrix. This is because sin(n pi) / (n pi) is zero for all n except n=0.
• However, this property cannot simply be extrapolated to the non-uniform case: a similar set of sinc functions, discretized over the input samples, will yield a nontrivial matrix. Hence, the contributions from surrounding samples will not be zero, and the reconstruction will no longer interpolate through the input samples.