Implement log2 function

This patch implements log2, which is necessary for amplitude modulation
(including amplitude LFO).

Also fixes up command line main.cc to track Dx7 init changes. (The
accuracy test is in that file).

cpp/src/core.gyp

@@ -11,6 +11,7 @@
         'fm_op_kernel.cc',
         'freqlut.cc',
         'lfo.cc',
+        'log2.cc',
         'patch.cc',
         'pitchenv.cc',
         'resofilter.cc',

cpp/src/log2.cc

@@ -0,0 +1,35 @@
+/*
+ * Copyright 2013 Google Inc.
+ * 
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * 
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ * 
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include <math.h>
+
+#include "synth.h"
+#include "log2.h"
+
+int32_t log2tab[LOG2_N_SAMPLES << 1];
+
+void Log2::init() {
+  const double mul = 1 / log(2);
+  for (int i = 0; i < LOG2_N_SAMPLES; i++) {
+    double y = mul * log(i + (LOG2_N_SAMPLES)) + (7 - LOG2_LG_N_SAMPLES);
+    log2tab[(i << 1) + 1] = (int32_t)floor(y * (1 << 24) + 0.5);
+  }
+  for (int i = 0; i < LOG2_N_SAMPLES - 1; i++) {
+    log2tab[i << 1] = log2tab[(i << 1) + 3] - log2tab[(i << 1) + 1];
+  }
+  log2tab[(LOG2_N_SAMPLES << 1) - 2] = (8 << 24) - log2tab[(LOG2_N_SAMPLES << 1) - 1];
+}
+

cpp/src/log2.h

@@ -0,0 +1,56 @@
+/*
+ * Copyright 2013 Google Inc.
+ * 
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * 
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ * 
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+class Log2 {
+ public:
+  Log2();
+
+  static void init();
+
+  // Q24 in, Q24 out
+  static int32_t lookup(uint32_t x);
+};
+
+#define LOG2_LG_N_SAMPLES 9
+#define LOG2_N_SAMPLES (1 << LOG2_LG_N_SAMPLES)
+
+#define LOG2_INLINE
+
+extern int32_t log2tab[LOG2_N_SAMPLES << 1];
+
+#ifdef __GNUC__
+static int clz(unsigned int x) {
+  return __builtin_clz(x);
+}
+// TODO: other implementations (including ANSI C)
+#endif
+
+#ifdef LOG2_INLINE
+inline
+int32_t Log2::lookup(uint32_t x) {
+  int exp = clz(x | 1);
+  unsigned int y = x << exp;
+
+  const int SHIFT = 31 - LOG2_LG_N_SAMPLES;
+  int lowbits = y & ((1 << SHIFT) - 1);
+  int y_int = (y >> (SHIFT - 1)) & ((LOG2_N_SAMPLES - 1) << 1);
+  int dz = log2tab[y_int];
+  int z0 = log2tab[y_int + 1];
+
+  int z = z0 + (((int64_t)dz * (int64_t)lowbits) >> SHIFT);
+  return z - (exp << 24);
+}
+#endif

cpp/src/main.cc

@@ -25,10 +25,13 @@
 #include "sawtooth.h"
 #include "sin.h"
 #include "exp2.h"
+#include "log2.h"
 #include "resofilter.h"
 #include "fm_core.h"
 #include "fm_op_kernel.h"
 #include "env.h"
+#include "patch.h"
+#include "controllers.h"
 #include "dx7note.h"
 
 using namespace std;
@@ -60,7 +63,22 @@ void test_sin_accuracy() {
     double err = fabs(y - yd);
     if (err > maxerr) maxerr = err;
   }
-  cout << "Max error: " << maxerr;
+  cout << "Max error: " << maxerr << endl;
+}
+
+void test_log_accuracy() {
+  double maxerr = 0;
+  for (int i = 0; i < 1000000; i++) {
+    uint32_t x = rand();
+    int32_t y = Log2::lookup(x);
+    double yd = (1 << 24) * log(x * (1.0 / (1 << 24))) / log(2);
+    double err = fabs(y - yd);
+    if (err > maxerr) {
+      maxerr = err;
+      cout << "x = " << x << ", y = " << y << ", yd = " << (int)yd << endl;
+    }
+  }
+  cout << "Max error: " << maxerr << endl;
 }
 
 void test_pure_accuracy() {
@@ -197,7 +215,12 @@ void mkdx7note(double sample_rate) {
   const int n_samples = 400 * 1024;
   WavOut w("/tmp/foo.wav", sample_rate, n_samples);
 
-  Dx7Note note(epiano, 57, 64);
+  Dx7Note note;
+  char unpacked_patch[156];
+  UnpackPatch(epiano, unpacked_patch);
+  note.init(unpacked_patch, 57, 64);
+  Controllers controllers;
+  controllers.values_[kControllerPitch] = 0x2000;
   int32_t buf[N];
 
   for (int i = 0; i < n_samples; i += N) {
@@ -207,7 +230,7 @@ void mkdx7note(double sample_rate) {
     if (i == n_samples / 2) {
       note.keyup();
     }
-     note.compute(buf);
+     note.compute(buf, 0, 0, &controllers);
     for (int j = 0; j < N; j++) {
       buf[j] >>= 2;
     }
@@ -236,17 +259,19 @@ int main(int argc, char **argv) {
   Sawtooth::init(sample_rate);
   Sin::init();
   Exp2::init();
+  Log2::init();
 
   //FmCore::dump();
   //test_sin_accuracy();
+  test_log_accuracy();
   //benchmark_fm_op();
   //test_pure_accuracy();
   //benchmark_sin();
   //int32_t freq = atoi(argv[1]);
   //cout << "Logfreq(" << freq << ") = " << Freqlut::lookup(freq) << endl;
 
-  //mkdx7note(sample_rate);
+  mkdx7note(sample_rate);
-  mksaw(sample_rate);
+  //mksaw(sample_rate);
   //test_ringbuffer();
   test_exp2();
   return 0;