// 4th order Linkwitz-Riley filters
//
// The filter is unstable for fast automation changes in the lower frequency range.
// Parameter interpolation and/or oversampling should fix this.
//
// The sum of the Linkwitz-Riley (Butterworth squared) HP and LP outputs, will result an all-pass filter
// at Fc and flat magnitude response - close to ideal for crossovers.
//
// Routing:
//
//  input-->lowpassA-->lowpassB/highpassB summed (allpass)-->low frequency output
//  input-->highpassA-->lowpassC-->mid frequency output
//  input-->highpassA-->highpassC-->high frequency output
//
// Input ch 3-6 disabled.

desc: LO=1+2  MID=3+4  HI=5+6

slider1:44.58<15.64,89.94,0.01>Frequency (Scale)
slider2:84.76<15.64,89.94,0.01>Frequency (Scale)

in_pin:L in
in_pin:R in
in_pin:xxx
in_pin:xxx
in_pin:xxx
in_pin:xxx
out_pin:L out 1
out_pin:R out 2
out_pin:L out 3
out_pin:R out 4
out_pin:L out 5
out_pin:R out 6

@init
sqrt2=sqrt(2);

@slider
// frequency slider scaling
tmpx = 16+slider1*1.20103;
tmpy = floor(exp(tmpx*log(1.059))*8.17742);

fc=tmpy;

wc=2*$pi*fc;
wc2=wc*wc;
wc3=wc2*wc;
wc4=wc2*wc2;
k=wc/tan($pi*fc/srate);
k2=k*k;
k3=k2*k;
k4=k2*k2;
sq_tmp1=sqrt2*wc3*k;
sq_tmp2=sqrt2*wc*k3;
a_tmp=4*wc2*k2+2*sq_tmp1+k4+2*sq_tmp2+wc4;

b1=(4*(wc4+sq_tmp1-k4-sq_tmp2))/a_tmp;
b2=(6*wc4-8*wc2*k2+6*k4)/a_tmp;
b3=(4*(wc4-sq_tmp1+sq_tmp2-k4))/a_tmp;
b4=(k4-2*sq_tmp1+wc4-2*sq_tmp2+4*wc2*k2)/a_tmp;

// low-pass
a0=wc4/a_tmp;
a1=4*wc4/a_tmp;
a2=6*wc4/a_tmp;
a3=a1;
a4=a0;

// high-pass
a0=k4/a_tmp;
a1=-4*k4/a_tmp;
a2=6*k4/a_tmp;
a3=a1;
a4=a0;


// frequency slider scaling
tmpxB = 16+slider2*1.20103;
tmpyB = floor(exp(tmpxB*log(1.059))*8.17742);

fcB=tmpyB;

wcB=2*$pi*fcB;
wc2B=wcB*wcB;
wc3B=wc2B*wcB;
wc4B=wc2B*wc2B;
kB=wcB/tan($pi*fcB/srate);
k2B=kB*kB;
k3B=k2B*kB;
k4B=k2B*k2B;
sq_tmp1B=sqrt2*wc3B*kB;
sq_tmp2B=sqrt2*wcB*k3B;
a_tmpB=4*wc2B*k2B+2*sq_tmp1B+k4B+2*sq_tmp2B+wc4B;

b1B=(4*(wc4B+sq_tmp1B-k4B-sq_tmp2B))/a_tmpB;
b2B=(6*wc4B-8*wc2B*k2B+6*k4B)/a_tmpB;
b3B=(4*(wc4B-sq_tmp1B+sq_tmp2B-k4B))/a_tmpB;
b4B=(k4B-2*sq_tmp1B+wc4B-2*sq_tmp2B+4*wc2B*k2B)/a_tmpB;

// low-pass
a0B=wc4B/a_tmpB;
a1B=4*wc4B/a_tmpB;
a2B=6*wc4B/a_tmpB;
a3B=a1B;
a4B=a0B;

// high-pass
a0B=k4B/a_tmpB;
a1B=-4*k4B/a_tmpB;
a2B=6*k4B/a_tmpB;
a3B=a1B;
a4B=a0B;

@sample
spl2=spl3=spl4=spl5=0;

tempx0=spl0;
tempx1=spl1;

tempy0=a0*tempx0+a1*xm1A+a2*xm2A+a3*xm3A+a4*xm4A-b1*ym1A-b2*ym2A-b3*ym3A-b4*ym4A;
xm4A=xm3A;
xm3A=xm2A;
xm2A=xm1A;
xm1A=tempx0;
ym4A=ym3A;
ym3A=ym2A;
ym2A=ym1A;
ym1A=tempy0;

tempy1=a0*tempx1+a1*xm1B+a2*xm2B+a3*xm3B+a4*xm4B-b1*ym1B-b2*ym2B-b3*ym3B-b4*ym4B;
xm4B=xm3B;
xm3B=xm2B;
xm2B=xm1B;
xm1B=tempx1;
ym4B=ym3B;
ym3B=ym2B;
ym2B=ym1B;
ym1B=tempy1;

out0=spl0-tempy0; // lo
out1=spl1-tempy1;
out2=tempy0; // hi
out3=tempy1;


tempy01=a0B*out0+a1B*xm1AB+a2B*xm2AB+a3B*xm3AB+a4B*xm4AB-b1B*ym1AB-b2B*ym2AB-b3B*ym3AB-b4B*ym4AB;
xm4AB=xm3AB;
xm3AB=xm2AB;
xm2AB=xm1AB;
xm1AB=out0;
ym4AB=ym3AB;
ym3AB=ym2AB;
ym2AB=ym1AB;
ym1AB=tempy01;

tempy11=a0B*out1+a1B*xm1BB+a2B*xm2BB+a3B*xm3BB+a4B*xm4BB-b1B*ym1BB-b2B*ym2BB-b3B*ym3BB-b4B*ym4BB;
xm4BB=xm3BB;
xm3BB=xm2BB;
xm2BB=xm1BB;
xm1BB=out1;
ym4BB=ym3BB;
ym3BB=ym2BB;
ym2BB=ym1BB;
ym1BB=tempy11;

out4=out0-tempy01; // lo
out5=out1-tempy11;
out6=tempy01; // hi
out7=tempy11;

spl0=out4+out6;
spl1=out5+out7;


tempy02=a0B*out2+a1B*xm1AC+a2B*xm2AC+a3B*xm3AC+a4B*xm4AC-b1B*ym1AC-b2B*ym2AC-b3B*ym3AC-b4B*ym4AC;
xm4AC=xm3AC;
xm3AC=xm2AC;
xm2AC=xm1AC;
xm1AC=out2;
ym4AC=ym3AC;
ym3AC=ym2AC;
ym2AC=ym1AC;
ym1AC=tempy02;

tempy12=a0B*out3+a1B*xm1BC+a2B*xm2BC+a3B*xm3BC+a4B*xm4BC-b1B*ym1BC-b2B*ym2BC-b3B*ym3BC-b4B*ym4BC;
xm4BC=xm3BC;
xm3BC=xm2BC;
xm2BC=xm1BC;
xm1BC=out3;
ym4BC=ym3BC;
ym3BC=ym2BC;
ym2BC=ym1BC;
ym1BC=tempy12;

spl2=out2-tempy02; // lo
spl3=out3-tempy12;
spl4=tempy02; // hi
spl5=tempy12;


spl0=min(max(spl0,-1),1);
spl1=min(max(spl1,-1),1);
spl2=min(max(spl2,-1),1);
spl3=min(max(spl3,-1),1);
spl4=min(max(spl4,-1),1);
spl5=min(max(spl5,-1),1);

@gfx 0 46
gfx_x=10;gfx_y=10;
gfx_r=gfx_b=0;
gfx_g=gfx_a=1;

gfx_drawchar($'F');
gfx_drawchar($' ');
gfx_drawchar($'=');
gfx_drawchar($' ');
gfx_drawnumber(tmpy,0);
gfx_drawchar($' ');
gfx_drawchar($'H');
gfx_drawchar($'z');

gfx_x=10;gfx_y=30;
gfx_drawchar($'F');
gfx_drawchar($' ');
gfx_drawchar($'=');
gfx_drawchar($' ');
gfx_drawnumber(tmpyB,0);
gfx_drawchar($' ');
gfx_drawchar($'H');
gfx_drawchar($'z');
