Add anti-aliasing/bandlimit for looped wav sample (NOT Fourier transform)

Question

How to build antialiasing interpolation using c++ code? I have a simple 4096 or 1024 buffer. Of course when I play this at high frequencies I get aliasing issues. to avoid this, the signal must be limited by the bandwidth at high frequencies. Roughly speaking, the 'sawtooth' wave at high frequencies should looks like a regular sine. And that is what I want to get so that my sound didn't sound dirty like you moving knobs in your old FM/AM radio in your car.

I know how to build bandlimited square,triangle,sawtoth with Fourier transform. So my question is only about wavetable

Answer 1

Found solution in the AudioKit sources. One buffer will be split into 10 buffers/octaves. So when you play a sound, you don't play your original wave, but play a sample that was prepared for a specific octave.

Import to your project WaveStack.hpp

namespace AudioKitCore
{

    // WaveStack represents a series of progressively lower-resolution sampled versions of a
    // waveform. Client code supplies the initial waveform, at a resolution of 1024 samples,
    // equivalent to 43.6 Hz at 44.1K samples/sec (about 23.44 cents below F1, midi note 29),
    // and then calls initStack() to create the filtered higher-octave versions.
    // This provides a basis for anti-aliased oscillators; see class WaveStackOscillator.

    struct WaveStack
    {
        // Highest-resolution rep uses 2^maxBits samples
        static constexpr int maxBits = 10;  // 1024

        // maxBits also defines the number of octave levels; highest level has just 2 samples
        float *pData[maxBits];

        WaveStack();
        ~WaveStack();

        // Fill pWaveData with 1024 samples, then call this
        void initStack(float *pWaveData, int maxHarmonic=512);

        void init();
        void deinit();

        float interp(int octave, float phase);
    };

}

WaveStack.cpp

#include "WaveStack.hpp"
#include "kiss_fftr.h"

namespace AudioKitCore
{

    WaveStack::WaveStack()
    {
        int length = 1 << maxBits;                  // length of level-0 data
        pData[0] = new float[2 * length];           // 2x is enough for all levels
        for (int i=1; i<maxBits; i++)
        {
            pData[i] = pData[i - 1] + length;
            length >>= 1;
        }
    }

    WaveStack::~WaveStack()
    {
        delete[] pData[0];
    }

    void WaveStack::initStack(float *pWaveData, int maxHarmonic)
    {
        // setup
        int fftLength = 1 << maxBits;
        float *buf = new float[fftLength];
        kiss_fftr_cfg fwd = kiss_fftr_alloc(fftLength, 0, 0, 0);
        kiss_fftr_cfg inv = kiss_fftr_alloc(fftLength, 1, 0, 0);

        // copy supplied wave data for octave 0
        for (int i=0; i < fftLength; i++) pData[0][i] = pWaveData[i];

        // perform initial forward FFT to get spectrum
        kiss_fft_cpx spectrum[fftLength / 2 + 1];
        kiss_fftr(fwd, pData[0], spectrum);

        float scaleFactor = 1.0f / (fftLength / 2);

        for (int octave = (maxHarmonic==512) ? 1 : 0; octave < maxBits; octave++)
        {
            // zero all harmonic coefficients above new Nyquist limit
            int maxHarm = 1 << (maxBits - octave - 1);
            if (maxHarm > maxHarmonic) maxHarm = maxHarmonic;
            for (int h=maxHarm; h <= fftLength/2; h++)
            {
                spectrum[h].r = 0.0f;
                spectrum[h].i = 0.0f;
            }

            // perform inverse FFT to get filtered waveform
            kiss_fftri(inv, spectrum, buf);

            // resample filtered waveform
            int skip = 1 << octave;
            float *pOut = pData[octave];
            for (int i=0; i < fftLength; i += skip) *pOut++ = scaleFactor * buf[i];
        }

        // teardown
        kiss_fftr_free(inv);
        kiss_fftr_free(fwd);
        delete[] buf;
    }

    void WaveStack::init()
    {
    }

    void WaveStack::deinit()
    {
    }

    float WaveStack::interp(int octave, float phase)
    {
        while (phase < 0) phase += 1.0;
        while (phase >= 1.0) phase -= 1.0f;

        int nTableSize = 1 << (maxBits - octave);
        float readIndex = phase * nTableSize;
        int ri = int(readIndex);
        float f = readIndex - ri;
        int rj = ri + 1; if (rj >= nTableSize) rj -= nTableSize;

        float *pWaveTable = pData[octave];
        float si = pWaveTable[ri];
        float sj = pWaveTable[rj];
        return (float)((1.0 - f) * si + f * sj);
    }

}

Then use it in this way:

//wave and outputWave should be float[1024];
void getSample(int octave, float* wave, float* outputWave){

    uint_fast32_t impulseCount = 1024;

    if (octave == 0){
        impulseCount = 737;
    }else if (octave == 1){
        impulseCount = 369;
    }
    else if (octave == 2){
        impulseCount = 185;
    }
    else if (octave == 3){
        impulseCount = 93;
    }
    else if (octave == 4){
        impulseCount = 47;
    }
    else if (octave == 5){
        impulseCount = 24;
    }
    else if (octave == 6){
        impulseCount = 12;
    }
    else if (octave == 7){
        impulseCount = 6;
    }
    else if (octave == 8){
        impulseCount = 3;
    }
    else if (octave == 9){
        impulseCount = 2;
    }

    //Get sample for octave
    stack->initStack(wave, impulseCount);

    for (int i = 0; i < 1024;i++){
        float phase = (1.0/float(1024))*i;

        //get interpolated wave and apply volume compensation
        outputWave[i] = stack->interp(0, phase) / 2.0;
    }

}

Then... when 10 buffers is ready. You can use them when playing a sound. Using this code you can get index of buffer/octave depending to your frequency

uint_fast8_t getBufferIndex(const float& frequency){

    if (frequency >= 0 && frequency < 40){
        return 0;
    }
    else if (frequency >= 40 && frequency < 80){
        return 1;
    }else if (frequency >= 80 && frequency < 160){
        return 2;
    }else if (frequency >= 160 && frequency < 320){
        return 3;
    }else if (frequency >= 320 && frequency < 640){
        return 4;
    }else if (frequency >= 640 && frequency < 1280){
        return 5;
    }else if (frequency >= 1280 && frequency < 2560){
        return 6;
    }else if (frequency >= 2560 && frequency < 5120){
        return 7;
    }else if (frequency >= 5120 && frequency < 10240){
        return 8;
    }else if (frequency >= 10240){
        return 9;
    }
    return 0;
}

So if I know that my note frequency 440hz. Then for this note I'm getting wave in this way:

float notInterpolatedSound[1024];
float interpolatedSound[1024];

uint_fast8_t octaveIndex = getBufferIndex(440.0);
getSample(octaveIndex, notInterpolatedSound, interpolatedSound);

//tada! 

ps. the code above is a low pass filter. I also tried to do sinc interpolation. But sinc worked for me very expensive and not exactly. Although maybe I did it wrong.