Writing a sinc-interpolating playback object

    christripledot
    Oct 10 2011 | 2:03 pm
    Hello,
    After a lot of reading and head-scratching, I've begun to write a sample playback engine using windowed sinc interpolation. The idea is, it will work a bit like [groove~] (taking a signal to determine playback speed), but it will be click-triggerable and have built-in ADSRs and a variety of looping options.
    The ADSR and loop logic is all done, but the interpolator is giving me grief. Maybe somebody on here is savvy enough to see where I'm going wrong...
    When I pass integer values into my "interp(playhead)" function, everything is hunky-dory. But when I pass fractional values I get hideous ringing and the interpolated values are WAY off.
    First off, I populate a table with the right-hand wing of a sinc function. I then multiply this table with the right-hand wing of a windowing function (Blackman-Harris in my case). I'm pretty damn sure that I'm populating and windowing my table correctly. There are (SINC_POINTS + 1) lobes in it.
    Here is my naive implementation so far:
    #define SINC_RES 512   // resolution of each lobe of the sinc table
#define SINC_POINTS 16 // number of samples to use for interpolation
                       // on either side of the sample in question

float interp(double playhead) {

  long playheadint = floor(playhead);
  double playheadfrac = playhead - playheadint;

  double accum = 0.0;

  // Loop through every sample we're interested in...
  int i;
  for(i = -SINC_POINTS; i < SINC_POINTS + 1; i ++) {

    long thissample = playheadint + i;

    //  HACK: To handle start and end points,
    //  we just duplicate the start and end samples.
    //  POSSIBLE FIX: mirror start/end of input?
    if(thissample < 0) thissample = 0;
    if(thissample >= frames) thissample = frames - 1;

    double moo = samplebuffer[thissample];

    // Get index into sinc table
    // This could be moved outside the loop but it's left here for clarity
    long sincindex = round(playheadfrac * SINC_RES);

    // Offset index relative to centre sample
    sincindex += abs(i) * SINC_RES;

    moo *= sinctable[sincindex];
    accum += moo;
  }

  return (float)accum;
}
    If anybody can help figure out where I'm going wrong, I'd be extremely grateful. Here's what I'm trying to do:
    1. Conceptually, replace each input sample with a sinc function. 2. Take the fractional part of my desired 'sub-sample' and scale this so it can be used as an offset into the first (or main) lobe of the sinc function. 3. Take the relevant value from the sinc function and multiply this by the input sample (as determined by the integer part). 4. Do the same for an arbitrary number of samples on either side of the sample we're interested in, but offset the index into the sinc function accordingly; i.e. if we're looking at one sample away from the main sample, we have to offset the sinc function index by one zero-crossing. 5. Sum the results of the above table lookups.
    Is this the right approach? I think my code above reflects this, but maybe I've just gone code-blind...
    P.S. Maybe I should mention that the eventual goal for this is to create a vari-speed sample player. I don't think this ought to affect the results of my simple function, but perhaps there are other issues to worry about further down the line? Some of the source I've seen for resampling libraries processes the source in chunks, rather than updating on a per-sample basis. However, I suspect this may be for performance reasons. For now I'd just like a slow-but-easy-to-follow piece of code that works...
    If debugging the above is too much of an ask, does anyone have a similar interpolation function they wouldn't mind sharing? One that reports the interpolated value at an arbitrary location in an input buffer?
    Thanks for looking,
    Chris
    • Timothy Place
      Oct 10 2011 | 4:50 pm
      I haven't looked at the details of your code to know exactly what might be going wrong, but I will mention that in Max 6 you can have a sinc-based interpolation in the groove~ object by turning on a new @resample attribute. FWIW, in our case we are windowing the sinc with a Kaiser function.
      There are also some features in buffer~ to make generating the sinc and windowing it super-easy.
      We generally aren't sharing Max 6 details just yet, but will be sharing more soon. I thought I'd give you a heads-up given that it may impact your (really interesting and cool!) project.
      Cheers, Tim
    • christripledot
      Oct 10 2011 | 5:44 pm
      Thanks, Tim! That's wicked news.
      I'm still hoping someone with more knowledge than me might chime in though, 'cos I'd like to trigger this at audio rate, and I'm interested in more advanced looping than groove~ currently offers (fwd, reverse, pingpong). I'll share once I've got to the bottom of this! :)
    • christripledot
      Oct 11 2011 | 1:35 pm
      A kind soul pointed out my error. Because I'm only using half a sinc function to save memory, I'm mirroring the index into the table. But I was only mirroring the integer part of the index.
      This is the offending line:
      // Offset index relative to centre sample
sincindex += abs(i) * SINC_RES;
      It should be:
      sincindex = abs(sincindex - (i * SINC_RES));
      I can't test at work, but I think that's cracked it. I'll share the external once it's tidied up a bit.