// Music 220a - Final Project - Romain Michon
// 12/11/2012
// Plate physical model based on a nonlinear feedback delay network. This model was primilarly implemented in effect.lib (Faust distribution) by Julius Smith and Romain Michon. This version is slightly different than the one available in Faust. Plotting the block diagram (faust2svg) will probably help you understand how it works. Sorry the comments in french :D... (I wrote this file for master thesis)

declare name "nlfFdnRev";

import("filter.lib");
import("music.lib");

  N = 4; // ordre du reseau de ligne de retard recursives (valeur multiple de 2 inférieure à 16)
 NB =  3; // nombre de duree de resonance controllable de bande de frequence (trois ou plus)
BSO =  1; // ordre des filtres passes-tout et passes-bas pour le partage de bandes (entier positif impaire)

nonl = 0.1; // taux de non-linearites par deffaut

//utilisation de la fonction fdnrev1_demo (reimplementation de fdnrev0_demo de effect.lib)
process = fdnrev1_demo(N,NB,BSO);

//re-declaration de fdnrev0 de effect.lib utilisant le filtre passe-tout non-lineaire allpassnn (filter.lib) au niveau du signal recursif des chaines de retards
fdnrev1(MAXDELAY, delays, BBSO, freqs, durs, theta, loopgainmax, nonl) 
  = (bus(2*N) :> bus(N) : delaylines(N)) ~ 
    (delayfilters(N,freqs,durs) : feedbackmatrix(N)) 
with {
  N = count(delays);
  NB = count(durs);
  delayval(i) = take(i+1,delays);
  dlmax(i) = MAXDELAY;
  delaylines(N) = par(i,N,(delay(dlmax(i),(delayval(i)-1))));
  delayfilters(N,freqs,durs) = par(i,N,filter(i,freqs,durs));
  feedbackmatrix(N) = hadamard(N);
  vbutterfly(n) = bus(n) <: (bus(n):>bus(n/2)) , ((bus(n/2),(bus(n/2):par(i,n/2,*(-1)))) :> bus(n/2));

  filter(i,freqs,durs) = filterbank(BBSO,freqs) : par(j,NB,*(g(j,i)))
                         :> *(loopgainmax) / sqrt(N) : nonlinallpass
  with {
    dur(j) = take(j+1,durs);
    n60(j) = dur(j)*SR; // decay time in samples
    g(j,i) = exp(-3.0*log(10.0)*delayval(i)/n60(j));
    //filtre passe-tout non-lineaire, nonl est le taux de non-linearites (valeur entre -1 et 1)
    s = nonl*PI;
    nonlinallpass(x) = x:allpassnn(3,(s*x,s*x*x,s*x*x*x)); // filter.lib
  };
};

//*************************************************************************************************
// La suite du code n'est pas commentee (re-declaration a l'identique des fonctions de effect.lib)
//*************************************************************************************************

// ---------- prime_power_delays -----

prime_power_delays(N,pathmin,pathmax) = par(i,N,delayvals(i)) with {
  Np = 16;
  primes = 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53;
  prime(n) = primes : selector(n,Np); // math.lib

  // Prime Power Bounds [matlab: floor(log(maxdel)./log(primes(53)))]
  maxdel=8192; // more than 63 meters at 44100 samples/sec & 343 m/s
  ppbs = 13,8,5,4, 3,3,3,3, 2,2,2,2, 2,2,2,2; // 8192 is enough for all
  ppb(i) = take(i+1,ppbs);

  // Approximate desired delay-line lengths using powers of distinct primes:
  c = 343; // soundspeed in m/s at 20 degrees C for dry air
  dmin = SR*pathmin/c;
  dmax = SR*pathmax/c;
  dl(i) = dmin * (dmax/dmin)^(i/float(N-1)); // desired delay in samples
  ppwr(i) = floor(0.5+log(dl(i))/log(prime(i))); // best prime power
  delayvals(i) = prime(i)^ppwr(i); // each delay a power of a distinct prime
};

//--------------------- stereo_reverb_tester --------------------
// Handy test inputs for reverberator demos below.

stereo_reverb_tester(revin_group,x,y) = inx,iny with {
  ck_group(x) = revin_group(vgroup("[1] Input Config",x));
  mutegain = 1 - ck_group(checkbox("[1] Mute Ext Inputs
         [tooltip: When this is checked, the stereo external audio inputs are disabled (good for hearing the impulse response or pink-noise response alone)]"));
  pinkin = ck_group(checkbox("[2] Pink Noise 
         [tooltip: Pink Noise (or 1/f noise) is Constant-Q Noise (useful for adjusting the EQ sections)]"));

  impulsify = _ <: _,mem : - : >(0);
  imp_group(x) = revin_group(hgroup("[2] Impulse Selection",x));
  pulseL =  imp_group(button("[1] Left
         [tooltip: Send impulse into LEFT channel]")) : impulsify;
  pulseC =  imp_group(button("[2] Center
         [tooltip: Send impulse into LEFT and RIGHT channels]")) : impulsify;
  pulseR =  imp_group(button("[3] Right
         [tooltip: Send impulse into RIGHT channel]")) : impulsify;
  
  inx = x*mutegain + (pulseL+pulseC) + pn;
  iny = y*mutegain + (pulseR+pulseC) + pn;
  pn = 0.1*pinkin*component("oscillator.lib").pink_noise;
};

fdnrev1_demo(N,NB,BBSO,x,y) = stereo_reverb_tester(revin_group,(x*inputGain*strPos),(y*inputGain*(1-strPos)))
	  <: fdnrev1(MAXDELAY,delays,BBSO,freqs,durs,theta,
	             loopgainmax,nonl) 
          :> *(gain) : echo
with {
  echo = +~(fdelay(65536, (int(hslider("echoDur", 0, 0, 1000, 0.10)*millisec)-1 : smooth(0.999))) * (hslider("echoFeedback", 0, 0, 100, 0.1)/100.0 : smooth(0.999)));
  MAXDELAY = 8192; // sync w delays and prime_power_delays above
  defdurs = (2.6,2.6,2.6); // NB default durations (sec)
  deffreqs = (500,1000);   // NB-1 default crossover frequencies (Hz)
  deflens = (56.3,63.0); // 2 default min and max path lengths

  fdn_group(x)  = vgroup("FEEDBACK DELAY NETWORK (FDN) REVERBERATOR, ORDER 16
  [tooltip: See Faust's effect.lib for documentation and references]", x);

  freq_group(x)  = fdn_group(hgroup("[1] Band Crossover Frequencies", x));
  t60_fbm_group(x)  = fdn_group(hgroup("[2] Feedback Control", x));
  t60_group(x)  = t60_fbm_group(hgroup("[1] Band Decay Times (T60)", x));
  fbm_group(x)  = t60_fbm_group(hgroup("[2] FBM", x));
  path_group(x)  = fdn_group(hgroup("[3] Room Dimensions", x));
  revin_group(x)  = fdn_group(hgroup("[4] Input Controls", x));
  nonl_group(x) = revin_group(hgroup("[4] Nonlinearity",x));
  quench_group(x) = revin_group(hgroup("[3] Reverb State",x));

  nonl = hslider("nonlinearity",0.999, -0.999, 0.999, 0.001);
  strPos = hslider("strikePosition",0.5,0,1,0.01) : smooth(0.999);
  loopgainmax = 0.99999; 
  pathmin = hslider("minAcRayL",0.6, 0.1, 63, 0.1);
  pathmax = hslider("maxAcRayL",9.4, 0.1, 63, 0.1);
  durvals(i) = hslider("bandDecayT%i",take(i+1,defdurs), 0.1, 10, 0.1);
  durs = par(i,NB,durvals(NB-1-i));
  theta = 0.25*PI*fbm_group(hslider("FBM", 1, 0, 1, 0.01));

  freqvals(i) = hslider("bandCross%i",take(i+1,deffreqs), 100, 10000, 1);
  freqs = par(i,NB-1,freqvals(i));

  delays = prime_power_delays(N,pathmin,pathmax);
inputGain = hslider("inputGain", 0.9, 0, 1, 0.01);
  gain = hslider("outputLevel", -40, -70, 40, 0.1) : db2linear;
};