Tutorial: BLEPs (using PolyBLEPs but extensible)

DSP, Plugin and Host development discussion.
RELATED
PRODUCTS

Post

Everyone loves BLEPs, and I've been meaning to write an example / tutorial on this stuff for a while now.. but the problem is that you need three different things: the actual oscillator code, the BLEP generation or data, and buffer management.

But I figured, if we use 2-sample PolyBLEPs then buffering becomes a trivial one sample delay thing (that we still do here, but no ring-buffers or anything), and the BLEP data is no longer an issue. The exact same concepts works, we just need less "dirty details" that make writing tutorials more painful. I'll probably give up and provide the "dirty details" in a later post anyway, but .. yeah ..

This should be beginner friendly, almost "your first oscillator" level.

However: this is still a "real" BLEP oscillator: you just replace the "one sample delay" that's held in two variables, with a proper ring-buffer, and you can use longer BLEPs without any additional trouble.. in fact it's a simplified and cleaned up version of code that uses polyphase FIR BLEPs.

For this example.. I support constant frequency, and constant mix between saw/pulse, but really these could come from buffer just as well. The PWM does come from a buffer, since I added a fix for more robust PWM, so it kinda made sense. The "broken" version is left as an #ifdef just so you can observe how it appears to works about 99% of time, but then occasionally with certain types of modulation there can be glitches.

Anyway, first we need the BLEPs, and here we use two-sample polynomial filter.

These appear to be what are most commonly suggested:

Code: Select all

// The first sample BLEP
static inline float poly3blep0(float t)
{
    // these are just sanity checks
    // correct code doesn't need them
    if(t < 0) return 0;
    if(t > 1) return 1;

    float t2 = t*t;
    return t * t2 - 0.5f * t2 * t2;
}

// And second sample as wrapper, optimize if you want.
static inline float poly3blep1(float t)
{
    return -poly3blep0(1-t);
}
[/size]

How I derived it (if you actually care): Build C1 continuous impulse from two cubic "smooth steps" which on their own looks like p(t)=3*t^2-2*t^3 and map t : [0,1] -> [0,1] in a smooth way. Then integrate the poly to P(t)=t^3-0.5*t^4 which gives the "BLEP for the first sample" and then the second sample is 1-P(1-t) for symmetry, except we remove the trivial step, so just -P(1-t), which makes sense because.. well, it works out that way in math, it works and doesn't sound totally horrid.

Now we can go on with the basic oscillator, which needs a bit of state. So I wrapped it inside the minimum struct that makes it possible to test this thing, then copy-pasted it here.

I wrote the code for this (it's fresh), but I tested it and fixed the bugs.

It's only the COMMENTS that make it long, and those are rest of this tutorial. :)

I suggest copying to an editor that does syntax highlight.

Code: Select all

struct Oscillator {

// these are internal state, not parameters
float phase;
float blepDelay;
float widthDelay;
int   pulseStage;

// convenience constructor
Oscillator() { reset(); }

// the initialization code
void reset()
{
   // set current phase to 0 
   phase = 0;
   // blep delay "buffer" to zero
   blepDelay = 0;
   // previous pulseWidth to any valid value
   widthDelay = 0.5f;
   // only 2 stages here, set to first
   pulseStage = 0;
}

// the freq should be normalized realFreq/samplerate 
// the mix parameter is 0 for saw, 1 for pulse
// the pwm is audio-rate buffer, since that needs treatment
// rest are pretty trivial to turn into audio-rate sources
void render(float freq, float mix, float * pwm,
            float * bufOut, int nsamples)
{

   for(int i = 0; i < nsamples; ++i)
   {
      // the BLEP latency is 1 sample, so first
      // take the delayed part from previous sample
      float out = blepDelay;

      // then reset the delay so we can build into it
      blepDelay = 0;

      // then proceed like a trivial oscillator
      phase += freq;

      // load and sanity-check the new PWM
      // ugly things would happen outside range
      float pulseWidth = pwm[i];

      if(pulseWidth > 1) pulseWidth = 1;
      if(pulseWidth < 0) pulseWidth = 0;

      // Then replace the reset logic: loop until we
      // can't find the next discontinuity, so during
      // one sample we can process any number of them!
      while(true)
      {
         // Now in order of the stages of the wave-form
         // check for discontinuity during this sample.

         // First is the "pulse-width" transition.
         if(pulseStage == 0)
         {
            // if we didn't hit the pulse-width yet
            // break out of the while(true) loop
            if(phase < pulseWidth) break;

#ifdef NO_PWM_FIX
            // otherwise solve transition: when during
            // this sample did we hit the pw-border..
            // t = (1-x) from: phase + (x-1)*freq = pw
            float t = (phase - pulseWidth) / freq;
#else
            // the above version is fine when pw is constant
            // and it's what we use for the reset part, but
            // for pw modulation, t could end up outside [0,1]
            // and that will sound pretty ugly, so use lerp:
            //   phase + (x-1)*freq = (1-x)*pwOld + x*pw
            // and again t = (1 - x), and hopefully..
            float t = (phase - pulseWidth) 
                    / (widthDelay - pulseWidth + freq);
#endif
            // so then scale by pulse mix
            // and add the first sample to output
            out       += mix * poly3blep0(t);

            // and second sample to delay
            blepDelay += mix * poly3blep1(t);

            // and then we proceed to next stage
            pulseStage = 1;
         }

         // then whether or not we did transition, if stage
         // is at this point 1, we process the final reset
         if(pulseStage == 1)
         {
            // not ready to reset yet?
            if(phase < 1) break;

            // otherwise same as the pw, except threshold 1
            float t = (phase - 1) / freq;
            
            // and negative transition.. normally you would
            // calculate step-sizes for all mixed waves, but
            // here both saw and pulse go from 1 to 0.. so
            // it's always the same transition size!
            out       -= poly3blep0(t);
            blepDelay -= poly3blep1(t);

            // and then we do a reset (just one!)
            pulseStage = 0;
            phase -= 1;
         }

         // and if we are here, then there are possibly
         // more transitions to process, so keep going
      }

      // When the loop breaks (and it'll always break)
      // we have collected all the various BLEPs into our
      // output and delay buffer, so add the trivial wave
      // into the buffer, so it's properly delayed
      //
      // note: using pulseStage instead of pw-comparison
      // avoids inconsistencies from numerical inaccuracy
      blepDelay += (1-mix)*phase
                 + mix*(pulseStage ? 1.f : 0.f);

      // and output is just what we collected, but
      // let's make it range -1 to 1 instead
      bufOut[i] += (2*out - 1);

      // and store pulse-width delay for the lerp
      widthDelay = pulseWidth;

   }
}

}; // struct Oscillator

[/size]
The reason it's structure this way, is that it makes it a lot easier to understand and alot easier to edit and extend. It's also the structure that actually works, when you want to do something fancier than basic saws and pulse.

And the exact same render-method, with comments removed.. so it's easier to see the structure:

Code: Select all

void render(float freq, float mix, float * pwm,
            float * bufOut, int nsamples)
{

   for(int i = 0; i < nsamples; ++i)
   {
      float out = blepDelay;
      blepDelay = 0;

      phase += freq;

      float pulseWidth = pwm[i];

      if(pulseWidth > 1) pulseWidth = 1;
      if(pulseWidth < 0) pulseWidth = 0;

      while(true)
      {
         if(pulseStage == 0)
         {
            if(phase < pulseWidth) break;
#ifdef NO_PWM_FIX
            float t = (phase - pulseWidth) / freq;
#else
            float t = (phase - pulseWidth) 
                    / (widthDelay - pulseWidth + freq);
#endif
            out       += mix * poly3blep0(t);
            blepDelay += mix * poly3blep1(t);
            pulseStage = 1;
         }

         if(pulseStage == 1)
         {
            if(phase < 1) break;

            float t = (phase - 1) / freq;
            out       -= poly3blep0(t);
            blepDelay -= poly3blep1(t);
            pulseStage = 0;
            phase -= 1;
         }
      }
      blepDelay += (1-mix)*phase
                 + mix*(pulseStage ? 1.f : 0.f);
      bufOut[i] += (2*out - 1);
      widthDelay = pulseWidth;

   }
}
[/size]

Note that most of the time the loop does nothing: a few integer comparisons, one floating point comparison, and then break out. We only do extra work when there are discontinuities on the current sample, otherwise it's identical to naive!

If you really want to, you can use this with the same terms that I used for the filter that got more popular that it should... but really it's more of a clean starting point for something fancier (just like that filter was supposed to be).

Post

Oh and.. please let me know if the extension to "proper" BLEPs is something that I should cover in a future post. Essentially I'm not sure if it's trivial enough, or whether anyone even cares. Same thing with things like triangles: it's also quite simple, but I'm not going to bother if nobody cares in the first place.

Post

cheers, i'm out. no desire to make plugins for fussy people, happy with the osc methods i have, but ty :)
you come and go, you come and go. amitabha neither a follower nor a leader be tagore "where roads are made i lose my way" where there is certainty, consideration is absent.

Post

Interesting stuff, thanks for sharing. :)

Post

Really interesting, thanks Teemu ... nothing is trivial.

And thanks for all your contributions.

Post

wow, thanks mystran!

Post

Oh and forgot the audio demo:

http://www.signaldust.com/files/polyeb.wav .. rendered with the exact code as-is, with vsync-driven 8-bit control path (except I cheat and do audio-rate PWM with real sine here to better demonstrate how it sounds). This tune was composed by writing hex-bytes into C-code directly. Don't ask. :D

Observe that while it's not that great, the oscillator is still clean enough that you can clearly hear the control noise. Final fade is audio-editor though. :)

Post

Firstly, the osc code looks very sensible, thanks!

Secondly, though, the one thing that bothers me about BLEPs in general, is how to figure out the polynomial in the first place, and how to modify it to sound just even a tad different, a solution for a saw might work but not as well as something like, for example a shape like tanh(saw). Also, different waveshapes have essentially different algorithms - i.e. I've seen saw and tri BLEP code before and they differ quite a bit. That being said - I find your thread heading quite interisting - what is your definition - pertaining to your code, that would describe it as being extensible? Also, if you could explain the latter a bit more - that would be awesome.
mystran wrote: How I derived it (if you actually care): Build C1 continuous impulse from two cubic "smooth steps" which on their own looks like p(t)=3*t^2-2*t^3 and map t : [0,1] -> [0,1] in a smooth way. Then integrate the poly to P(t)=t^3-0.5*t^4 which gives the "BLEP for the first sample" and then the second sample is 1-P(1-t) for symmetry, except we remove the trivial step, so just -P(1-t), which makes sense because.. well, it works out that way in math, it works and doesn't sound totally horrid.
mystran wrote: This tune was composed by writing hex-bytes into C-code directly. Don't ask. Very Happy
+100 geek points - that's just insane! :clap:

Post

As for .. stuff that BLEPs can't do: they can solve the issue of aliasing discontinuities in (reasonable number of) finite derivatives. They won't solve all problems (not even all interesting ones), but they solve a large subset of discontinuity problems, and if you first cancel the naive-step of the tanh(saw) oscillator, you might already get quite significant quality improvement... or just generate oscillator, and shape it later.

As for extensibility: The triangle-code would be identical to the pulse-code... well fixed 0.5 duty cycle if you don't want PWM triangles.. but really you just change the actual BLEP functions to deal with derivative instead of step discontinuity (and scale by frequency because slopes and hence the discontinuities get steeper at higher frequencies.. and in general Nth derivative you'd scale by freq^N, though in practice going much further than 3 or 4 might nor be worth it). No changes in algorithm what-so-ever: just different naive wave-form and different BLEP (instead of "step type" we need "slope-change" type).

The code is structured as a state-machine. It's a trivial two-state machine here, but if you wanted to add more stages, or even do analytic ring-modulation of multiple waves, the state set and transition would get more complex, yet the same idea would work (well, the ring-mod would need quite a bit of work, but the idea works): start from the current state, see if there is a state-transition during this sample period. If there's a transition, then solve and insert corrections and repeat until next "transition" is the "end of current sample."

The fact that the transitions are handled in order and completely at once means that as soon as the transition has been processed, you can forget about it's existence.. and that's really powerful, because a version that actually tries to track "currently playing BLEPs" gets really messy if you have a lot of transitions that can randomly overlap; here's is irrelevant, and you don't need to worry about any of that, so adding features is a lot easier.

Somewhat less extreme example: handle regular symmetric triangle and varying pulse-width at the same time.. tracking the two in separate stage-variables is probably easier, but really the "full state-machine states" are just the valid combinations of those.

Anyway, I could add triangle, but I want to check it for correctness first then, so won't do it today.

As for sounding different: you can still put processing after the oscillator core. If you start with anti-aliased saw, then when you proceed to wave-shaping it, you can still get much better results than if you started with a naive saw.. even if you weren't technically "properly" anti-aliased anymore (and oversampling can now work much better with lower factor, since you only need to combat the aliasing of the shaper, rather than the whole thing!).

Post

Oh and the key concept: you choose exactly one type of BLEP for all "step type" discontinuities. Doesn't matter where that discontinuity is, if there's step, then it's the "step type" BLEP. Similarly for all slope changes.. you use "slope type" BLEP. And so on.

The reason I'm not posting polynomials to use for slope-type right here, right now, is that I want to check a few alternatives and see which ones results in the least aliasing in practice, so I don't end up recommending something stupid. :)

After all that polynomial I used for the steps, wasn't quite the first one I tried (which was still better than naive though, but only slightly).. it was only the second try that I got it right. ;)

Post

Ah, nice implementation - how much aliasing reduction does a one sample blep give you?

Cheers,
- Clemens

Post

cheppner wrote:Ah, nice implementation - how much aliasing reduction does a one sample blep give you?
Well, mm.. don't have a picture to post, but basically a bit cleaner than DPW. You can take the sound-file I posted and throw it at an analyzer, rendered 44.1kHz without oversampling, no filters either.

Also.. technically the BLEP is 2 samples and the true discontinuity lies between the two sampling instants. But yeah, it was more about the implementation structure than about PolyBLEPs specifically.. since the implementation works as-is with longer BLEPs. :)

Post

I have ported your tutorial code to REAPER JS, I'm sure (hope) you don't mind:

Code: Select all

desc:Tutorial: BLEPs (using PolyBLEPs but extensible)

slider1:440<20,20000,1>Frequency (Hz)
slider2:0.0<0.0,1.0,0.01>Waveform Mix
slider3:0.5<0.1,0.9,0.01>Pulse Width

@init

// The first sample BLEP
function poly3blep0(t) 
  local(t2)
( 
  // these are just sanity checks
  // correct code doesn't need them
  t < 0 ? 0 :
  t > 1 ? 1 : (

    t2 = t*t;
    t * t2 - 0.5 * t2 * t2;
  );
);

// And second sample as wrapper, optimize if you want.
function poly3blep1(t)
( 
  -poly3blep0(1-t);
);

// the initialization code
function reset()
  // these are internal state, not parameters
  instance(phase, blepDelay, widthDelay, pulseStage)
(
  // set current phase to 0 
  phase = 0; 
  // blep delay "buffer" to zero
  blepDelay = 0;
  // previous pulseWidth to any valid value
  widthDelay = 0.5;
  // only 2 stages here, set to first
  pulseStage = 0;
);

// the freq should be normalized realFreq/samplerate
// the mix parameter is 0 for saw, 1 for pulse
// the pwm is audio-rate buffer, since that needs treatment
// rest are pretty trivial to turn into audio-rate sources
function render(freq, mix, pwm)
  instance(phase, blepDelay, widthDelay, pulseStage)
  local(out, pulseWidth, t, break)
( 
  // the BLEP latency is 1 sample, so first
  // take the delayed part from previous sample
  out = blepDelay;

  // then reset the delay so we can build into it
  blepDelay = 0;

  // then proceed like a trivial oscillator
  phase += freq; 

  // load and sanity-check the new PWM
  // ugly things would happen outside range
  pulseWidth = pwm;

  pulseWidth > 1 ? pulseWidth = 1;
  pulseWidth < 0 ? pulseWidth = 0;

  // Then replace the reset logic: loop until we
  // can't find the next discontinuity, so during
  // one sample we can process any number of them!
  break = 0;
  while(
    // Now in order of the stages of the wave-form
    // check for discontinuity during this sample.

    // First is the "pulse-width" transition.
    pulseStage == 0 ? (
      // if we didn't hit the pulse-width yet
      // break out of the while(true) loop
      phase < pulseWidth ? break = 1 : (

/* #ifdef NO_PWM_FIX
        // otherwise solve transition: when during
        // this sample did we hit the pw-border..
        // t = (1-x) from: phase + (x-1)*freq = pw
        t = (phase - pulseWidth) / freq;
#else */
        // the above version is fine when pw is constant
        // and it's what we use for the reset part, but
        // for pw modulation, t could end up outside [0,1]
        // and that will sound pretty ugly, so use lerp:
        //   phase + (x-1)*freq = (1-x)*pwOld + x*pw
        // and again t = (1 - x), and hopefully..
        t = (phase - pulseWidth)
          / (widthDelay - pulseWidth + freq);
// #endif
        // so then scale by pulse mix
        // and add the first sample to output
        out       += mix * poly3blep0(t);

        // and second sample to delay
        blepDelay += mix * poly3blep1(t);

        // and then we proceed to next stage
        pulseStage = 1;
      );
    );

    // then whether or not we did transition, if stage
    // is at this point 1, we process the final reset
    !break && pulseStage == 1 ? (
      // not ready to reset yet? 
      phase < 1 ? break = 1 : (

        // otherwise same as the pw, except threshold 1
        t = (phase - 1) / freq;

        // and negative transition.. normally you would
        // calculate step-sizes for all mixed waves, but
        // here both saw and pulse go from 1 to 0.. so
        // it's always the same transition size!
        out       -= poly3blep0(t);
        blepDelay -= poly3blep1(t);

        // and then we do a reset (just one!)
        pulseStage = 0; 
        phase -= 1; 
      );
    );

    // and if we are here, then there are possibly
    // more transitions to process, so keep going
    !break;
  );

  // When the loop breaks (and it'll always break)
  // we have collected all the various BLEPs into our
  // output and delay buffer, so add the trivial wave
  // into the buffer, so it's properly delayed
  //
  // note: using pulseStage instead of pw-comparison
  // avoids inconsistencies from numerical inaccuracy
  blepDelay += (1-mix)*phase
                + mix*(pulseStage ? 1 : 0);

  // and output is just what we collected, but
  // let's make it range -1 to 1 instead
  out = (2*out - 1);

  // and store pulse-width delay for the lerp
  widthDelay = pulseWidth;

  out;
);

osc.reset();

@slider

gain = 0.25;
freq = slider1 / srate;
mix = slider2;
pwm = slider3;

@sample

spl0 = spl1 = gain * osc.render(freq, mix, pwm);
BTW, you could also use 0.5*t^2 instead of t^3 - 0.5*t^4. This will give you even less aliasing, although it does roll off a bit more near Nyquist:

Image
t^3 - 0.5*t^4

Image
0.5*t^2

Post

mystran wrote:Well, mm.. don't have a picture to post, but basically a bit cleaner than DPW.
Well, I do:

Image
DPW (order 2)

Post

Hi,
I would love to see a 44.1 wave file with a 1000 hz saw generated by this algo.
Thanks, Chris

EDIT: Not needed anymore, have generated by myself. :wink:

Post Reply

Return to “DSP and Plugin Development”