sinc interpolation in MATLAB

UPDATE: The examples given here are meant to give mathematical insight into how sinc interpolation works by using a finite-time APPROXIMATION. Sinc interpolation necessarily involves a circular convolution, which is not a finite computation in the time domain. If you actually need to do sinc interpolation, use the interpft function. It does an FFT, pads the FFT with zeros, and does an IFFT. Consequently, it is VERY FAST. Moreover, using an FFT (or DFT, in general) is the only way to use a finite computation to do a sinc interpolation (or circular convolution in general). Just make sure that the resampled time vector you use has a length that is an integer multiple of your original time vector. Also make sure that it lands on the same points as your original time vector (i.e., it should only add new points between old points).

I've also created a interpftw that does the same job as interpft but allows you to reconstruct samples from different aliasing windows. Whereas interpft pads the end of the FFT with zeros, interpftw pads both the beginning and the end. In other words, the former is a low-pass filter and the latter is a bandpass filter.

If you search Google for sinc interpolation in MATLAB, many pages will reference the sinc_interp example from John Loomis. Unfortunately, I've found few sites that recognize that the function is not meant to do general purpose sinc interpolation. That is, it makes a few assumptions about the sampling rates that may not be evident to the average user.

So, I'm giving some of my students this example:

% 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")

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

    for i=1:length(u)
        y( i ) = sum( x .* sinc( (1/T)*(u(i) - s) ) );
    end

    % Make sure y is same shape as u (row->row, col->col)
    y = reshape(y, size(u));
end

Here's a vectorized (i.e., MUCH FASTER) version:

% Ideally "resamples" x vector from s to u by sinc interpolation
function y = sinc_interp(x,s,u)
    % Interpolates x sampled at "s" uniformly spaced instants
    % Output y is sampled at "u" uniformly spaced instants
    % ("s" for "sampled" and "u" for "upsampled")
    % (consequently, length(x)=length(s))

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

    % The entries of this matrix are each u-s permutation.
    % It will be used to generate the sinc transform that will
    % be convolved below with the input signal to do the
    % interpolation.
    %
    % (recall that u(:) will be a column vector regardless
    % of the row-ness of u. So u(:) is a row, and s(:) is a
    % column)
    sincM = repmat( u(:), 1, length(s) ) ...
           - repmat( s(:)', length(u), 1 );

    % * Sinc is the inverse Fourier transform of the boxcar in
    % the frequency domain that was used to filter out the
    % ambiguous copies of the signal generated from sampling.
    % * That sinc, which is now sampled at length(u) instants,
    % is convolved with the input signal becuse the boxcar was
    % multipled with its Fourier transform.
    % So this multiplication (which is a matrix transformation
    % of the input vector x) is an implementation of a
    % convolution.
    % (reshape is used to ensure y has same shape as upsampled u)
    y = reshape( sinc( sincM/T )*x(:) , size(u) );
end

My function sinc_interp resamples the data in the x vector. The original time vector is given by the s vector and the new time vector is given by the u vector.

Alternatively, you can try this more advanced one that lets you pick which aliasing window you want to use to reconstruct (i.e., it's an ideal bandpass filter rather than just a low-pass filter).

% Ideally "resamples" x vector from s to u by sinc interpolation
function y = sinc_interp(x,s,u,N)
    % Interpolates x sampled sampled at "s" instants
    % Output y is sampled at "u" instants ("u" for "upsampled")
    % Optionally, uses the Nth sampling window where N=0 is DC
    %     (so non-baseband signals have N = 1,2,3,...)

    if nargin < 4
        N = 0;
    end

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

    for i=1:length(u)
        y( i ) = ...
            sum( x .* ...
                ( (N+1)*sinc( ((N+1)/T)*(u(i) - s) ) - ...
                  N*sinc( (N/T)*(u(i) - s) ) ) );
    end
end

Here's a vectorized (i.e., MUCH FASTER) version:

% Ideally "resamples" x vector from s to u by sinc interpolation
function y = sinc_interp(x,s,u,N)
    % Interpolates x sampled sampled at "s" instants
    % Output y is sampled at "u" instants ("u" for "upsampled")
    % Optionally, uses the Nth sampling window where N=0 is DC
    % (so non-baseband signals have N = 1,2,3,...)

    if nargin < 4
        N = 0;
    end

    % 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'(N+1)/T )*x';
    y = x*( (N+1)*sinc( sincM*(N+1)/T ) - N*sinc( sincM*N/T ) );
end