-1

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

user1195202
  • 1,423
  • 1
  • 16
  • 20
  • 1
    You might door worse than starting with a half-band filter (https://en.wikipedia.org/wiki/Half-band_filter). Building IIR or FIR filters is easy from an algorithmic point of view - they're just weighted sums. Finding the weights on the other hands rather harder - many filter topologies are analytically insoluble, so are solved by numerical methods. There are loads of Matab based filter solvers out there. – marko Mar 24 '19 at 23:32
  • Thanks for the answer. I tried to find a code for this half-band filter. But looks like it's not working good. Also I tried to reproduce wave with sum of sinc impulses, but I still get some harmonics. Do you have some example in c/c++? – user1195202 Mar 31 '19 at 10:00
  • I think you need to improve the question. It's not at all clear what problem you have. – Alexey Frunze Mar 31 '19 at 10:33
  • Perhaps you would care to add that to the question? – marko Mar 31 '19 at 22:06
  • Note that we prefer a technical style of writing here. We gently discourage greetings, hope-you-can-helps, thanks, advance thanks, notes of appreciation, regards, kind regards, signatures, please-can-you-helps, chatty material and abbreviated txtspk, pleading, how long you've been stuck, voting advice, meta commentary, etc. Just explain your problem, and show what you've tried, what you expected, and what actually happened. – double-beep May 27 '19 at 11:01

1 Answers1

0

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.

user1195202
  • 1,423
  • 1
  • 16
  • 20