Skip to content

Instantly share code, notes, and snippets.

@endolith
Created October 19, 2011 01:05
Show Gist options
  • Save endolith/1297227 to your computer and use it in GitHub Desktop.
Save endolith/1297227 to your computer and use it in GitHub Desktop.
Perfect sinc interpolation in Matlab and Python
% 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
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
@hahazhar
Copy link

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?

@bicyclesonthemoon
Copy link

bicyclesonthemoon commented Oct 4, 2022

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.

@fschwar4
Copy link

fschwar4 commented Jun 24, 2023

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment