Audio interpolation (Catmull-Rom/Cubic Hermite Spline)

Matlab example for Cubic Hermite Spline audio interpolation for sample-rate conversion or variable-speed playback. It is computationally efficient and audio quality is good.

% *****************************************************************
% Audio interpolation using Cubic Hermite Spline (Catmull-Rom)
%
% The audio file samples describe a continuous-waveform. 
% The task is to interpolate between samples, for example in sample rate
% conversion or variable-speed playback.
% The cubic spline approach is one particular trade-off between 
% quality and efficiency.
% 
% For an efficient C-implementation, refer to fluid_rvoice_dsp.c
% in the fluidsynth project on sourceforge.
% Based on work by Olli Niemitalo:
% http://yehar.com/blog/wp-content/uploads/2009/08/deip.pdf page 11

% ******************************************************************
% build coefficient matrix
% ******************************************************************
% ideally, the interpolation coefficients would be recalculated for any possible
% value of ptrFrac below.
% example input file: http://www.elisanet.fi/mnentwig/webroot/common/in.wav
% Note: This example is written for clarity, not speed.
nPoly=128;
x=((0:nPoly-1)/nPoly)'; % column vector

c1=(x .* (-0.5 + x .* (1 - 0.5 .* x)));
c2=(1.0 + x .* x .* (1.5 .* x - 2.5));
c3=(x .* (0.5 + x .* (2.0 - 1.5 .* x)));
c4=(0.5 .* x .* x .* (x - 1.0));

% ******************************************************************
% load input data
% ******************************************************************
[y, fs, bits]=wavread('in.wav');
[nSamples, nChan]=size(y);
inputdata=transpose(y(:, 1)); % first channel if stereo

% prepare output data storage
outIndex=1; outputData=[]; nAlloc=0;

% ******************************************************************
% "ptr" indicates a position in the input data. 
% Loop through the whole range until end-of-data.
% ******************************************************************
ptr=2;
while ptr < nSamples-2

    % "rate" determines the number of input samples per output sample
    % rate=1.1; % faster constant-speed playback
    % variable playback speed (arbitrary example)
    rate=1+0.3*sin(outIndex/40000);;

    % ******************************************************************
    % Interpolate samples
    % ******************************************************************
    % ptr is the index into the input data, needed for the next output sample.
    % Usually, it lies betwen points, for example ptr=1234.56

    % integer part (1234)
    ptrBase=floor(ptr);
    % fractional part [0..1[, for example 0.56
    ptrFrac=ptr-ptrBase;

    % look up the filter for the fractional part
    fIndex=floor(ptrFrac*nPoly)+1;

    % input samples
    s1=inputdata(ptrBase-1);
    s2=inputdata(ptrBase);
    s3=inputdata(ptrBase+1);
    s4=inputdata(ptrBase+2);

    % performance optimization. Reallocate some extra space.
    if outIndex>=nAlloc
        nAlloc=nAlloc+2000;
        outputData(nAlloc)=0;
    end

    % ******************************************************************
    % Calculate output sample.
    % ******************************************************************
    outputData(outIndex)=s1*c1(fIndex)+s2*c2(fIndex)+s3*c3(fIndex)+s4*c4(fIndex);

    % that's it. Next output sample
    outIndex=outIndex+1;

    % advance the index in the source data
    ptr=ptr+rate;
end

% performance optimization part II, chop extra memory
outputData=outputData(1:outIndex-1)';

% ******************************************************************
% Write output to file
% ******************************************************************
% wavwrite('out.wav', outputData, fs, bits); % OCTAVE version
wavwrite(outputData, fs, bits, 'out.wav'); % MATLAB version