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

Rate this code snippet:

Rating: 4 | Votes: 3

posted by Markus Nentwig
Markus received his Dipl. Ing. degree in electrical engineering / communications in 1999. Work interests include RF transceiver system design, implementation, modeling and verification. He works as senior architect for Renesas Mobile Europe in Finland.

mnentwig wrote:

11/4/2010

Thanks for the comment!

But keep in mind that the error is frequency dependent (you won't see any interpolation error for a DC signal) and delay dependent (again, no error for zero delay).

If you're interested in resampling, I strongly recommend to have a look at linear filters instead of the spline "trickstery". It works very well for the human ear, but I'd be careful in communications applications or the like. For a "proper" interpolator I'd decide on a maximum bandwidth and design an equiripple filter (Parks McClellan) or the like.

If you need a reference delay with practically zero error, try this:
http://www.elisanet.fi/~d635415/webroot/MatlabOctaveBlocks/mn_FFT_delay.m
The only catch is that it treats the signal as periodic. Which is the same as having the "repeat" button pushed on the CD player if I process a whole song, almost a philosophical question...

Sign in

Search code

See Also

DSP Code Sharing > Audio interpolation (Catmull-Rom/Cubic Hermite Spline)

Audio interpolation (Catmull-Rom/Cubic Hermite Spline)

Audio interpolation (Catmull-Rom/Cubic Hermite Spline)

oliviert wrote:

mnentwig wrote:

gilgamash wrote: