Prototype of nonlinear ladder filter

This patch contains a test implementation of a nonlinear ladder filter,
closely following Antti Huovilainen's 2004 paper. It's hidden behind an
ifdef, though, as it's not ready for prime time.

cpp/src/resofilter.cc

@@ -14,8 +14,16 @@
  * limitations under the License.
  */
 
+// Resonant filter implementation. This closely follows "Non-Linear
+// Digital Implementation of the Moog Ladder Filter" by Antti
+// Huovilainen, 2004.
+
+// The full implementation requires both a tuning table and 2x
+// oversampling, neither of which are present yet, but we'll get there. 
+
 #include "synth.h"
 #include "freqlut.h"
+#include "exp2.h"
 #include "resofilter.h"
 
 double this_sample_rate;
@@ -27,6 +35,9 @@ void ResoFilter::init(double sample_rate) {
 ResoFilter::ResoFilter() {
   for (int i = 0; i < 4; i++) {
     x[i] = 0;
+#if defined(NONLINEARITY)
+    w[i] = 0;
+#endif
   }
 }
 
@@ -51,23 +62,52 @@ void ResoFilter::process(const int32_t **inbufs, const int32_t *control_in,
   }
   const int32_t *ibuf = inbufs[0];
   int32_t *obuf = outbufs[0];
-  int x0 = x[0];
-  int x1 = x[1];
-  int x2 = x[2];
-  int x3 = x[3];
+  int32_t x0 = x[0];
+  int32_t x1 = x[1];
+  int32_t x2 = x[2];
+  int32_t x3 = x[3];
+#if defined(NONLINEARITY)
+  int32_t w0 = w[0];
+  int32_t w1 = w[1];
+  int32_t w2 = w[2];
+  int32_t w3 = w[3];
+  int32_t yy0 = yy;
+#endif;
   for (int i = 0; i < n; i++) {
     alpha += delta_alpha;
     k += delta_k;
     int32_t signal = ibuf[i];
+#if defined(NONLINEARITY) 
+    int32_t fb = ((int64_t)k * (int64_t)(x3 + yy0)) >> 25;
+    yy0 = x3;
+    int32_t rx = signal - fb;
+    int32_t trx = Tanh::lookup(rx);
+    x0 = x0 + ((((int64_t)(trx - w0) * (int64_t)alpha)) >> 24);
+    w0 = Tanh::lookup(x0);
+    x1 = x1 + ((((int64_t)(w0 - w1) * (int64_t)alpha)) >> 24);
+    w1 = Tanh::lookup(x1);
+    x2 = x2 + ((((int64_t)(w1 - w2) * (int64_t)alpha)) >> 24);
+    w2 = Tanh::lookup(x2);
+    x3 = x3 + ((((int64_t)(w2 - w3) * (int64_t)alpha)) >> 24);
+    w3 = Tanh::lookup(x3);
+#else
     int32_t fb = ((int64_t)k * (int64_t)x3) >> 24;
     x0 = x0 + ((((int64_t)(signal - fb - x0) * (int64_t)alpha)) >> 24);
     x1 = x1 + ((((int64_t)(x0 - x1) * (int64_t)alpha)) >> 24);
     x2 = x2 + ((((int64_t)(x1 - x2) * (int64_t)alpha)) >> 24);
     x3 = x3 + ((((int64_t)(x2 - x3) * (int64_t)alpha)) >> 24);
+#endif
     obuf[i] = x3;
   }
   x[0] = x0;
   x[1] = x1;
   x[2] = x2;
   x[3] = x3;
+#if defined(NONLINEARITY)
+  w[0] = w0;
+  w[1] = w1;
+  w[2] = w2;
+  w[3] = w3;
+  yy = yy0;
+#endif
 }

cpp/src/resofilter.h

@@ -29,6 +29,10 @@ class ResoFilter : Module {
 			   const int32_t *control_last, int32_t **outbufs);
  private:
   int32_t x[4];
+#if defined(NONLINEARITY)
+  int32_t w[4];
+  int32_t yy;
+#endif
 };
 
 #endif  // SYNTH_RESOFILTER_H_