0
$\begingroup$

I'm currently trying to implement a windowed sinc interpolation. I've already written some code for that, but it only seems to work for cases where phaseInc <= 1.0 (phaseInc = outSampleRate / inSampleRate). By "not working" I mean that something is causing a lot of distortion or aliasing. I hope that the comments in the code do explain what I'm trying to do:

// SINC_FILT_THRESH = 0.8, SINC_WINDOW_SIZE = 16

float sincStep = phaseInc > SINC_FILT_THRESH ? SINC_FILT_THRESH / phaseInc : 1.00f;
/* ^^ this should in theory detect the transition between phaseInc <= 1.0 and phaseInc > 1.0
 * but with a little headroom for the transition band. "sincStep" then shall be used later
 * to make the sinc function wider to narrow down the rectangle in frequency domain
 * The code below seems to work as long as sincStep = 1.0 */

int i = 0;
do {
    float sampleSum = 0.0f;
    float kernelSum = 0.0f;
    for (int wi = -SINC_WINDOW_SIZE + 1; wi <= SINC_WINDOW_SIZE; wi++) {
        float sincIndex = float(wi) * sincStep - phase;
        float windowIndex = (float(wi) - phase) / SINC_WINDOW_SIZE;

        float s = sinc(sincIndex);
        /* ^^ this is a normalized sinc function with zeroes at whole numbers */

        float w = windowFunc(windowIndex);
        /* ^^ windowFunc = 0.5 + (0.5 * cos(pi * t))
         * I know there is stuff like blackman window, but I'd rather like to start with
         * something simple */

        float kernel = s * w;
        sampleSum += kernel * fetchBuffer[i + wi + SINC_WINDOW_SIZE - 1];
        kernelSum += kernel;
    }
    phase += phaseInc;
    int istep = static_cast<int>(phase);
    phase -= static_cast<float>(istep);
    i += istep;
    /* ^^ This simply advanced phase and i by the appropriate amount.
     * I used it with a simple linear interpolator and didn't have problems with that */

    *outData++ = sampleSum / kernelSum;
    /* ^^ Store the result in the output buffer with unity gain */
} while (--numSamples > 0);

So, as I've mentioned it only seems to work right in cases where sincStep = 1.0 and I don't know why. My intention is to widen the sinc function to lower the cutoff. Doesn't seem to work quite as I want.

What I do know is that there are no memory access violation problems (according to fsanitize) and my fetchBuffer always has enough data (so no out of bounds thing).

I'm sorry if I'm missing something quote obvious but I really don't want to just copy someone elses code. I'd like to actually understand it. Most examples however just seem to focus on the equations and not on the actual thing with all the possibly faulty boiler plate code.

Thanks in advance!

EDIT: I want to use that interpolation for resampling different waves at different pitch for music. I thought I could add this for clarification for my problem: spectogram

$\endgroup$
  • $\begingroup$ Hi! Do you have access to Matlab / Octave (free) ? $\endgroup$ – Fat32 Jan 11 at 20:56
  • $\begingroup$ Well, I think we have Matlab at university, so yes. $\endgroup$ – theIpatix Jan 11 at 21:02
  • $\begingroup$ So, being quite lazy in debugging (other people's) C/C++ code, I would kindly ask whether you were able to (successfully) implement your algorithm in Matlab, before moving it into C++ ? $\endgroup$ – Fat32 Jan 11 at 21:11
  • $\begingroup$ I guess I could do that, but frankly enough I have never used Matlab before so I'd guess it's probably gonna be more buggy than this. $\endgroup$ – theIpatix Jan 11 at 21:21
  • 1
    $\begingroup$ Still wanted to thank you for the advice, even though I solved the issue. Perhaps Learning some Matlab for these kinds of things seems like a decent idea. $\endgroup$ – theIpatix Jan 11 at 23:36
0
$\begingroup$

So, sleeping a night before posting the question on stackexchange won't help but trying out random things AFTER posting the question results in magic results.

It seems like the following line was faulty:

float sincIndex = float(wi) * sincStep - phase;

The solution was to change it to this:

float sincIndex = (float(wi) - phase) * sincStep;
$\endgroup$

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.