Tutorial: BLEPs (using PolyBLEPs but extensible)

[mystran]

Thanks! I definitely will try this, when I'll give longer Bleps a try. The workflow you write down sounds like convolution, but with an IR that is dependent on the size of the phase resets.

Well, most of the time you would do convolution in "gather" mode, where you multiply each input sample with one kernel tap, add them all together and that's your output sample.

But we can transpose the whole thing into "scatter" where we take one input sample, multiply the whole kernel with that sample and then add all these taps into the output buffer one by one.

For regular filtering, gather is generally preferable as more efficient (with some special cases where scatter might make sense), but for BLEPs the scatter convolution is perfect, because even though the input sampling rate is theoretically infinite (ie. continuous time), only the occasional "input sample" here and there actually produces non-zero residue... so BLEPs are convolution really, just in a very optimized form.

[9b0]

My previous code had still issues. I was not checking for the case, when both oscillators reset in between the same two samples correctly. Now there's no wrong BLEPS, the only noise that appears is really just aliasing when the master frequency is really high, but this is already acceptable, and better than any hard-sync algorithm I was able to find online.

Sound demo:

LUA code (with comments for each line, and also the polyBlep function):

Code: Select all

local function dsp_oscSaw(self,F1,F2)

    local function polyBlep(blep,d,a,u)
        local pos={1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8}
        local d2=d*d
        local d3=d*d*d
        local d4=d*d*d*d
        blep[pos[a  ]]=(blep[pos[a  ]]           + 0.00029*d + 0.00737*d2 + 0.00617*d3 + 0.03817*d4)*u
        blep[pos[a+1]]=(blep[pos[a+1]] + 0.05254 + 0.18783*d + 0.22663*d2 + 0.14955*d3 - 0.11656*d4)*u
        blep[pos[a+2]]=(blep[pos[a+2]] - 0.5     + 0.62351*d + 0.02409*d2 - 0.31670*d3 + 0.11656*d4)*u
        blep[pos[a+3]]=(blep[pos[a+3]] - 0.05254 + 0.18839*d - 0.25816*d2 + 0.16102*d3 - 0.03871*d4)*u
    end

    local F={F1,F2}
    local pmax=phase[2] -- remember phase 2's maximum value
    for i=1,2 do
        inc[i]=min(F[i]*sRR,0.5) -- calculate the incremental
        phase[i]=phase[i]+inc[i] -- update the phase
        flip[i]=math.floor(phase[i]) -- if phase>=1 then flip=1
    end

    if flip[2]==1 then
        phase[2]=phase[2]-1 -- we reset osc2
        d[2]=phase[2]/inc[2] -- calculate the distance from the current sample and where the function crossed 1 (how far we are into the discontinuity)
        if flip[1]==1 then -- if both osc's reset inbetween the same 2 samples
            d[1]=(phase[1]-1)/inc[1] -- calculate the distance of osc1's discontinuity from the current sample
            if d[1]<d[2] then -- if osc1 resets after osc2 (the distance is smaller), we have to calculate how much osc1 travels after the reset
                polyBlep(blep[2],d[2],blepIndex,1) -- calculate the blep for the reset
                pmax=d[1]*inc[2] -- this is going to be the phase value where osc1 will reset osc2 again (in the same sample)
                phase[2]=d[1]*inc[2] -- this is the phase of osc2 at the intersample discontinuity
            end
        else
            polyBlep(blep[2],d[2],blepIndex,1) -- calculate the blep
        end
    end
    if flip[1]==1 then -- if the slave did not flip, but the master did
        phase[1]=phase[1]-1 -- reset phase 1
        d[1]=phase[1]/inc[1] -- calculate the intersample position of the phase crossing 1
        phase[2]=d[1]*inc[2] -- reset phasse 2 based on the new value of phase 1
        scale=pmax-phase[2]+inc[2] -- calculate the scaling factor for the blep based on the new value of phase 2
        polyBlep(blep[1],d[1],blepIndex,scale) -- calculate the blep
    end
    y=z[2]-blep[1][blepIndex]-blep[2][blepIndex] -- calculate the output
    for i=1,2 do
        blep[i][blepIndex]=0 -- reset the blep
        flip[i]=0 -- set flips to 0 (who knows)
    end
    z[2]=z[1] -- sample delay
    z[1]=phase[2] -- another delay
    blepIndex=(blepIndex%8)+1 -- increment the blep index
    return y*2-1
end

[mystran]

Now there's no wrong BLEPS, the only noise that appears is really just aliasing when the master frequency is really high, but this is already acceptable, and better than any hard-sync algorithm I was able to find online.

Something to keep in mind is that if you allow the oscillator frequency to get high enough, then you can have multiple BLEP transitions within the same sampling interval even without hardsync. With a straight saw-wave it only really happens once the wavelength is less than the sampling interval (that is, frequency is higher than the sampling rate), but this can still be an issue with hard-sync slave and more so with narrow pulse widths and such (where the pulse can go high, then low, or vice versa within a single sampling interval).

This is still solvable, you basically just have to structure your code in such a way that it solves for the next transition in a loop (and keeps adding arbitrary number of BLEPs on top of each other) until there can be no more transition during the current sampling interval, at which point with a good enough BLEP kernel you should be able to sweep the oscillators past Nyquist and have them completely filtered out (well, down to whatever attenuation your filter can do). Even then, it probably makes sense to cap the frequencies to some value anyway in order to make sure that the CPU cost of having to insert dozens of BLEPs within a single sampling interval just to produce silence doesn't get completely out of hand.

[Urs]

Hehehe, reading this I realise how very unaccustomed I have become to terminology using the master and slave analogy for such kind of relation. I usually use "MainX" and just "X", sometimes something akin to "SecondaryX" or "FollowingX" if it needs to be made clear, and I have not experienced any downsides of doing so. Like, doing so wasn't an ordeal, while others - hopefully - still know what I'm talking about.

That said, I fully appreciate what you guys write, and if I ever find the time to get back to the topic, I'd love to share my own implementation thoughts on it.

[mystran]

Hehehe, reading this I realise how very unaccustomed I have become to terminology using the master and slave analogy for such kind of relation.

I think hard-sync is a voluntary relationship between consenting adult oscillators.

[9b0]

Now there's no wrong BLEPS, the only noise that appears is really just aliasing when the master frequency is really high, but this is already acceptable, and better than any hard-sync algorithm I was able to find online.

Something to keep in mind is that if you allow the oscillator frequency to get high enough, then you can have multiple BLEP transitions within the same sampling interval even without hardsync.

Thanks! I was aware of this, and I put a min function in front of my incremented values, so they can't get higher than 0.5.

I've implemented the lookup table method for the bleps, so when a polyBlep is called, I'm not calculating the polynomials anymore, refer to a table. I've created tables for 2 points, 4 points, 6 points, and 8 points since this page made implementing the residuals dead easy (I found this yesterday): https://ryukau.github.io/filter_notes/p ... idual.html
- Thanks to the owner of the Git

I've also allowed negative frequencies now, and everything seems to work fine. I'm worried a bit, that maybe I'm scaling the Blep for the osc1 reset wrong a bit, but not sure if it's only a normal amount of aliasing, or me being dumb.

[9b0]

Hehehe, reading this I realise how very unaccustomed I have become to terminology using the master and slave analogy for such kind of relation. I usually use "MainX" and just "X", sometimes something akin to "SecondaryX" or "FollowingX" if it needs to be made clear, and I have not experienced any downsides of doing so. Like, doing so wasn't an ordeal, while others - hopefully - still know what I'm talking about.

That said, I fully appreciate what you guys write, and if I ever find the time to get back to the topic, I'd love to share my own implementation thoughts on it.

Thanks, Urs! I've felt pretty uncomfortable from this too. Still, I was used to the terminology... in my latest code I did reference the oscillators and phases and such with indexing already (and got used to it).

[Gordonjcp]

Something to keep in mind is that if you allow the oscillator frequency to get high enough, then you can have multiple BLEP transitions within the same sampling interval ... more so with narrow pulse widths and such (where the pulse can go high, then low, or vice versa within a single sampling interval).

So for example in the "Alpha Juno" code, if I wanted to antialias the PWM gate where it's extremely likely that you can miss a very narrow pulse (because at 48kHz you start to increment the counter by more than 1 every sample at the low hundreds of Hz output frequency), you'd need to check if the counter wrapped and then see if there was a need for a "blep up" before putting in the "blep down" again?

I guess then I'd detect if it had wrapped while the pulse was less than the phase increment, and then reverse back over my own tracks to find the "up" within the sample?

[mystran]

Something to keep in mind is that if you allow the oscillator frequency to get high enough, then you can have multiple BLEP transitions within the same sampling interval ... more so with narrow pulse widths and such (where the pulse can go high, then low, or vice versa within a single sampling interval).

So for example in the "Alpha Juno" code, if I wanted to antialias the PWM gate where it's extremely likely that you can miss a very narrow pulse (because at 48kHz you start to increment the counter by more than 1 every sample at the low hundreds of Hz output frequency), you'd need to check if the counter wrapped and then see if there was a need for a "blep up" before putting in the "blep down" again?

IMHO the best way to handle BLEP oscillators is to treat them as state machines where every sample you check for the next state transition and then if there is one, repeat until there is none.

For example, if we have a pulse wave, then we might have states "low" and "high" with "low" being the initial state after a reset. Every sample you increment the phase and check if it has crossed the current PWM threshold. If the threshold has been crossed, that's a state transition, you solve for the exact transition time and insert a BLEP, then set "high" as the current oscillator state.

Then we repeat the same thing for the new state, this time checking if we're past the point where the phase should be reset. If that threshold has been crossed as well, then again we solve for the exact time offset and insert a BLEP, subtract the reset threshold from the phase and set "low" as the current oscillator state.. and then we repeat, again checking if the PWM threshold has been met again.. and so on, until we get to the point where we observe that there aren't any more transitions during the current sample.

Next sample, we look at which state we were left in and carry on doing the same thing we did the previous sample, again looping until there are no more transitions. Now, for most samples (at reasonable oscillator frequency) there's no transitions at all. Sometimes there's a single transition. But if we have a highish frequency and narrow pulse-width, we might get two transitions. If we allow a ridiculous frequency, then we might have 500 transitions; that'd obviously be very slow to process, so you want to cap the frequency at some sensible maximum, but the point is.. if you loop for "are we at the next transition already" until there's no more transitions for the current sample, it "just works" and with large enough kernel you can run at 44.1kHz and have the oscillator at 1MHz and it's outputting a crapton of BLEPs every sample and the output is the expected silence (and in practice a bit of noise).

The code in the very first post of this thread uses this strategy: there's a while(true) loop for solving BLEPs and in each stage of the oscillator, we break out of the loop if the next transition cannot be during the current sample. It doesn't care about how many BLEPs we need, it just keeps inserting them until it's gone through all the possible transitions and it just works(tm).

For hard-sync, the only complication is that every time you have a "natural" transition, you check whether that natural transition or the sync transition comes first. If it's the natural transition, you proceed normally. If it's the sync transition, you perform sync, then check for natural transitions again after the reset.

[Gordonjcp]

The code in the very first post of this thread uses this strategy: there's a while(true) loop for solving BLEPs and in each stage of the oscillator, we break out of the loop if the next transition cannot be during the current sample. It doesn't care about how many BLEPs we need, it just keeps inserting them until it's gone through all the possible transitions and it just works(tm).

I'm going to have to go and stare at that some more until I get it back into my brain again, probably tomorrow because it's too near midnight for that sort of thing.

How does "know" that a transition happens mid-sample then? Assume the pulse width is set to 127, so it's one "clock" wide. If the phase accumulator is running slowly enough that one "clock" is always larger than a sample then it has to play that 1-clock (really 1-sample) click of the narrowest possible pulse. But at some point it may skip right over - the counter may go from 254 to 2, say, and the pulse ought to have happened somewhere in the middle, but we skipped entirely over it.

[mystran]

How does "know" that a transition happens mid-sample then? Assume the pulse width is set to 127, so it's one "clock" wide. If the phase accumulator is running slowly enough that one "clock" is always larger than a sample then it has to play that 1-clock (really 1-sample) click of the narrowest possible pulse. But at some point it may skip right over - the counter may go from 254 to 2, say, and the pulse ought to have happened somewhere in the middle, but we skipped entirely over it.

You don't skip over anything if you always check against the next expected transition.

At the beginning of the sample, you know the current state and you know the current phase. Based on the current state, you know the next expected transition. After incrementing the phase, compare the phase at the end of the sample with the next expected transition; if it's larger, then solve backwards to figure out when we crossed the transition threshold. Process transition and repeat.

Square wave example (with phase ranging [0,1] for simplicity). If we're in state 0 (low) then we know that the next transition (ignoring sync here) will be the transition to state 1 (high) where the phase increases past 0.5 and after incrementing the phase, we just compare whether phase is now larger than 0.5? Suppose that's the case, we solve how large a fraction of the phase increment did we actually need to hit the 0.5 exactly, that fraction is the sub-sample position for the transition where you insert BLEP. Great, we're in state 1 now. Next we compare the phase against the next expected transition at 1 (=natural reset). Suppose our phase is also larger than that, so again we solve exactly how much of the phase increment did we need to hit that value exactly, then we insert another BLEP there. Because this is a reset, we also subtract the threshold (=1 here) from the phase value at the end of the sample and go back to state 0. The next transition is 0.5 again, so another comparison .. and so on..

At some point (most of the time that's the very first check) we conclude that "we're not yet at the next transition" and we break out of the loop. With all the BLEPs already processed, all you then need to do is sample the naive waveform (ie. -1 if state 0, +1 if state 1, don't even need to look at the phase for a pulse). The naive waveform might "skip" some stages, but that's fine, we did all the BLEPs for the transitions and that's what matters.

ps. Note that if your "state" is not a simple sequence of stages, there might be multiple possible "next" transitions, but that's fine: you just have to figure out which one is the next. If those transitions can move, then sometimes you have to solve a bunch of stuff to figure out which one happens first, worst-case solve backwards for multiple possible transitions and compare them against each other.. but as long as you only ever allow the oscillator to move from one polynomial segment to a different polynomial segment through the earliest possible transition, you'll be ok.

[jokona]

AUTO-ADMIN: Non-MP3, WAV, OGG, SoundCloud, YouTube, Vimeo, Twitter and Facebook links in this post have been protected automatically. Once the member reaches 5 posts the links will function as normal.

This is already acceptable, and better than any hard-sync algorithm I was able to find online.

Hey and greetings from Hungary,
I translated your LUA code to FAUST to enhance the available resources online on hardsync with polyblep. The advantage is that you can try it online by pasting the code to faustide.grame.fr/ .
ATTENTION: Before pressing the start icon to run the code, set the bottom editor from the Diagram tab to the Plot tab. The Diagram tab takes ages to compile with this code, but if you set it to the Plot tab, the code should compile fast.

Another advantage is that you can view a spectrogram of the output. To configure the spectrogram after running the code, I recommend setting the plot mode to continuous, the samples to, say, 8192, the fft size to 4096, and keeping the default fft overlap of 2. Don't forget to tick the Draw spectrogram tickbox. To actually view the spectrogram, navigate to the top left corner of the Plot tab at the bottom of the editor, which should say Oscilloscope by default. No matter what it says, click the text to cycle to different plot formats until you get to the spectrogram. (If you did not tick the draw spectrogram tickbox, you only get a spectroscope plot but not a spectrogram.) I recommend playing with the freqs and observing what happens in the spectrogram.

Let me paste the code and then add a few comments:

Code: Select all (#)

import("stdfaust.lib");
SR = fconstant(int fSamplingFreq, <math.h>);

freq1 = hslider("master", 240, 20, 20000,0.01);
freq2 = hslider("slave", 240, 20, 20000,0.01);
// *** THIS IS THE DSP. SINCE THE ORGINAL CODE IS STATEFUL, WHICH IS NOT EASY TO CODE IN FAUST PER SE, WE RECURSE OVER THE MAIN LOOP WITH THE STATES EXPLICITELY PASSED AS INPUT AND RETURNED AS OUTPUT. THIS WAY THE STATES ARE "MUTABLE". ***
process =  ((freq1, freq2, _,_,_,_,_, _,_,_,_,_, _,_,_,_,_, _,_,_,_,_, _ : main_loop) ~ (_,_,_,_,_, _,_,_,_,_, _,_,_,_,_, _,_,_,_,_, _)) : !,!,!,!,!, !,!,!,!,!, !,!,!,!,!, !,!,!,!,!, !,_;

// *** HERE IS THE MAIN CODE. I DEVIDED IT INTO SECTIONS. SEE BELOW FOR CODE OF THE SECTIONS ***
//params : freq1, freq2
//+ states : phase1, phase2, z1, z2, blepindex, blep1_0 ... blep1_7, blep2_0 ... blep2_7. Note that input states get prefix in and output states get prefix out as things like phase = phase + 1; don't work in faust. Instead we use outphase = inphase + 1; to mimic mutability. There is a lot of clutter due to having to do this.
//+ out  : 2*y-1
main_loop(freq1, freq2, inphase1, inphase2, inz1, inz2, inblepindex, inblep1_0, inblep1_1, inblep1_2, inblep1_3, inblep1_4, inblep1_5, inblep1_6, inblep1_7, inblep2_0, inblep2_1, inblep2_2, inblep2_3, inblep2_4, inblep2_5, inblep2_6, inblep2_7) =  outphase1, outphase2, outz1, outz2, outblepindex, outblep1_0, outblep1_1, outblep1_2, outblep1_3, outblep1_4, outblep1_5, outblep1_6, outblep1_7, outblep2_0, outblep2_1, outblep2_2, outblep2_3, outblep2_4, outblep2_5, outblep2_6, outblep2_7, out with
{
    // SECTION A
    out_section_a = freq1, freq2, inphase1, inphase2 : section_a;
    a_pmax      = out_section_a : _,!,!,!,!,!,!;
    a_phase1    = out_section_a : !,_,!,!,!,!,!;
    a_phase2    = out_section_a : !,!,_,!,!,!,!;
    a_inc1      = out_section_a : !,!,!,_,!,!,!;
    a_inc2      = out_section_a : !,!,!,!,_,!,!;
    a_flip1     = out_section_a : !,!,!,!,!,_,!;
    a_flip2     = out_section_a : !,!,!,!,!,!,_;

    a_blep2_0 = (inblepindex, inblep2_0, inblep2_1, inblep2_2, inblep2_3, inblep2_4, inblep2_5, inblep2_6, inblep2_7 : selectblep);
    a_blep2_1 = (index, inblep2_0, inblep2_1, inblep2_2, inblep2_3, inblep2_4, inblep2_5, inblep2_6, inblep2_7 : selectblep) with {index = inblepindex + 1 : %(8);};
    a_blep2_2 = (index, inblep2_0, inblep2_1, inblep2_2, inblep2_3, inblep2_4, inblep2_5, inblep2_6, inblep2_7 : selectblep) with {index = inblepindex + 2 : %(8);};
    a_blep2_3 = (index, inblep2_0, inblep2_1, inblep2_2, inblep2_3, inblep2_4, inblep2_5, inblep2_6, inblep2_7 : selectblep) with {index = inblepindex + 3 : %(8);};

    // SECTION B
    out_section_b = a_pmax, a_phase1, a_phase2, a_inc1, a_inc2, a_flip1, a_flip2, a_blep2_0, a_blep2_1, a_blep2_2, a_blep2_3 : section_b;
    b_pmax      = out_section_b : _,!,!,!,!,!;
    b_phase2    = out_section_b : !,_,!,!,!,!;
    b_blep2_0   = out_section_b : !,!,_,!,!,!;
    b_blep2_1   = out_section_b : !,!,!,_,!,!;
    b_blep2_2   = out_section_b : !,!,!,!,_,!;
    b_blep2_3   = out_section_b : !,!,!,!,!,_;

    b_phase1 = a_phase1;
    b_inc1   = a_inc1;
    b_inc2   = a_inc2;
    b_flip1  = a_flip1;
    b_flip2  = a_flip2;

    b_blep1_0 = (inblepindex, inblep1_0, inblep1_1, inblep1_2, inblep1_3, inblep1_4, inblep1_5, inblep1_6, inblep1_7 : selectblep);
    b_blep1_1 = (index, inblep1_0, inblep1_1, inblep1_2, inblep1_3, inblep1_4, inblep1_5, inblep1_6, inblep1_7 : selectblep) with {index = inblepindex + 1 : %(8);};
    b_blep1_2 = (index, inblep1_0, inblep1_1, inblep1_2, inblep1_3, inblep1_4, inblep1_5, inblep1_6, inblep1_7 : selectblep) with {index = inblepindex + 2 : %(8);};
    b_blep1_3 = (index, inblep1_0, inblep1_1, inblep1_2, inblep1_3, inblep1_4, inblep1_5, inblep1_6, inblep1_7 : selectblep) with {index = inblepindex + 3 : %(8);};
    
    // SECTION C
    out_section_c = b_pmax, b_phase1, b_phase2, b_inc1, b_inc2, b_flip1, b_blep1_0, b_blep1_1, b_blep1_2, b_blep1_3 : section_c;
    c_phase1  = out_section_c : _,!,!,!,!,!;
    c_phase2  = out_section_c : !,_,!,!,!,!;
    c_blep1_0 = out_section_c : !,!,_,!,!,!;
    c_blep1_1 = out_section_c : !,!,!,_,!,!;
    c_blep1_2 = out_section_c : !,!,!,!,_,!;
    c_blep1_3 = out_section_c : !,!,!,!,!,_;
    
    c_blep2_0 = b_blep2_0;
    c_blep2_1 = b_blep2_1;
    c_blep2_2 = b_blep2_2;
    c_blep2_3 = b_blep2_3;

    bleps1 = (inblepindex, c_blep1_0, c_blep1_1, c_blep1_2, c_blep1_3, inblep1_0, inblep1_1, inblep1_2, inblep1_3, inblep1_4, inblep1_5, inblep1_6, inblep1_7 : set4bleps);
    bleps2 = (inblepindex, c_blep2_0, c_blep2_1, c_blep2_2, c_blep2_3, inblep2_0, inblep2_1, inblep2_2, inblep2_3, inblep2_4, inblep2_5, inblep2_6, inblep2_7 : set4bleps);

    blep1_0 = bleps1 : _,!,!,!,!,!,!,!;
    blep1_1 = bleps1 : !,_,!,!,!,!,!,!;
    blep1_2 = bleps1 : !,!,_,!,!,!,!,!;
    blep1_3 = bleps1 : !,!,!,_,!,!,!,!;
    blep1_4 = bleps1 : !,!,!,!,_,!,!,!;
    blep1_5 = bleps1 : !,!,!,!,!,_,!,!;
    blep1_6 = bleps1 : !,!,!,!,!,!,_,!;
    blep1_7 = bleps1 : !,!,!,!,!,!,!,_;

    blep2_0 = bleps2 : _,!,!,!,!,!,!,!;
    blep2_1 = bleps2 : !,_,!,!,!,!,!,!;
    blep2_2 = bleps2 : !,!,_,!,!,!,!,!;
    blep2_3 = bleps2 : !,!,!,_,!,!,!,!;
    blep2_4 = bleps2 : !,!,!,!,_,!,!,!;
    blep2_5 = bleps2 : !,!,!,!,!,_,!,!;
    blep2_6 = bleps2 : !,!,!,!,!,!,_,!;
    blep2_7 = bleps2 : !,!,!,!,!,!,!,_;

    y = inz2 - (inblepindex, blep1_0, blep1_1, blep1_2, blep1_3, blep1_4, blep1_5, blep1_6, blep1_7 : selectblep) - (inblepindex, blep2_0, blep2_1, blep2_2, blep2_3, blep2_4, blep2_5, blep2_6, blep2_7 : selectblep); // calculate the output

    outphase1 = c_phase1;
    outphase2 = c_phase2;
    outz2 = inz1;                       // sample delay
    outz1 = c_phase2;                   // another delay
    outblepindex = (inblepindex+1) % 8; // increment the blep index
    
    outbleps1 = (inblepindex, 0, blep1_0, blep1_1, blep1_2, blep1_3, blep1_4, blep1_5, blep1_6, blep1_7 : setblep);
    outbleps2 = (inblepindex, 0, blep2_0, blep2_1, blep2_2, blep2_3, blep2_4, blep2_5, blep2_6, blep2_7 : setblep);
    
    outblep1_0 = outbleps1 : _,!,!,!,!,!,!,!;
    outblep1_1 = outbleps1 : !,_,!,!,!,!,!,!;
    outblep1_2 = outbleps1 : !,!,_,!,!,!,!,!;
    outblep1_3 = outbleps1 : !,!,!,_,!,!,!,!;
    outblep1_4 = outbleps1 : !,!,!,!,_,!,!,!;
    outblep1_5 = outbleps1 : !,!,!,!,!,_,!,!;
    outblep1_6 = outbleps1 : !,!,!,!,!,!,_,!;
    outblep1_7 = outbleps1 : !,!,!,!,!,!,!,_;

    outblep2_0 = outbleps2 : _,!,!,!,!,!,!,!;
    outblep2_1 = outbleps2 : !,_,!,!,!,!,!,!;
    outblep2_2 = outbleps2 : !,!,_,!,!,!,!,!;
    outblep2_3 = outbleps2 : !,!,!,_,!,!,!,!;
    outblep2_4 = outbleps2 : !,!,!,!,_,!,!,!;
    outblep2_5 = outbleps2 : !,!,!,!,!,_,!,!;
    outblep2_6 = outbleps2 : !,!,!,!,!,!,_,!;
    outblep2_7 = outbleps2 : !,!,!,!,!,!,!,_;

    out   = y*2-1;
};

// params: inblep0, ... inblep3, d, u
// + out: outblep0, ... outblep3
// note: select pos based on "a" (a, a+1, a+2, a+3) %8
polyBlep(inblep0, inblep1, inblep2, inblep3, d, u) = outblep0, outblep1, outblep2, outblep3 with
{
    d2=d*d;
    d3=d*d*d;
    d4=d*d*d*d;
    outblep0 =(inblep0           + 0.00029*d + 0.00737*d2 + 0.00617*d3 + 0.03817*d4)*u;
    outblep1 =(inblep1 + 0.05254 + 0.18783*d + 0.22663*d2 + 0.14955*d3 - 0.11656*d4)*u;
    outblep2 =(inblep2 - 0.5     + 0.62351*d + 0.02409*d2 - 0.31670*d3 + 0.11656*d4)*u;
    outblep3 =(inblep3 - 0.05254 + 0.18839*d - 0.25816*d2 + 0.16102*d3 - 0.03871*d4)*u;
};


// *** I DEVIDED THE CODE INTO SMALL SECTIONS TO KEEP MYSELF SANE. (HARD TO TRANSLATE STATEFUL CODE TO THE FUNCTIONAL FAUST.) HERE ARE THE SECTIONS.***

// SECTION A
// params: freq1, freq2, inphase1, inphase2
// out: pmax, outphase1, outphase2, inc1, inc2, flip1, flip2
section_a(freq1, freq2, inphase1, inphase2) = pmax, outphase1, outphase2, inc1, inc2, flip1, flip2 with
{
    pmax = inphase2; // remember phase 2's maximum value

    inc1 = freq1/SR, 0.5 : min;   // calculate the incremental
    outphase1 = inphase1 + inc1;  // update the phase
    flip1 = outphase1 : floor;    // if phase>=1 then flip=1
    
    inc2 = freq2/SR, 0.5 : min;   // calculate the incremental
    outphase2 = inphase2 + inc2;  // update the phase
    flip2 = outphase2 : floor;    // if phase>=1 then flip=1
};

// SECTION B
// params: inpmax, inphase1, inphase2, inc1, inc2, flip1, flip2, inblep2_0, inblep2_1, inblep2_2, inblep2_3
// out: outpmax, outphase2, outblep2_0, outblep2_1, outblep2_2, outblep2_3
section_b(inpmax, inphase1, inphase2, inc1, inc2, flip1, flip2, inblep2_0, inblep2_1, inblep2_2, inblep2_3) = outpmax, outphase2, outblep2_0, outblep2_1, outblep2_2, outblep2_3 with
{
    cond_flip2 = flip2 == 1;
    case_flip2 = outpmax, outphase2, outblep2_0, outblep2_1, outblep2_2, outblep2_3 with
    {
        temp0_outphase2 = inphase2 - 1;  // we reset osc2
        d2 = temp0_outphase2 / inc2; // calculate the distance from the current sample and where the function crossed 1 (how far we are into the discontinuity)

        cond_flip1 = flip1 == 1;
        case_flip1 =  outpmax, outphase2, outblep2_0, outblep2_1, outblep2_2, outblep2_3 with
        {
            d1 = (inphase1-1)/inc1;   // calculate the distance of osc1's discontinuity from the current sample
            cond_d1_sm_d2 = d1 < d2;  // if osc1 resets after osc2 (the distance is smaller), we have to calculate how much osc1 travels after the reset
            case_d1_sm_d2 =  outpmax, outphase2, outblep2_0, outblep2_1, outblep2_2, outblep2_3 with
            {
                polyblep_out = inblep2_0, inblep2_1, inblep2_2, inblep2_3, d2, 1 : polyBlep;
                outblep2_0 = polyblep_out : _,!,!,!;
                outblep2_1 = polyblep_out : !,_,!,!;
                outblep2_2 = polyblep_out : !,!,_,!;
                outblep2_3 = polyblep_out : !,!,!,_;

                outpmax   = d1*inc1;               // this is going to be the phase value where osc1 will reset osc2 again (in the same sample)
                outphase2 = d1*inc2;               // this is the phase of osc2 at the intersample discontinuity
            };

            case_not_d1_sm_d2 = outpmax, outphase2, outblep2_0, outblep2_1, outblep2_2, outblep2_3 with
            {
                outpmax = inpmax;
                outphase2 = temp0_outphase2;
                outblep2_0 = inblep2_0;
                outblep2_1 = inblep2_1;
                outblep2_2 = inblep2_2;
                outblep2_3 = inblep2_3;
            };

            out_case_d1_sm_d2 = case_not_d1_sm_d2, case_d1_sm_d2 : zip6 : par(i, 6, select2(cond_d1_sm_d2));
            
            outpmax    = out_case_d1_sm_d2 : _,!,!,!,!,!;
            outphase2  = out_case_d1_sm_d2 : !,_,!,!,!,!;
            outblep2_0 = out_case_d1_sm_d2 : !,!,_,!,!,!;
            outblep2_1 = out_case_d1_sm_d2 : !,!,!,_,!,!;
            outblep2_2 = out_case_d1_sm_d2 : !,!,!,!,_,!;
            outblep2_3 = out_case_d1_sm_d2 : !,!,!,!,!,_;
        };
        case_not_flip1 = outpmax, outphase2, outblep2_0, outblep2_1, outblep2_2, outblep2_3 with
        {
            polyblep_out = inblep2_0, inblep2_1, inblep2_2, inblep2_3, d2, 1 : polyBlep;
            outblep2_0 = polyblep_out : _,!,!,!;
            outblep2_1 = polyblep_out : !,_,!,!;
            outblep2_2 = polyblep_out : !,!,_,!;
            outblep2_3 = polyblep_out : !,!,!,_;

            outpmax   = inpmax;
            outphase2 = temp0_outphase2;
        };

        out_case_flip1 = case_not_flip1, case_flip1 : zip6 : par(i, 6, select2(cond_flip1));
        
        outpmax    = out_case_flip1 : _,!,!,!,!,!;
        outphase2  = out_case_flip1 : !,_,!,!,!,!;
        outblep2_0 = out_case_flip1 : !,!,_,!,!,!;
        outblep2_1 = out_case_flip1 : !,!,!,_,!,!;
        outblep2_2 = out_case_flip1 : !,!,!,!,_,!;
        outblep2_3 = out_case_flip1 : !,!,!,!,!,_;

    };
    case_not_flip2 = outpmax, outphase2, outblep2_0, outblep2_1, outblep2_2, outblep2_3 with
    {
        outpmax = inpmax;
        outphase2 = inphase2;
        outblep2_0 = inblep2_0;
        outblep2_1 = inblep2_1;
        outblep2_2 = inblep2_2;
        outblep2_3 = inblep2_3;
    };

    out_case_flip2 = case_not_flip2, case_flip2 : zip6 : par(i, 6, select2(cond_flip2));    
    outpmax    = out_case_flip2 : _,!,!,!,!,!;
    outphase2  = out_case_flip2 : !,_,!,!,!,!;
    outblep2_0 = out_case_flip2 : !,!,_,!,!,!;
    outblep2_1 = out_case_flip2 : !,!,!,_,!,!;
    outblep2_2 = out_case_flip2 : !,!,!,!,_,!;
    outblep2_3 = out_case_flip2 : !,!,!,!,!,_;
};

// SECTION C
// params:pmax, inphase1, inphase2, inc1, inc2, flip1, inblep1_0, inblep1_1, inblep1_2, inblep1_3
// out: outphase2, outblep1_0, outblep1_1, outblep1_2, outblep1_3
section_c(pmax, inphase1, inphase2, inc1, inc2, flip1, inblep1_0, inblep1_1, inblep1_2, inblep1_3) = outphase1, outphase2, outblep1_0, outblep1_1, outblep1_2, outblep1_3 with
{
    cond_flip1 = flip1 == 1;
    case_flip1 = outphase1, outphase2, outblep1_0, outblep1_1, outblep1_2, outblep1_3 with
    {
        outphase1 = inphase1-1;                     // reset phase 1
        d1 = outphase1/inc1;                    // calculate the intersample position of the phase crossing 1
        outphase2 = d1*inc2;                    // reset phasse 2 based on the new value of phase 1
        scale = pmax - outphase2 + inc2;            // calculate the scaling factor for the blep based on the new value of phase 2
        polyblep_out = inblep1_0, inblep1_1, inblep1_2, inblep1_3, d1, scale : polyBlep;
        outblep1_0 = polyblep_out : _,!,!,!;
        outblep1_1 = polyblep_out : !,_,!,!;
        outblep1_2 = polyblep_out : !,!,_,!;
        outblep1_3 = polyblep_out : !,!,!,_;
    };

    case_not_flip1 = outphase1, outphase2, outblep1_0, outblep1_1, outblep1_2, outblep1_3 with
    {
        outphase1  = inphase1;
        outphase2  = inphase2;
        outblep1_0 = inblep1_0;
        outblep1_1 = inblep1_1;
        outblep1_2 = inblep1_2;
        outblep1_3 = inblep1_3;
    };

    out_case_flip1 = case_not_flip1, case_flip1 : zip6 : par(i, 6, select2(cond_flip1));
    
    outphase1  = out_case_flip1 : _,!,!,!,!,!;
    outphase2  = out_case_flip1 : !,_,!,!,!,!;
    outblep1_0 = out_case_flip1 : !,!,_,!,!,!;
    outblep1_1 = out_case_flip1 : !,!,!,_,!,!;
    outblep1_2 = out_case_flip1 : !,!,!,!,_,!;
    outblep1_3 = out_case_flip1 : !,!,!,!,!,_;
};

// UTIL FOR HANDLING ROUTING AND ENABLING ME TO HANDLE 
zip6(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6) = a1, b1, a2, b2, a3, b3, a4, b4, a5, b5, a6, b6;
selectblep(blepindex, blep0, blep1, blep2, blep3, blep4, blep5, blep6, blep7) = blep0, blep1, blep2, blep3, blep4, blep5, blep6, blep7 : ba.selectn(8, blepindex);
setblep(blepindex, inserted_value, inblep0, inblep1, inblep2, inblep3, inblep4, inblep5, inblep6, inblep7) = outblep0, outblep1, outblep2, outblep3, outblep4, outblep5, outblep6, outblep7 with
{
    outblep0 = inblep0, inserted_value : select2(blepindex == 0);
    outblep1 = inblep1, inserted_value : select2(blepindex == 1);
    outblep2 = inblep2, inserted_value : select2(blepindex == 2);
    outblep3 = inblep3, inserted_value : select2(blepindex == 3);
    outblep4 = inblep4, inserted_value : select2(blepindex == 4);
    outblep5 = inblep5, inserted_value : select2(blepindex == 5);
    outblep6 = inblep6, inserted_value : select2(blepindex == 6);
    outblep7 = inblep7, inserted_value : select2(blepindex == 7);
};
set4bleps(blepindex, val0, val1, val2, val3, inblep0, inblep1, inblep2, inblep3, inblep4, inblep5, inblep6, inblep7) = outblep0, outblep1, outblep2, outblep3, outblep4, outblep5, outblep6, outblep7 with
{

    idx(outindex) = 1*(blepindex == outindex) + 2*((blepindex + 1 % 8)  == outindex) + 3*((blepindex + 2 % 8)  == outindex) + 4*((blepindex + 3 % 8)  == outindex);

    outblep0 = inblep0, val0, val1, val2, val3 : ba.selectn(5, idx(0));
    outblep1 = inblep1, val0, val1, val2, val3 : ba.selectn(5, idx(1));
    outblep2 = inblep2, val0, val1, val2, val3 : ba.selectn(5, idx(2));
    outblep3 = inblep3, val0, val1, val2, val3 : ba.selectn(5, idx(3));
    outblep4 = inblep4, val0, val1, val2, val3 : ba.selectn(5, idx(4));
    outblep5 = inblep5, val0, val1, val2, val3 : ba.selectn(5, idx(5));
    outblep6 = inblep6, val0, val1, val2, val3 : ba.selectn(5, idx(6));
    outblep7 = inblep7, val0, val1, val2, val3 : ba.selectn(5, idx(7));
};

This code is not intended for reading, it is simply a translation of the LUA code from the post I reply to, into Faust. This way you can listen to it at https://faustide.grame.fr/ (https://faustide.grame.fr/) . The comments are from the original LUA code. The reason I recommend reading the OP's LUA code instead to understand the algorithm is that the original LUA code is stateful, where the states are mutable and change over time, whereas in Faust - which is a functional language - this is not possible per se. E.g. phase[1] = phase[1] -1 works in LUA, but we cannot write phase1 = phase1 -1; in faust and get the same result. To mimic mutable, stateful behaviour in this functional language I did something like outphase1 = inphase1 -1; etc rendering the code harder to read. The states are also explicitly passed to a recursion operator, further worsening readability. But again this is needed to mimic mutable, stateful code. Furthermore, in faust there are no arrays and no real "lazy evaluation" ifs, so I needed to hack around to mimic these (see the util section at the bottom if you really must). All in all, I recommend not to read the Faust code, read the LUA code from the OP if you want to understand the algorithm. Paste the Faust code into the IDE at https://faustide.grame.fr/ (https://faustide.grame.fr/) , however and play around with it!

[Marvinh]

What happens when we use C2 C3 etc .
More filtered, but the delay has to become longer ?


https://wikimedia.org/api/rest_v1/media ... 26afe6c289

[jokona]

What happens when we use C2 C3 etc .
More filtered, but the delay has to become longer ?

Hi Marvinh,

Is this a question regarding some code, or are you perhaps asking about the window length of higher-order polyblep? (I see your link containing some functions S0 through S6.) Could you please be a bit more specific?

Best,
jokona

[Marvinh]

Oh it’s about this step in general

“Build C1 continuous impulse from two cubic "smooth steps" which on their own looks like p(t)=3*t^2-2*t^3”

In reference to higher order bleeps I was interested in how the formulas were made so I looked up “continuous steps” and saw there were higher orders of the original function in the first post by mystran .

They use S instead C in the higher order formulas

Question being do we need longer delays for higher orders ?