Help with Runge-Kutta Filter

[Ichad.c]

Can't seem to get my head around this filter on musicdsp.org. The CSound notation is confusing me - my head feels like a cat in tumble-dryer!

Here's the code.

Code: Select all

;---------------------------------------------------------
; Runge-Kutta Freaky Filter
;---------------------------------------------------------
       instr  7

idur   =      p3            ; Duration
iamp   =      p4            ; Amplitude
kfco2  zkr    p5            ; Filter cutoff
kq1    zkr    p6            ; Q
ih     =      .001          ; Diff eq step size
ipanl  =      sqrt(p8)      ; Pan left
ipanr  =      sqrt(1-p8)    ; Pan right
ifqc   =      cpspch(p9)    ; Pitch to frequency
kpb    zkr    p10           ; Pentic bounce frequency
kpa    zkr    p11           ; Pentic bounce amount
kasym  zkr    p12           ; Q assymmetry amount
kasep  zkr    p13           ; Q asymmetry separation


kdclck linseg 0, .02, 1, idur-.04, 1, .02, 0 ; Declick envelope
kfco1  expseg 1, idur, .1
kfco   =      kfco1*kfco2
kq     =      kq1*kfco^1.2*.1
kfc    =      kfco/8/sr*44100

ay     init   0
ay1    init   0
ay2    init   0
ay3    init   0
axs    init   0
avxs   init   0
ax     vco   1, ifqc, 2, 1, 1, 1        ; Square wave

  ; R-K Section 1
  afdbk =      kq*ay/(1+exp(-ay*3*kasep)*kasym)   ; Only oscillate in one
direction
  ak11  =      ih*((ax-ay1)*kfc-afdbk)
  ak21  =      ih*((ax-(ay1+.5*ak11))*kfc-afdbk)
  ak31  =      ih*((ax-(ay1+.5*ak21))*kfc-afdbk)
  ak41  =      ih*((ax-(ay1+ak31))*kfc-afdbk)
  ay1   =      ay1+(ak11+2*ak21+2*ak31+ak41)/6

  ; R-K Section 2
  ak12  =      ih*((ay1-ay2)*kfc)
  ak22  =      ih*((ay1-(ay2+.5*ak12))*kfc)
  ak32  =      ih*((ay1-(ay2+.5*ak22))*kfc)
  ak42  =      ih*((ay1-(ay2+ak32))*kfc)
  ay2   =      ay2+(ak12+2*ak22+2*ak32+ak42)/6

  ; Pentic bounce equation
  ax3    =      -.1*ay*kpb
  aaxs   =      (ax3*ax3*ax3*ax3*ax3+ay2)*1000*kpa     ; Update acceleration

  ; R-K Section 3
  ak13  =      ih*((ay2-ay3)*kfc+aaxs)
  ak23  =      ih*((ay2-(ay3+.5*ak13))*kfc+aaxs)
  ak33  =      ih*((ay2-(ay3+.5*ak23))*kfc+aaxs)
  ak43  =      ih*((ay2-(ay3+ak33))*kfc+aaxs)
  ay3   =      ay3+(ak13+2*ak23+2*ak33+ak43)/6

  ; R-K Section 4
  ak14  =      ih*((ay3-ay)*kfc)
  ak24  =      ih*((ay3-(ay+.5*ak14))*kfc)
  ak34  =      ih*((ay3-(ay+.5*ak24))*kfc)
  ak44  =      ih*((ay3-(ay+ak34))*kfc)
  ay    =      ay+(ak14+2*ak24+2*ak34+ak44)/6

aout   =      ay*iamp*kdclck*.07 ; Apply amp envelope and declick

       outs   aout*ipanl, aout*ipanr  ; Output the sound

       endin

These are the parts I don't get:

Code: Select all

kfco1  expseg 1, idur, .1
kfco   =      kfco1*kfco2

What is kfco1?

Code: Select all

afdbk =      kq*ay/(1+exp(-ay*3*kasep)*kasym)   ; Only oscillate in one
direction

WTF?

Code: Select all

; Pentic bounce equation
  ax3    =      -.1*ay*kpb
  aaxs   =      (ax3*ax3*ax3*ax3*ax3+ay2)*1000*kpa     ; Update acceleration

Another WTF?

Can anybody explain this in plain c++?

Regards
Andrew

[Ichad.c]

Or any other filter example that uses the RK algo?

[Richard_Synapse]

Looks like a diode ladder filter, weird indeed. Might be easiest to dig up the paper this was based on (if that exists).

Richard

[Ichad.c]

Looks like a diode ladder filter, weird indeed. Might be easiest to dig up the paper this was based on (if that exists).

Richard

Yip, I also think it's a ladder. As far as I can figure - Runge-Kutta(4) is the equation solver. Also seems that it upsamples to solve it's equation(?).
And is seems like unbuffered stages(?) like a 4th-order filter - which is uncommon in itself. Usually 2nd order stages are cascade or 1st orders with cascade + feedbck. My observations might be totally whack though.

Andrew

[mystran]

Diode ladder has 4 stages like a transistor ladder, but unbuffered because of the way the transistors are connected. On top of that, the transistor pairs are still non-linear, so the stages are connected unbuffered through non-linear stages, which makes it a pain to discretize properly.

That's probably why it's using RK4. I experimented with similar stuff for my "Dolphin bassline" but ended up using a 2nd order method with just enough oversampling that it sorta-kinda-works since the "upsampling" that RK4 does doesn't really help with aliasing, so on top of that you need to oversample anyway, and it gets a bit costly (if you're doing full non-linearities).

[Ichad.c]

Yeah, it seems that RK4 is sluggish, but I still don't get the CSound stuff, will still like to implement it - even if it is for educational purposes.

Thanks for the insight though.

Regards
Andrew