Perfect sinc interpolation in Matlab and Python

sinc_interp.m

% From http://phaseportrait.blogspot.com/2008/06/sinc-interpolation-in-matlab.html

% Ideally "resamples" x vector from s to u by sinc interpolation
function y = sinc_interp(x,s,u)
    % Interpolates x sampled sampled at "s" instants
    % Output y is sampled at "u" instants ("u" for "upsampled")
    % (EXPECTS x, s, and u to be ROW VECTORS!!)

    % Find the period of the undersampled signal
    T = s(2)-s(1);

    % When generating this matrix, remember that "s" and "u" are
    % passed as ROW vectors and "y" is expected to also be a ROW
    % vector. If everything were column vectors, we'd do.
    %
    % sincM = repmat( u, 1, length(s) ) - repmat( s', length(u), 1 );
    %
    % So that the matrix would be longer than it is wide.
    % Here, we generate the transpose of that matrix.
    sincM = repmat( u, length(s), 1 ) - repmat( s', 1, length(u) );

    % Equivalent to column vector math:
    % y = sinc( sincM'/T )*x';
    y = x*sinc( sincM/T );
end

sinc_interp.py

def sinc_interp(x, s, u):
    """
    Interpolates x, sampled at "s" instants
    Output y is sampled at "u" instants ("u" for "upsampled")
    
    from Matlab:
    http://phaseportrait.blogspot.com/2008/06/sinc-interpolation-in-matlab.html        
    """
    
    if len(x) != len(s):
        raise Exception, 'x and s must be the same length'
    
    # Find the period    
    T = s[1] - s[0]
    
    sincM = tile(u, (len(s), 1)) - tile(s[:, newaxis], (1, len(u)))
    y = dot(x, sinc(sincM/T))
    return y

What if T is not fixed. s(n)-s(n-1) is chaning all the time?

newaxis is not defined

@NadavB. They're all numpy functions (including newaxis). Presumably the author is simply omitting a from numpy import *.

It can be rewritten as follows:

import numpy as np

def sinc_interp(x, s, u):
    if len(x) != len(s):
        raise ValueError('x and s must be the same length')
    
    # Find the period    
    T = s[1] - s[0]
    
    sincM = np.tile(u, (len(s), 1)) - np.tile(s[:, np.newaxis], (1, len(u)))
    y = np.dot(x, np.sinc(sincM/T))
    return y

Great, perfect for variable sample rate conversion - works like a charm! Thanks!

I tried this but when i applied for out of boundary sample points this method is not working.

Does anybody has experience with extending this method with a windowed sinc function (Lanczos interpolation)?

sincM = tile(u, (len(s), 1)) - tile(s[:, newaxis], (1, len(u))) can be rewritten as just sincM = u - s[:, None].

Note however that this implementation scales as u.size * s.size but a much faster one could be written using the convolution theorem and the FFT.

sincM = tile(u, (len(s), 1)) - tile(s[:, newaxis], (1, len(u))) can be rewritten as just sincM = u - s[:, None].

Note however that this implementation scales as u.size * s.size but a much faster one could be written using the convolution theorem and the FFT.

How could a faster one be written using convolution theorem and FFT?

This is for python 2 or 3?
What types should be x, s, u?
I tried with List[float]:

x = [1]*25 + [0]*75
s = np.arange(0.0, 1.0, 0.01).tolist()
u = np.arange(0.0, 1.0, 0.0001).tolist()

y = sinc_interp(x, s, u)

just to test this, to understand,
but got error:
TypeError: list indices must be integers or slices, not tuple

Any simple usage example?

ETA: now I see that this should be np.array.
It would be nice to have it hinted somehow.

thanks for the code! for simple speedup, based on this Gist see here - also with a FFT based implementation.