Sinc interpolation problems

[christripledot]

Hello all,

After a lot of reading and head-scratching, I'm trying to write a sample playback engine using windowed sinc interpolation. 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:

Code: Select all

#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 = input[thissample];

    // Get index into sinc table
    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.

Thanks for looking,

Chris

[jupiter8]

Just a thought: shouldn't the sinc table be heavily oversampled for this to work ?

[christripledot]

I'm not sure what you mean. The table has 512 entries per 'lobe', and its frequency is Nyquist - meaning that there is a lobe, or zero-crossing, for every sample in my input signal.

The resolution could be increased more but it makes no difference to the erroneous output. Once I get it working I'll incorporate linear interpolation between each table entry.

The way I understand it, this is an oversampled table already, but I freely admit I feel pretty out of my depth here.

[jupiter8]

Ah ok i just quickly glanced thru it and thought you only had 16 sinc points in total.

[christripledot]

No worries! Thanks for your suggestion all the same.

Wow, 122 views so far? Lots of interest, but not much input... Either I'm embarassingly stupid, or this is a tough topic for a lot of people! I'm sorry to just dive in and say "debug my code," it's probably seen as poor form. Anyway...

I just can't see what I'm doing wrong here. Perhaps I just haven't understood the process correctly? 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...

Cheers!
- Chris

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?

[Borogove]

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?

Olli Niemitalo's "Elephant Paper" describes a lot of resampling interpolation algorithms with different quality/performance tradeoffs. If you haven't already, you should implement straight linear interpolation to make sure that the rest of your resampling engine works properly.

[christripledot]

That's a damn good paper, thanks! It's got a few hairy equations in it, but I'll try and puzzle something out.

In the meantime, I have linear and cubic interpolation working just fine. I'm sure there's just something silly in my sinc algo that's screwing things up. I can't see it though!

[sonigen]

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.

Because you are using the fractional part to index into the sinc, and in fact index in the full (or half?) length of it.

If you have 512 tap FIR, and are using 32 samples to interpolate, then you should using the fractional position to index in by 16 samples.

position = position to evaluate
idx = floor(position) = integer position
frac = position-idx = fractional part of position
firlen = length of FIR filter kernel, the sinc
sampletaps = number of source samples being used
ostaps = firlen / sampletaps = number of FIR taps per source sample
firoffset = floor(frac * ostaps) = offset into FIR

x = 0.0;
for (int i = 0; i < sampletaps; i++)
{
x += sample * sinc[i*sampletaps + firoffset];
}

Something like that.

You can simplify things by storing the FIR in polyphase form. It'll simplify the indexing and make it better regards the CPU cache.

[christripledot]

Chris, you have just spared me a lot of pain! Thank you so much!!!

[christripledot]

x += sample * sinc[i*sampletaps + firoffset];

Did you mean:

x += sample * sinc[i*ostaps + firoffset];

?

[christripledot]

Because you are using the fractional part to index into the sinc, and in fact index in the full (or half?) length of it.

Sorry, I've done a bit of re-reading and I'm confused. If you are saying that I'm mistakenly using the fractional part to index into the whole sinc table, I'm not. Check it out:

Code: Select all

// Get index into sinc table
long sincindex = round(playheadfrac * SINC_RES); 

SINC_RES is the number of table entries per lobe (per input sample), not the length of the whole sinc.

Now I'm even more confused :S

...but still ever so grateful for your input!

[christripledot]

Ahhhhhhh. I get it now. I stumbled upon this page: http://www.nicholson.com/rhn/dsp.html and reverse-engineered the BASIC code.

Thanks for your time and patience

[sonigen]

Ah sorry i thought SINC_RES was the length of your table.

So the length of the table is SINC_RES * SINC_POINTS

That makes more sense!

In which case i think you want this...

sincindex = abs(round((i+playaheadfrac) * SINC_RES))

reason being that because you're only using half the SINC kernel, and you're mirroring the lookup, you need to mirror both the fractional position and the integer part. You were only mirroring the integer part and not the fractional, for example...

if integer pos = -4, and fractional offset is 0.3, then the actual offset is

-4 + 0.3 == -3.7

But your doing

abs(4) + 0.3 = 4.3

When it should be abs(-3.7) = 3.7

[christripledot]

Yes, I sussed that out!

I'd tried to implement that fix before, but somehow buggered it up and assumed I must have been missing something else. All sorted now.

P.S. Thanks for the polyphase tip; it makes parallelism a bit easier too!

Rice and peas