|
|
|
09-15-2015, 03:31 PM
|
#1
|
Human being with feelings
Join Date: Jul 2014
Location: Là bas les huîtres (FR)
Posts: 424
|
JSFX Multimode State Variable Filter
I was digging some of my old c+- code from my phd era... and i found this filter that can be of interest if one needs an easy to use and versatile filter.
It is a multimode State Variable Filter with frequency and resonance controls and simultaneous low-pass, band-pass and high-pass outputs.
As far as i remember, i designed it mainly to emulate the analog SEM filter.
Unlike the "Chamberlin" filter, it is stable and accurate. It even uses sort of "zero delay feedback"... well... it is rather a 1/32th of sample delay feedback.
I just finished to translate it to jsfx. I think that it can still be perfected in many ways.
The @init section contains all the functions (control rate and audio rate processings)
The @slider and @sample sections show how to use them (which is quite easy).
Code:
desc:TiaR state variable filter
/*
state variable filter T.Rochebois 1997
Integrateur Euler
Itérations matricielles
Factorisation/Réduction
feedback intermédiaire TSA/AXIS/IEF/Paris Sud Orsay
*/
slider1:sl_Fc=400<10,20000>Fc
slider2:sl_Q=1<0.02,1>1/Q
slider3:sl_mode=<0,2,1{LP,BP,HP}>Mode
// ___________________________________________________________________
@init
TRF_coef = 2 *$pi/(32 * srate);
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
/*
Itération matricielle réduite
*/
// S Source D Destination
// 4 x 5+
function TRF_iter(S* D*)local(b2)(
b2 = S.b * S.b;
D.a = b2 + S.a * (2 - S.a);
D.b = S.b *(1 + S.c - S.a);
D.c = S.c * S.c - b2;
);
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// Kontrol Rate Processing
// 25x 25+
function TRF_kProc(Fc _Q)
instance(F D a b c )
local(tmp b2 na nb nc) (
D = _Q;
F = Fc * TRF_coef;
a = F * F;
tmp = 1 - a - D * F;
b = F * tmp + F;
c = tmp * tmp - a; // 5 x 4 +
TRF_iter(this, this.A); // 4 x 5+
TRF_iter(this.A, this); // 4 x 5+
TRF_iter(this, this.A); // 4 x 5+
TRF_iter(this.A, this); // 4 x 5+
TRF_iter(this, this.A); // 4 x 5+
);
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// Audio Rate Processing
// 5 x 6 +
function TRF_aProc(x)
instance(lp bp hp a b c)
local( x_lp)(
x_lp = x - lp;
lp += a * x_lp + b * bp;
hp = x - D * bp - lp;
bp = b * x_lp + c * bp;
lp;
);
// ___________________________________________________________________
@slider
fl.TRF_kProc(sl_Fc, sl_Q); // fl : left filter
fr.TRF_kProc(sl_Fc, sl_Q); // fr : right filter
// ___________________________________________________________________
@sample
fl.TRF_aProc(spl0);
fr.TRF_aProc(spl1);
sl_mode === 0 ? (
spl0 = fl.lp; //low pass
spl1 = fr.lp;
) : sl_mode === 1 ? (
spl0 = fl.bp; //band pass
spl1 = fr.bp;
) : (
spl0 = fl.hp; //high pass
spl1 = fr.hp;
);
Note:
As i was googling some info on "zero" delay feedback filters...
It is not a TPT (topology preserving transform). It does not preserve topology as it adds some intermediary feedback.
It is more like http://levien.com/ladder.pdf applied to the Variable State Filter (instead of the moog), with many optimisations in the parameter matrix iterations.
Last edited by Smashed Transistors; 09-15-2015 at 04:14 PM.
|
|
|
09-16-2015, 12:32 PM
|
#2
|
Human being with feelings
Join Date: Jan 2007
Location: mcr:uk
Posts: 3,891
|
Thanks! Another lesson to learn.
|
|
|
09-16-2015, 02:15 PM
|
#3
|
Human being with feelings
Join Date: Jul 2014
Location: Là bas les huîtres (FR)
Posts: 424
|
I've found a few notes i wrote about this, and it is clearer for me now...
I reworked the code a little (exp cutoff control, interpolators).
I added a ramp LFO that modulates the cutoff and Q... basically it was to test the filters' stability... but produces a nice old school filter effect.
Code:
desc:TiaR state variable filter
/*
state variable filter T.Rochebois 1997
Integrateur Euler
Iterations matricielles
Factorisation/Reduction
feedback intermediaire TSA/AXIS/IEF/Paris Sud Orsay
Added
cutoff control a la MIDI (ie pitch = 69 <-> 440Hz)
linear interpolation of coefficients
LFO modulation for demo and test
*/
slider1:sl_Pitch1=110<0,150>Pitch1
slider2:sl_Pitch0=12<0,150>Pitch0
slider3:sl_Q1=1<0.01,1>1/Q
slider4:sl_Q0=0.11<0.01,1>1/Q
slider5:sl_rate=4<0,20>Rate
slider10:sl_mode=<0,2,1{LP,BP,HP}>Mode
// ___________________________________________________________________
// Library P A R T
// ___________________________________________________________________
@init
function TRF_softSat(x)( x*=2/3; x = min(1,max(-1,x)); x*(1.5 -0.5*x*x); );
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// Global init for all instances
function TRF()
local(i Fc)(
TRF_coef = 2 *$pi/(64 * srate);
KRATE === 0 ? (KRATE = 16; _KRATE = 1 / KRATE);
TRF_P2F = ad; ad += 150;
TRF_P2FD = ad; ad += 150;
i = 0;
loop(150,
Fc = 440 * 2 ^((i-69)*(1/12));
TRF_P2F[i] = 2*sin(0.5*Fc * TRF_coef);
Fc = 440 * 2 ^((i+1-69)*(1/12));
TRF_P2FD[i] = 2*sin(0.5*Fc * TRF_coef) - TRF_P2F[i];
i += 1;
);
);
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// S Source D Destination
// 4 x 5+
function TRF_iter(S* D*)local(b2)(
b2 = S.b * S.b;
D.a = b2 + S.a * (2 - S.a);
D.b = S.b * (1 + S.c - S.a);
D.c = S.c * S.c - b2;
);
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// Kontrol Rate Processing
// 33 x 25 +
function TRF_kProc(pitch _Q)
instance(F D a b c da db dc)
local(tmp b2 na nb nc p0) (
D = _Q;
p0 = pitch|0;
F = TRF_P2F[p0] + (pitch - p0) * TRF_P2FD[p0];
this.N.a = F * F;
tmp = 1 - this.N.a - D * F;
this.N.b = F * tmp + F;
this.N.c = tmp * tmp - this.N.a;
TRF_iter(this.N, this.A); // 4 x 5+
TRF_iter(this.A, this.N); // 4 x 5+
TRF_iter(this.N, this.A); // 4 x 5+
TRF_iter(this.A, this.N); // 4 x 5+
TRF_iter(this.N, this.A); // 4 x 5+
da = (this.A.a - a)*_KRATE; // prepare linear interpolators
db = (this.A.b - b)*_KRATE;
dc = (this.A.c - c)*_KRATE;
);
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// Audio Rate Processing
// 5 x 9 +
function TRF_aProc(x)
instance(lp bp hp a b c da db dc)
local( x_lp)(
a += da; b += db; c += dc; // linear interpolators
x_lp = x - lp;
lp += a * x_lp + b * bp;
hp = x - D * bp - lp;
bp = b * x_lp + c * bp;
lp;
);
// ___________________________________________________________________
// Usage P A R T
// ___________________________________________________________________
//init
TRF();
// ___________________________________________________________________
@slider
dp = sl_rate * KRATE/srate;
// ___________________________________________________________________
@sample
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// Kontrol Rate Processing
k >= KRATE ? ( // KRATE = 16 this section is ran evey 16 samples
p -= dp; // saw tooth LFO
p < 0 ? p += 1;
// LFO modulated pitch and Q
pitch = sl_Pitch0 + p * (sl_Pitch1 - sl_Pitch0);
_Q = sl_Q0 + p * (sl_Q1 - sl_Q0 );
// control rate processing for left and right filter
fl.TRF_kProc(pitch, _Q);
fr.TRF_kProc(pitch, _Q);
k = 0;
);
k += 1;
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// Audio Rate Processing
fl.TRF_aProc(spl0);
fr.TRF_aProc(spl1);
sl_mode === 0 ? ( spl0 = fl.lp; spl1 = fr.lp;)
: sl_mode === 1 ? ( spl0 = fl.bp; spl1 = fr.bp;)
: ( spl0 = fl.hp; spl1 = fr.hp;);
//soft clipping
spl0 = TRF_softSat(spl0); spl1 = TRF_softSat(spl1);
Last edited by Smashed Transistors; 09-16-2015 at 02:27 PM.
|
|
|
09-16-2015, 02:43 PM
|
#4
|
Human being with feelings
Join Date: Oct 2013
Location: Seattle, WA
Posts: 876
|
Thanks! ...now to figure out how it works.
*adding to checklist - bone up on matrix math*
|
|
|
09-16-2015, 03:10 PM
|
#5
|
Human being with feelings
Join Date: Oct 2013
Location: Seattle, WA
Posts: 876
|
Something I wanted to throw out there - the classic cubic clipper is an approximation of tanh(). Approximating sin() is a little smoother and generates less higher harmonics. Might be worth considering.
This is one of the ways I write it :
Code:
function cubic6(x) (
(abs(x)>1.4142135623) ? sign(x)*0.9428090416 : (x-(x*x*x)/6);
);
Of course, if heavily driven it sounds best with oversampling...
|
|
|
09-17-2015, 02:06 AM
|
#6
|
Human being with feelings
Join Date: Jul 2008
Location: The Netherlands
Posts: 3,653
|
Quote:
Originally Posted by Smashed Transistors
basically it was to test the filters' stability...
|
I have only quickly tested (using REAPER's paramater modulation), but the filter seems stable enough. And it seems rather light on CPU. So thanks for sharing.
|
|
|
09-18-2015, 03:17 AM
|
#7
|
Human being with feelings
Join Date: Jul 2014
Location: Là bas les huîtres (FR)
Posts: 424
|
Quote:
Originally Posted by SaulT
Thanks! ...now to figure out how it works.
*adding to checklist - bone up on matrix math*
|
I did that too... refreshing memories, almost forgotten formulas.
I'll eventually take some time to write some explanatory doc... including for myself.
I like the classic cubic clipper because it has a zero derivative at saturation point ... and also because i can remember it.
|
|
|
09-18-2015, 03:23 AM
|
#8
|
Human being with feelings
Join Date: Jul 2014
Location: Là bas les huîtres (FR)
Posts: 424
|
Quote:
Originally Posted by Tale
I have only quickly tested (using REAPER's paramater modulation), but the filter seems stable enough. And it seems rather light on CPU. So thanks for sharing.
|
Thanks Tale, and thanks for sharing your implementation of the TPT SVF
|
|
|
09-19-2015, 12:32 AM
|
#9
|
Human being with feelings
Join Date: Jul 2008
Location: The Netherlands
Posts: 3,653
|
Quote:
Originally Posted by Smashed Transistors
Thanks Tale, and thanks for sharing your implementation of the TPT SVF
|
You're welcome. BTW, it is actually Robin Schmidt's implementation, but I did port it to JSFX.
|
|
|
09-21-2015, 04:30 PM
|
#10
|
Human being with feelings
Join Date: Oct 2013
Location: Seattle, WA
Posts: 876
|
Is there any chance you could explain what is happening in the TRF_iter() portion of the code? The paper wasn't nearly as much of a help as I was hoping. It seems like this should be the area where the matrix is exponentiated, but I'm not following exactly how that's happening. Any explanation would be greatly appreciated!
|
|
|
09-22-2015, 11:24 AM
|
#11
|
Human being with feelings
Join Date: Jul 2014
Location: Là bas les huîtres (FR)
Posts: 424
|
I've had a look at my old notes and figured out some stuff.
The starting point is the Chamberlin filter:
Code:
For every new input sample x, the state variables are updated:
n n-1
lp <- lp + F * bp
bp <- F * x - F * lp + (1 -F*D - F*F) * bp
F being related to the normalized frequency and D being a resonance factor (1/2Q).
Note: the Chamberlin filter is sort of "zero delay feedback" filter as it uses a backward and a forward Euler integrator. The drawback is that it is not stable for all F and D.
The Chamberlin filter can be decribed as a transition matrix M:
Code:
x lp bp
x 1 0 0 constant input
lp 0 1 F bp feeds the lp integrator
bp F -F 1 - F*D - F*F double feedback on bp
If you use M^k, instead of M^1, it's just like if you iterated the filter k times... its just like if you upsampled the filter k times with a constant input.
When you exponentiate M to k=2^i by hand (or with Maxima) you quickly see that the form of matrix M^k remains
Code:
x lp bp
x 1 0 0
lp a 1-a b
bp b -b c
an iteration of the matrix exponentiation can be reduced to:
Code:
i+1 i
a <- b^2 + a * (2 - a)
b <- b * (1 + c - a)
c <- c^2 - b^2
So it's easy to get M^2^i.
With only 6 iterations on coefficients a,b and c, you get M^64.
If you use M^64 instead of M you've got the equivalent of a 64x
upsampling (with a constant input). You get better bandpass and lowpass
filter response and the filter is stable (which is not the case
of the Chamberlin filter).
The audio rate processing is still quite simple:
Code:
n n-1
lp <- a * x + (1 - a) * lp + b * bp;
bp <- b * x - b * lp + c * bp;
|
|
|
09-22-2015, 11:33 AM
|
#12
|
Human being with feelings
Join Date: Jul 2014
Location: Là bas les huîtres (FR)
Posts: 424
|
Some more experiments
In many papers and books i've read about the "Step Invariant Transform" and the matrix exponential exp(M).
Starting with M the "analog" matrix expresses the relationships between the variables and their derivatives.
Code:
x lp bp
x' 0 0 0 constant "step" input
bp' F -F*D -F double feedback on bp
lp' 0 F 0 bp feeds the lp integrator
We get the digital transition matrix as
Code:
exp(M) = I + M + M^2/2! + M^3/3!...
Numerically it converges even faster than the M^2^i, and the calc iterations can be optimised.
But, i experience unstable behaviour when i use less than 8-9 terms.
So even if it's theoretically better, i'd prefer the previous method... or maybe i got something wrong.
__________________________________________________ _____________
I've just tried to use a "Ramp invariant" version of the method: instead of a constant input, i consider that the input varies linearly. It resulted in very slight improvements of the high pass filter.
__________________________________________________ _____________
Last and for fun, i edited the code to audio modulate the filters:
Change line 20
Code:
slider5:sl_rate=4<0,10000>Rate
And line 99 add
Code:
KRATE = 1;_KRATE = 1/KRATE;
just BEFORE
Vintage dissonant effects à la Beaubourg (Vangelis)...
Last edited by Smashed Transistors; 09-22-2015 at 02:49 PM.
|
|
|
09-23-2015, 03:42 PM
|
#13
|
Human being with feelings
Join Date: Oct 2013
Location: Seattle, WA
Posts: 876
|
Thank you for the explanation, it really helps. It's making more sense to me now, the question for me is whether I can take this concept and run any further with it. E.g. if there's a way to apply saturation within the algorithm (is "analytically" the correct term?) rather than just tacking it on after the fact, e.g. similar to the tanh stages in a Moog.
|
|
|
09-23-2015, 07:53 PM
|
#15
|
Human being with feelings
Join Date: Jul 2012
Location: inside a man
Posts: 84
|
holy crap batman. this thing is tits! i'm sure it must be possible to scale the rate logarithmicly so that it would align to pitch values? even still, this is one of the coolest digital filters ive messed with.
|
|
|
09-23-2015, 08:22 PM
|
#16
|
Human being with feelings
Join Date: Jul 2012
Location: inside a man
Posts: 84
|
ALSO: have you tried setting the rate to 0 and messing with the sliders? i cant quite figure out what its doing but i really like it.
|
|
|
09-24-2015, 01:53 AM
|
#17
|
Human being with feelings
Join Date: Jul 2012
Location: inside a man
Posts: 84
|
i figured it out, and made a little utility that follows midi pitch and can control rate and filter cutoff by linking the sliders with parameter modulation. i figured it was more flexible that way, plus this is my very first js effect or even attempt at scripting
Code:
// author: kamalmanzukie
// date: 24.9.2015
desc: midi/param mod utility for TiaR's gr8 filter
slider1:70<0,150,1>slidername
slider2:0.5<0,10000,0.01>bob
@init
// -- note_to_freq
_md1_d12=1/12;
@slider
_md0_value = slider1;
_md3_value = slider2;
@sample
@block
// -- note_to_freq
while (midirecv(_md1_offset,_md1_msg1,_md1_msg23) ?
(
_md1_status=_md1_msg1&240;
(_md1_status== 8*16||_md1_status==9*16) ?
(
_md1_note=_md1_msg23&127;
_md1_frequency=440*2^((_md1_note-69)*_md1_d12);
):(
midisend(_md1_offset,_md1_msg1,_md1_msg23);
);
);
);
_md2_midi = _md1_frequency;
// -- midi_to_freq
_md2_frequency = (exp(0.9172010659072+_md2_midi*0.0542122581767108208)*8.17742)|0;
slider2=_md1_frequency;
sliderchange(slider2);
slider1=_md1_note;
sliderchange(slider1)
@gfx
@serialize
|
|
|
09-24-2015, 01:18 PM
|
#18
|
Human being with feelings
Join Date: Jul 2014
Location: Là bas les huîtres (FR)
Posts: 424
|
Quote:
Originally Posted by SaulT
Thank you for the explanation, it really helps. It's making more sense to me now, the question for me is whether I can take this concept and run any further with it. E.g. if there's a way to apply saturation within the algorithm (is "analytically" the correct term?) rather than just tacking it on after the fact, e.g. similar to the tanh stages in a Moog.
|
Hi SaulT,
here is my first attempt to do such a thing.
The non linearities are evaluated once per sample. I do not use them to modify the internal variables, i use them to change the integrators behavior
(non linear input gains). The aliasing seem low to me (as the saturated signals are integrated they are transformed into triangle waves... just like in the analog filters).
See the comments for details.
Maybe other input functions can be more interesting...
Code:
desc:Ze Nasty Little Filter 02
slider1:sl_pitch1=120<0,150>Pitch1
slider2:sl_pitch0=60<0,150>Pitch0
slider3:sl_D1=1<0.001,1>D1
slider4:sl_D0=1<0.001,1>D0
slider5:sl_rate=4<0,20>rate
slider6:sl_mode=0<0,2,1{LP,BP,HP}>Mode
slider7:sl_drive=0<-12,36>Drive (dB)
slider8:sl_outGain=0<-36,12>Output Gain (dB)
slider9:sl_satType=0<0,3,1{None,Sat3,Tanh,Sat1}>Sat Type
/*
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
Template Filter Structure
I use a Chamberlin filter as the filter template.
It is unstable under certain conditions but it offers good resonance
and frequency behavior (sort of "zero delay feedback" thanks to a
combo of backward and forward Euler integrators).
Matrix exponentiations are equivalent to step invariant oversampling.
They solve the stability issues and enforce the "ZDF" properties of
the filter.
-1_____________________________
/___________________ \ Filter
-D / \ \ feedbacks
v HP BP \ \ LP
--> + --->F1-- + -->[Z-1]----->F2-- + --->[Z-1]----->
^ / ^ /
\_________/ \__________/ Integrator
feedbacks
F1: BP integrator gain, F2: LP integrator gain, D = 1/2Q
for a Chamberlin filter F1 = F2 = F = 2sin(pi*f/2) f:normalized freq
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
Non linearities
I model non linearities at integrator inputs as non linear gains.
I modulate the integrator gains by non linear functions that depend on
their inputs:
F1 = F * g(HP)
F2 = F * g(BP)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
Matrix iterations
The matrix is exponentiated to obtain "step invariant oversampled"
matrixes.
Note: The iterations change the structure of the filter. Note that the
LP integrator is fed with the input signal. This is not a topology
preserving method.
Matrix for Chamberlin SVF Iterated Matrix Matrix iteration
i.e. Matrix0 M <- M^2
x LP BP x LP BP a <- a^2 - b*d;
x 1 0 0 x 1 0 0 b <- b * (a+c);
LP 0 1 F2 LP (1-a) a d c <- c^2 - b*d;
BP F1 -F1 1-DF1-F1F2 BP b -b c d <- d * (a+c);
|
V
LP = (1-a) * x + a * LP + d * BP
BP = b * x - b * LP + c * BP
*/
//____________________________________________________________________
@init
ZNF_iter = 2;
ZNF_coef = (2^(-ZNF_iter-69/12))*440*$pi/srate;
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// Non linearities seen as gains
// g = cubic soft sat divided by x --> bell shaped
// todo: tabulate
function tanh(x) ( x = exp(2*x); (x - 1) / (x + 1) );
function ZNF_gSat1(x)( abs(x) > 1 ? 1 / abs(x) : 1; );
function ZNF_gSat3(x)( abs(x) > 1.5 ? 1 / abs(x) : 1 - x * x * (4/27); );
function ZNF_gTanh(x)( x == 0 ? 1 : tanh(x) / x; );
function ZNF_g(x)(
sl_satType === 0 ? 1
: sl_satType === 1 ? ZNF_gSat3(x)
: sl_satType === 2 ? ZNF_gTanh(x)
: sl_satType === 3 ? ZNF_gSat1(x)
;
);
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// audio rate process
function ZNF_proc(x pitch _2Q)
local(ac bd x_LP)
instance(LP BP HP F F1 F2 a b c d)(
// once tabulated --> 3 adds 3 muls
F = 2 * sin(min($pi/2,ZNF_coef*2^(pitch*(1/12)))); // todo: tabulate
F1 = F * ZNF_g(HP); // gain for BP integrator (its input is HP)
F2 = F * ZNF_g(BP); // gain for LP integrator (its input is BP)
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// Matrix init and iterations
// 2 + 3i adds 1 + 5i muls
a = 1;
b = F1;
c = 1 - F1 * (_2Q + F2);
d = F2;
loop(ZNF_iter,
ac = a + c;
bd = b * d;
a = a*a - bd;
b = b * ac;
c = c*c - bd;
d = d * ac;
);
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// Filter calc
// 6 adds 5 muls
x_LP = x - LP;
LP = x - a * x_LP + d * BP;
HP = x - _2Q * BP - LP;
BP = b * x_LP + c * BP;
);
//____________________________________________________________________
function sat3(x)local(x2)(
x = min(215/104, max(-215/104,x));
x2 = x*x;
x * (1
+ x2 * ( ( -411008.0/1987675.0)
+ x2 * ( ( 11581599744.0/459401384375.0)
+ x2 *( -5061276073984.0/4247165798546875.0))));
);//____________________________________________________________________
@slider
dp = -sl_rate/srate; // positive rate for down ramp
gain = 2^(sl_drive/6); // input gain
_gain = 2^(sl_outGain/6); // output gain
//____________________________________________________________________
@sample
p += dp; p -= (p >= 1); p += (p < 0); // lfo phase inc and modulo
pitch = sl_pitch0 + p * (sl_pitch1 - sl_pitch0);
_2Q = sl_D0 + p * ( sl_D1 - sl_D0);
l.ZNF_proc(gain * spl0, pitch,_2Q);
r.ZNF_proc(gain * spl1, pitch,_2Q);
sl_mode === 0 ? (spl0 = l.LP; spl1 = r.LP; )
: sl_mode === 1 ? (spl0 = l.BP; spl1 = r.BP; )
: sl_mode === 2 ? (spl0 = l.HP; spl1 = r.HP; );
spl0 *= _gain;
spl1 *= _gain;
spl0 = sat3(spl0);
spl1 = sat3(spl1);
Last edited by Smashed Transistors; 09-24-2015 at 04:18 PM.
|
|
|
09-24-2015, 01:20 PM
|
#19
|
Human being with feelings
Join Date: Jul 2014
Location: Là bas les huîtres (FR)
Posts: 424
|
Quote:
Originally Posted by kamalmanzkie
i figured it out, and made a little utility that follows midi pitch and can control rate and filter cutoff by linking the sliders with parameter modulation. i figured it was more flexible that way, plus this is my very first js effect or even attempt at scripting
|
Thanks for the trick Kamalmanzukie !
|
|
|
09-25-2015, 04:40 PM
|
#20
|
Human being with feelings
Join Date: Jul 2012
Location: inside a man
Posts: 84
|
its cool if you have two controllers, one for keytracking filter, one for notes. i do this with my volca keys, except my volca cant modulate at insane audio rate, only pretty good audio rate. but yeah, playing filters like they were notes, all for that. if i use this on my volca at similar rates it sounds near identical, too.
|
|
|
09-25-2015, 05:58 PM
|
#21
|
Human being with feelings
Join Date: Oct 2013
Location: Seattle, WA
Posts: 876
|
I disengaged all the modulation to hear the sound of the filter itself, and I've got to say - this filter sounds filthy! Well done!
Big thanks for the explanations, I think I've finally got my head around it. After a little bit of work I've gone back and implemented a Moog filter using this same technique. There's work that needs to be done, for sure, and it doesn't have the filthy goodness that this does... but there's time. I'll publish it when I've done some polishing on it.
I like where you're going with the nonlinearities, I'll try a few different ones myself and see what I can come up with.
|
|
|
09-26-2015, 12:13 AM
|
#22
|
Human being with feelings
Join Date: Oct 2013
Location: Seattle, WA
Posts: 876
|
I have an observation to make. As I tweaked your filter to increase the number of exponentiations, I noticed that the filter became thinner, as if the Fc had changed. I then went back to the paper that I have been working off of to develop a Moog matrix filter...
http://levien.com/ladder.pdf
What I noticed was that the matrix had "a delta t" and "1 - a delta t". Now, when I implemented "a" and "1-a" in a non-exponentiated matrix it performed as expected, but when once I started squaring the matrix it thinned out the same way that yours did!
After a bit of thought, I went back to what you said earlier, that the point of squaring a matrix was to simulate oversampling. Well, oversampling is expanding the bandwidth by a certain factor... and since Fc is normalized, that would move it as well.
So I did a bit of experimenting, and found that for the ladder filter multiplying a by 1/oversampling factor put the results right back to where they should be and the filter got its sound back. I did the same for your filter and it worked out as well! In other words, dt is 1/oversampling.
If you want I can post the whole code, or I can just show the highlights...
Code:
slider10:ZNF_iter=2<1,8,1>Iterations
Code:
@init
correction = 1/ZNF_iter;
Code:
@slider
correction = 1/ZNF_iter;
Code:
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// Matrix init and iterations
// 2 + 3i adds 1 + 9i muls
a = 1;
b = F1*correction;
c = 1 - F1 * correction * (_2Q + F2 * correction);
d = F2 * correction;
loop(ZNF_iter,
ac = a + c;
bd = b * d;
a = a*a - bd;
b = b * ac;
c = c*c - bd;
d = d * ac;
);
...anyways, that seems to work. I won't have time to work on it for a few days, I think, but I wanted to post it where I was at to see what you thought.
|
|
|
09-27-2015, 03:56 AM
|
#23
|
Human being with feelings
Join Date: Jul 2014
Location: Là bas les huîtres (FR)
Posts: 424
|
Hi, SaulT
Quote:
Originally Posted by SaulT
I have an observation to make. As I tweaked your filter to increase the number of exponentiations, I noticed that the filter became thinner, as if the Fc had changed.
...
After a bit of thought, I went back to what you said earlier, that the point of squaring a matrix was to simulate oversampling. Well, oversampling is expanding the bandwidth by a certain factor... and since Fc is normalized, that would move it as well.
...
|
That's right. In fact, it the "oversampling factor" is 2^nb_iterations. In the latest code, i the correction is hidden in the pitch to Fc translation :
Code:
ZNF_coef = (2^(-ZNF_iter-69/12))*440*$pi/srate;
. One octave per iteration.
The state variable filter uses integrators while the Moog filter model is a chain of low pass filters (hence the 1-a factor) with feedback.
I think that the resonance of the Moog filter is quite tricky to simulate as the fed back signal path is quite long and the simulation needs accurate phase response.
I've found a thread on KVR http://www.kvraudio.com/forum/viewtopic.php?t=385262 that discusses the Moog filter exponentiation. They seem to have found an analytic solution to the Exp(M) method. That may be helpful if you want to make a zero delay fedback Moog Filter Hope you'll get some time to do it.
|
|
|
09-27-2015, 10:31 AM
|
#24
|
Human being with feelings
Join Date: Apr 2007
Posts: 372
|
Wow this filter sounds very very good !!!
Is it possible to have it non modulated ?
(although an midi triggered envelope would be amazing!
Thanks for you work on this !
Quote:
Originally Posted by Smashed Transistors
I was digging some of my old c+- code from my phd era... and i found this filter that can be of interest if one needs an easy to use and versatile filter.
It is a multimode State Variable Filter with frequency and resonance controls and simultaneous low-pass, band-pass and high-pass outputs.
As far as i remember, i designed it mainly to emulate the analog SEM filter.
Unlike the "Chamberlin" filter, it is stable and accurate. It even uses sort of "zero delay feedback"... well... it is rather a 1/32th of sample delay feedback.
I just finished to translate it to jsfx. I think that it can still be perfected in many ways.
The @init section contains all the functions (control rate and audio rate processings)
The @slider and @sample sections show how to use them (which is quite easy).
Code:
desc:TiaR state variable filter
/*
state variable filter T.Rochebois 1997
Integrateur Euler
Itérations matricielles
Factorisation/Réduction
feedback intermédiaire TSA/AXIS/IEF/Paris Sud Orsay
*/
slider1:sl_Fc=400<10,20000>Fc
slider2:sl_Q=1<0.02,1>1/Q
slider3:sl_mode=<0,2,1{LP,BP,HP}>Mode
// ___________________________________________________________________
@init
TRF_coef = 2 *$pi/(32 * srate);
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
/*
Itération matricielle réduite
*/
// S Source D Destination
// 4 x 5+
function TRF_iter(S* D*)local(b2)(
b2 = S.b * S.b;
D.a = b2 + S.a * (2 - S.a);
D.b = S.b *(1 + S.c - S.a);
D.c = S.c * S.c - b2;
);
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// Kontrol Rate Processing
// 25x 25+
function TRF_kProc(Fc _Q)
instance(F D a b c )
local(tmp b2 na nb nc) (
D = _Q;
F = Fc * TRF_coef;
a = F * F;
tmp = 1 - a - D * F;
b = F * tmp + F;
c = tmp * tmp - a; // 5 x 4 +
TRF_iter(this, this.A); // 4 x 5+
TRF_iter(this.A, this); // 4 x 5+
TRF_iter(this, this.A); // 4 x 5+
TRF_iter(this.A, this); // 4 x 5+
TRF_iter(this, this.A); // 4 x 5+
);
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// Audio Rate Processing
// 5 x 6 +
function TRF_aProc(x)
instance(lp bp hp a b c)
local( x_lp)(
x_lp = x - lp;
lp += a * x_lp + b * bp;
hp = x - D * bp - lp;
bp = b * x_lp + c * bp;
lp;
);
// ___________________________________________________________________
@slider
fl.TRF_kProc(sl_Fc, sl_Q); // fl : left filter
fr.TRF_kProc(sl_Fc, sl_Q); // fr : right filter
// ___________________________________________________________________
@sample
fl.TRF_aProc(spl0);
fr.TRF_aProc(spl1);
sl_mode === 0 ? (
spl0 = fl.lp; //low pass
spl1 = fr.lp;
) : sl_mode === 1 ? (
spl0 = fl.bp; //band pass
spl1 = fr.bp;
) : (
spl0 = fl.hp; //high pass
spl1 = fr.hp;
);
Note:
As i was googling some info on "zero" delay feedback filters...
It is not a TPT (topology preserving transform). It does not preserve topology as it adds some intermediary feedback.
It is more like http://levien.com/ladder.pdf applied to the Variable State Filter (instead of the moog), with many optimisations in the parameter matrix iterations.
|
Last edited by geoslake; 09-27-2015 at 10:43 AM.
|
|
|
09-27-2015, 07:48 PM
|
#25
|
Human being with feelings
Join Date: Oct 2013
Location: Seattle, WA
Posts: 876
|
Quote:
The state variable filter uses integrators while the Moog filter model is a chain of low pass filters (hence the 1-a factor) with feedback.
|
I'll go back to working on code, I don't mean to hijack your thread!
I've implemented a ZDF/TPT moog a la Zavalishin, but it breaks in anything but LP mode when oversampled, thus my interest in tackling it via this matrix method.
I will start a separate thread when I have some code I'm happy with.
Last edited by SaulT; 09-27-2015 at 09:45 PM.
|
|
|
09-28-2015, 01:08 PM
|
#26
|
Human being with feelings
Join Date: Jul 2014
Location: Là bas les huîtres (FR)
Posts: 424
|
@SaulT
I look forward to see where the Moog ZDF/TPT + matrix track will lead you.
Excuse my poor English, what do you mean by "it breaks in anything but LP mode when oversampled" ? I think that the Moog is essentially LP.
I've done some experimentations with the Moog filter eons ago, but i used a completely different trick to enhance the resonance. I'll dig it up.
@geoslake
i will certainly include it as an update of one of my synths
Meanwhile, here is another version, with a tempo synch'd LFO, with lots of sliders to control the LFO shape. I added an asymetrical saturation... quite nice for dirty sounds.
Code:
desc:Ze Nasty Little Filter 03
slider1:sl_pitch1=130<0,140>Pitch1
slider2:sl_pitch0=60<0,140>Pitch0
slider3:sl_D1=0.15<0.001,1>D1
slider4:sl_D0=1<0.001,1>D0
slider6:sl_mode=0<0,2,1{LP,BP,HP}>Mode
slider7:sl_drive=16<-12,36>Drive (dB)
slider8:sl_outGain=-13<-36,12>Output Gain (dB)
slider9:sl_satType=4<0,4,1{None,Sat3,Tanh,Sat1,Asymetrical}>Sat Type
slider11:sl_sw1=-0.5<-1,1>Swing 1
slider12:sl_sw2=-0.50<-1,1>Swing 2
slider13:sl_sw3=-0.50<-1,1>Swing 3
slider14:sl_sym=-0.9<-1,1>Triangle
slider15:sl_exp=0.2<-1,1>Exp
slider16:sl_sin=1<0,1>Sin
slider18:sl_rateScale=1<0.25,4>LFO rate scale
/*
T.Rochebois Sept 2015
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
Template Filter Structure
I use a Chamberlin filter as the filter template.
It is unstable under certain conditions but it offers good resonance
and frequency behavior (sort of "zero delay feedback" thanks to a
combo of backward and forward Euler integrators).
Matrix exponentiations are equivalent to step invariant oversampling.
They solve the stability issues and enforce the "ZDF" properties of
the filter.
-1_____________________________
/___________________ \ Filter
-D / \ \ feedbacks
v HP BP \ \ LP
--> + --->F1-- + -->[Z-1]----->F2-- + --->[Z-1]----->
^ / ^ /
\_________/ \__________/ Integrator
feedbacks
F1: BP integrator gain, F2: LP integrator gain, D = 1/2Q
for a Chamberlin filter F1 = F2 = F = 2sin(pi*f/2) f:normalized freq
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
Non linearities
I model non linearities at integrator inputs as non linear gains.
I modulate the integrator gains by non linear functions that depend on
their inputs:
F1 = F * g(HP)
F2 = F * g(BP)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
Matrix iterations
The matrix is exponentiated to obtain "step invariant oversampled"
matrixes.
Note: The iterations change the structure of the filter. Note that the
LP integrator is fed with the input signal. This is not a topology
preserving method.
Matrix for Chamberlin SVF Iterated Matrix Matrix iteration
i.e. Matrix0 M <- M^2
x LP BP x LP BP a <- a^2 - b*d;
x 1 0 0 x 1 0 0 b <- b * (a+c);
LP 0 1 F2 LP (1-a) a d c <- c^2 - b*d;
BP F1 -F1 1-DF1-F1F2 BP b -b c d <- d * (a+c);
|
V
LP = (1-a) * x + a * LP + d * BP
BP = b * x - b * LP + c * BP
*/
//____________________________________________________________________
@init
ZNF_iter = 2;
ZNF_coef = (2^(-ZNF_iter-69/12))*440*$pi/srate;
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// Non linearities seen as gains
// g = cubic soft sat divided by x --> bell shaped
// todo: tabulate
function tanh(x) ( x = exp(2*x); (x - 1) / (x + 1) );
function ZNF_gSat1(x)( abs(x) > 1 ? 1 / abs(x) : 1; );
function ZNF_gSat3(x)( abs(x) > 1.5 ? 1 / abs(x) : 1 - x * x * (4/27); );
function ZNF_gTanh(x)( x == 0 ? 1 : tanh(x) / x; );
function ZNF_gAsy(x)( x == 0 ? 1 : x>0 ? 1 : ZNF_gSat3(x); );
function ZNF_g(x)(
sl_satType === 0 ? 1
: sl_satType === 1 ? ZNF_gSat3(x)
: sl_satType === 2 ? ZNF_gTanh(x)
: sl_satType === 3 ? ZNF_gSat1(x)
: sl_satType === 4 ? ZNF_gAsy(x)
;
);
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// audio rate process
function ZNF_proc(x pitch _2Q)
local(ac bd x_LP)
instance(LP BP HP F F1 F2 a b c d)(
// once tabulated --> 3 adds 3 muls
F = 2 * sin(min($pi/2,ZNF_coef*2^(pitch*(1/12)))); // todo: tabulate
F1 = F * ZNF_g(HP); // gain for BP integrator (its input is HP)
F2 = F * ZNF_g(BP); // gain for LP integrator (its input is BP)
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// Matrix init and iterations
// 2 + 3i adds 1 + 5i muls
a = 1;
b = F1;
c = 1 - F1 * (_2Q + F2);
d = F2;
loop(ZNF_iter,
ac = a + c;
bd = b * d;
a = a*a - bd;
b = b * ac;
c = c*c - bd;
d = d * ac;
);
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
// Filter calc
// 6 adds 5 muls
x_LP = x - LP;
LP = x - a * x_LP + d * BP;
HP = x - _2Q * BP - LP;
BP = b * x_LP + c * BP;
);
//____________________________________________________________________
function sat3(x)local(x2)(
x = min(215/104, max(-215/104,x));
x2 = x*x;
x * (1
+ x2 * ( ( -411008.0/1987675.0)
+ x2 * ( ( 11581599744.0/459401384375.0)
+ x2 *( -5061276073984.0/4247165798546875.0))));
);
// sw ]-1 1[ typical ]-0.5, 0.5[
// sw = 0 => symetrical
// sw < 0 => short, long
// sw > 0 => long, short
function LFO_swing(x sw)(
x = x + x - 1;
x < sw ? (x+1)/(1+sw) : (x-sw)/(1-sw);
);
// sym [-1 1]
// sym = 0 => symetrical
// sym < 0 => short, long
// sym > 0 => long, short
function LFO_triangle(x sym)local(y)(
y = x + x - 1;
sym == 1 ? x
: sym == -1 ? 1 - x
: y <= sym ? (y+1)/(1+sym) : (1-y)/(1-sym);
);
// e [-1 1]
function LFO_exp(x e)(
e >= 0 ? x ^ (4*e + 1)
: 1 - (1 - x) ^ (-4*e + 1);
);
// a [0 1]
function LFO_sin(x a) local(y)(
y = x + x - 1;
y = 0.25 * y * (3 - y * y) + 0.5;
x + a * (y - x);
);
function LFO_val(y)(
y = LFO_swing(y, sl_sw1);
y = LFO_swing(y, sl_sw2);
y = LFO_swing(y, sl_sw3);
y = LFO_triangle(y, sl_sym);
y = LFO_exp(y, sl_exp);
LFO_sin(y, sl_sin);
);
//____________________________________________________________________
@block
play_state === 1 || play_state === 5 ?
beatPos = beat_position;
dBeatPos = tempo/(60 * srate);
//____________________________________________________________________
@slider
dp = -sl_rate/srate; // positive rate for down ramp
gain = 2^(sl_drive/6); // input gain
_gain = 2^(sl_outGain/6); // output gain
//____________________________________________________________________
@sample
beatPos += dBeatPos;
sl_rateScale>=1 ? p = (beatPos/4 - floor(beatPos/4)) * sl_rateScale // sync to 4/4
: sl_rateScale>=0.5 ? p = (beatPos/8 - floor(beatPos/8)) * sl_rateScale*2 // sync to 2*4/4
: p = (beatPos/16 - floor(beatPos/16)) * sl_rateScale*4; // sync to 4*4/4
p -= floor(p);
env = LFO_val(p);
pitch = sl_pitch0 + env * (sl_pitch1 - sl_pitch0);
_2Q = sl_D0 + env * ( sl_D1 - sl_D0);
l.ZNF_proc(gain * spl0, pitch,_2Q);
r.ZNF_proc(gain * spl1, pitch,_2Q);
sl_mode === 0 ? (spl0 = l.LP; spl1 = r.LP; )
: sl_mode === 1 ? (spl0 = l.BP; spl1 = r.BP; )
: sl_mode === 2 ? (spl0 = l.HP; spl1 = r.HP; );
spl0 *= _gain;
spl1 *= _gain;
spl0 = sat3(spl0);
spl1 = sat3(spl1);
@gfx
gfx.p = 0;
gfx.dp = 1/gfx_w;
gfx_r = 1; gfx_g = 1; gfx_b = 1;
gfx_x = -1;
loop(gfx_w,
gfx.y = LFO_val(gfx.p);
gfx_lineto(gfx_x + 1, gfx_h - gfx_h * gfx.y,1);
gfx.p += gfx.dp;
);
gfx_line(gfx_w*p, 0, gfx_w * p , gfx_h);
EDIT 02/10/15
I've found an analytic solution for the coefficients.. alleluia and thanks to Maxima. Hence, i think that the numerical iterative method is still the best way to calculate them, the analytic solution allows to check how accurate the numerical method is.
EDIT 02/10/15
Frequency response for the 3 outputs (2 iterations <-> x4 "virtual" oversampling)
Last edited by Smashed Transistors; 10-05-2015 at 12:47 PM.
|
|
|
Thread Tools |
|
Display Modes |
Linear Mode
|
Posting Rules
|
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts
HTML code is Off
|
|
|
All times are GMT -7. The time now is 05:06 AM.
|