Sinc Interpolation problem

[DummyDSP]

Hi everyone,
I have been trying to understand sinc interpolation for a while now and have coded several versions of it.
The problem is, as you might quess, it doesn't work.

I have tested the interpolation code with pure sine wave changing the frequency from 50hz to something bigger like 100hz.
The latest code I have been using is from here viewtopic.php?t=330731#

Code: Select all

#define SINC_WIDTH 32
#define SAMPLES_PER_ZERO_CROSSING 11500

float sincInterp(float playhead, float* input, int sourceLength){
        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_WIDTH + 1; i < SINC_WIDTH -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 >= sourceLength) thissample = sourceLength - 1;
            
            double moo = input[thissample];
            
            float sincindex = fabs(round((i+playheadfrac) * SAMPLES_PER_ZERO_CROSSING));
            moo *= sinc(sincindex);
            accum += moo;
        }
        
        return (float)accum;
}

sinc generation code (copied from somewhere I don't remember)

Code: Select all

void make_sinc(){
    int i;
    double temp,win_freq,win;
    win_freq = PI / SINC_WIDTH / SAMPLES_PER_ZERO_CROSSING;
    sinc_table[0] = 1.0;
    for (i=1;i<SINC_WIDTH * SAMPLES_PER_ZERO_CROSSING;i++)   {
        temp = (double) i * PI / SAMPLES_PER_ZERO_CROSSING;
        sinc_table[i] = sin(temp) / temp;
        win = 0.5 + 0.5 * cos(win_freq * i);
        sinc_table[i] *= win;
    }
}

sinc is a method where the value is got from a sinc table

Code: Select all

...
temp = fabs(x);
low = temp;          /* these are interpolation steps */
delta = temp - low;  /* and can be ommited if desired */
return linear_interp(sinc_table[low],sinc_table[low + 1],delta);
...

Using 50hz sine wave with rate = 1.1 the output waveform is like this.



And frequency spectrum like this.



We can clearly see that the interpolation is not perfect. In fact linear interpolation would be much better than this

It doesn't help If I change SINC_WIDTH or SAMPLES_PER_ZERO_CROSSING to bigger values. The quality stays the same.
I'm making some mistake in somewehere but I don't know where.

I don't wan't opinions about I should use some other interpolation methods. I just want to learn why my code doesn't work and after then I can move on and perhaps use some other method for interpolation

I would really appreciate your help

[kryptonaut]

My first guess would be something wrong with your sinc table, maybe in the treatment of SAMPLES_PER_ZERO_CROSSING.

I would try feeding a delta impulse into your algorithm, with a fractional offset to delay/advance it by a small amount, and check that you see something sync-like.

[DummyDSP]

My first guess would be something wrong with your sinc table, maybe in the treatment of SAMPLES_PER_ZERO_CROSSING.

I would try feeding a delta impulse into your algorithm, with a fractional offset to delay/advance it by a small amount, and check that you see something sync-like.

Thanks I updated the OP with sinc generation code I use.

[noizebox]

Got confused for a bit by the one-sided sinc generation. But looking at the expression for sincindex quickly, that looks a bit peculiar. I would expect sincindex to be roughly symmetrical when looking a the mix and max calculated for every sample, but with i+playheadfrac in there it def becomes non-symmetrical.

Also, kryptonauts suggestion is a good one.

[kryptonaut]

Ok, I think the problem is where you calculate sincindex, I reckon you should be using (i-playheadfrac) not (i+playheadfrac)

If you work through what happens when playhead is, say, 99.9 then you ought to see a weighting of nearly 1 given to the 100th sample, by virtue of sincindex being close to zero, but it doesn't work out that way. For playhead=100.1 you should see the same weighting for the 100th sample, and in this case you do. But by changing the sign of playheadfrac as suggested, I think you'll get the correct value for sincindex both times - both will be close to zero, and the nearer playhead is to 100 (either from above or below) the closer sincindex (for the 100th sample) will be to zero.

[DummyDSP]

Ok, I think the problem is where you calculate sincindex, I reckon you should be using (i-playheadfrac) not (i+playheadfrac)

If you work through what happens when playhead is, say, 99.9 then you ought to see a weighting of nearly 1 given to the 100th sample, by virtue of sincindex being close to zero, but it doesn't work out that way. For playhead=100.1 you should see the same weighting for the 100th sample, and in this case you do. But by changing the sign of playheadfrac as suggested, I think you'll get the correct value for sincindex both times - both will be close to zero, and the nearer playhead is to 100 (either from above or below) the closer sincindex (for the 100th sample) will be to zero.

Thanks man. We are getting there...

You mean like this?

Code: Select all

float sincindex = (abs(i)+playheadfrac) * SAMPLES_PER_ZERO_CROSSING;
if(i>0){
  sincindex = (abs(i-1)+(1-playheadfrac)) * SAMPLES_PER_ZERO_CROSSING;
}

This code improves things but not completely.

Screen Shot 2017-10-04 at 20.20.05.png

This is how I 'see' this whole situation
The red dot is the sample to be interpolated. black dots are initial sample points

sinc_interp.png

[kryptonaut]

No, I meant simply:

Code: Select all

float sincindex = fabs(round((i-playheadfrac) * SAMPLES_PER_ZERO_CROSSING));

Try working through it with playhead equal to, say, 99.9 and then 100.1, and see the sincindex calculated for samples 99, 100 and 101 - the sets of values will be nicely symmetric.

[DummyDSP]

No, I meant simply:

Code: Select all

float sincindex = fabs(round((i-playheadfrac) * SAMPLES_PER_ZERO_CROSSING));

Try working through it with playhead equal to, say, 99.9 and then 100.1, and see the sincindex calculated for samples 99, 100 and 101 - the sets of values will be nicely symmetric.

Thanks.
Now I get it

For others who are interested:
If we take look my drawing in previous post and assume playheadfrac = 0.1
Going LEFT from red dot
if i = 0 then sincindex should be 0.1
if i = -1 then sincindex should be 1.1
if i = -2 then sincindex should be 2.1
...

Going RIGHT from red dot
if i = 1 then sincindex should be 0.9
if i = 2 then sincindex should be 1.9
if i = 3 then sincindex should be 2.9
....

Magically this code produces the same result

Code: Select all

float sincindex = fabs(round((i-playheadfrac) * SAMPLES_PER_ZERO_CROSSING));

Your code makes things better but not even close to perfect.

The first image is sinc interpolator with SINC_WIDTH 32 (larger SINC_WIDTH didn't improve the situation)
The second image is simple linear interpolation

Screen Shot 2017-10-05 at 21.49.09.png

Screen Shot 2017-10-05 at 21.55.25.png

btw. the weird freq spectrums in first post were due generating the wav file in separate program and then importing it to Reaper. The process introduced additional error.
Now I use the code directly in my VST plugin.

[noizebox]

Why do you need to round sincindex before taking the absolute value of it? You're still doing a linear interpolation when looking up the value from the sinc table, right?

[kryptonaut]

Also, the limits of the 'for' loop in sincInterp() look a bit suspect.

Before trying to alter the sampling frequency, I would try just shifting everything by a fractional offset and see if you get a clean-looking result. Try plotting the output rather than the spectrum, and try plotting the difference between the output and the calculated expected output.

Also, if you try adding all the weights that are used when calculating an output sample, I suspect they will not sum to exactly 1 (because of the windowing function) - the effect will be that using this function to interpolate a DC signal will produce a slightly different DC signal. Clearly a linear interpolation of DC will be exactly right. So maybe it's understandable that at low frequencies a linear interpolation may outperform this function, with some kind of crossover happening as the frequency increases and linear interpolation becomes the worse option?

[DummyDSP]

Why do you need to round sincindex before taking the absolute value of it? You're still doing a linear interpolation when looking up the value from the sinc table, right?

There is no need to do that. It was just accidentally left there when copyin code from several places.
Removing it doesn't really improve things because SAMPLES_PER_ZERO_CROSSING is already quite large.

[DummyDSP]

Also, the limits of the 'for' loop in sincInterp() look a bit suspect.

Before trying to alter the sampling frequency, I would try just shifting everything by a fractional offset and see if you get a clean-looking result. Try plotting the output rather than the spectrum, and try plotting the difference between the output and the calculated expected output.

Also, if you try adding all the weights that are used when calculating an output sample, I suspect they will not sum to exactly 1 (because of the windowing function) - the effect will be that using this function to interpolate a DC signal will produce a slightly different DC signal. Clearly a linear interpolation of DC will be exactly right. So maybe it's understandable that at low frequencies a linear interpolation may outperform this function, with some kind of crossover happening as the frequency increases and linear interpolation becomes the worse option?

I changed the limits to

Code: Select all

for(int i = -SINC_WIDTH + 2 ; i < SINC_WIDTH; i ++) {

it didn't improve the quality although I think this is more correct

Now I have tried every trick I can imagine.
Maybe I have to admit that sinc interpolation is not as good as I thought. At least in this simple form of implementation. To increase the quality I have to do something else which I don't know at the moment.

Now the quality is 'good' if we don't mind the noise below the actual waveform frequency spectrum as you can see in these images. Maybe this is how sinc interpolation behaves. Who knows. Please tell us Who

Played C6 in Reaper
#define SINC_WIDTH 128
#define SAMPLES_PER_ZERO_CROSSING 10000

Screen Shot 2017-10-06 at 21.35.02.png

Screen Shot 2017-10-06 at 21.34.40.png