## Tutorial: BLEPs (using PolyBLEPs but extensible)

mystran
KVRAF

4949 posts since 11 Feb, 2006, from Helsinki, Finland
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 BLEPstatic 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);}`

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 parametersfloat phase;float blepDelay;float widthDelay;int   pulseStage;// convenience constructorOscillator() { reset(); }// the initialization codevoid 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 sourcesvoid 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`

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;   }}`

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).
<- plugins | forum
mystran
KVRAF

4949 posts since 11 Feb, 2006, from Helsinki, Finland
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.
<- plugins | forum
xoxos
Mr Entertainment

12018 posts since 29 Apr, 2002, from i might peeramid
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 xoxos.net free vst. neither a follower nor a leader be
where there is certainty, consideration is absent.
Tale
KVRist

478 posts since 12 Apr, 2010, from Lowlands of Holland
Interesting stuff, thanks for sharing.
Martinic Scanner Vibrato
peterst_de
KVRer

4 posts since 2 Dec, 2013
Really interesting, thanks Teemu ... nothing is trivial.

And thanks for all your contributions.
hibrasil
KVRian

764 posts since 23 Jun, 2002, from Huddersfield, UK
wow, thanks mystran!
My Website | WDL-OL | Web Audio Modules - WAMs | Oli Larkin Plugin's Facebook
Available for Audio Dev tuition via Skype (IPlug/JUCE/C++)
mystran
KVRAF

4949 posts since 11 Feb, 2006, from Helsinki, Finland
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.

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.
<- plugins | forum
KVRian

1033 posts since 8 Feb, 2012, from South - Africa
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!
mystran
KVRAF

4949 posts since 11 Feb, 2006, from Helsinki, Finland
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!).
<- plugins | forum
mystran
KVRAF

4949 posts since 11 Feb, 2006, from Helsinki, Finland
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.
<- plugins | forum
cheppner
KVRist

172 posts since 20 Aug, 2010
Ah, nice implementation - how much aliasing reduction does a one sample blep give you?

Cheers,
- Clemens
mystran
KVRAF

4949 posts since 11 Feb, 2006, from Helsinki, Finland
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.
<- plugins | forum
Tale
KVRist

478 posts since 12 Apr, 2010, from Lowlands of Holland
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 Mixslider3:0.5<0.1,0.9,0.01>Pulse Width@init// The first sample BLEPfunction 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 codefunction 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 sourcesfunction 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();@slidergain = 0.25;freq = slider1 / srate;mix = slider2;pwm = slider3;@samplespl0 = 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:

t^3 - 0.5*t^4

0.5*t^2
Martinic Scanner Vibrato
Tale
KVRist

478 posts since 12 Apr, 2010, from Lowlands of Holland
mystran wrote:Well, mm.. don't have a picture to post, but basically a bit cleaner than DPW.

Well, I do:

DPW (order 2)
Martinic Scanner Vibrato
Chris-S
KVRAF

2562 posts since 10 Nov, 2013, from Germany
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.
Next

Moderator: Moderators (Main)