From 0d217e93264077a05b388b77cf23c4f8f6a26a98 Mon Sep 17 00:00:00 2001
From: kris <kris.kennaway@gmail.com>
Date: Sat, 2 Jul 2022 20:53:02 +0100
Subject: [PATCH] Improved quality and speed

---
 encode_audio.py | 16 +++++------
 lookahead.pyx   | 74 +++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 82 insertions(+), 8 deletions(-)
 create mode 100644 lookahead.pyx

diff --git a/encode_audio.py b/encode_audio.py
index 3bb993a..0704594 100755
--- a/encode_audio.py
+++ b/encode_audio.py
@@ -133,14 +133,13 @@ def audio_bytestream(data: numpy.ndarray, step: int, lookahead_steps: int,
 
     clicks = 0
     min_lookahead_steps = lookahead_steps
-    while i < dlen // 20:
+    while i < dlen // 1:
         # XXX handle end of data cleanly
         if i >= next_tick:
             eta.print_status()
             next_tick = int(eta.i * dlen / 1000)
 
-        # XXX
-        if frame_offset >= 2045:  # XXX
+        if frame_offset >= 2043:  # XXX
             lookahead_steps = min_lookahead_steps + 120  # XXX parametrize
         else:
             lookahead_steps = min_lookahead_steps
@@ -150,14 +149,15 @@ def audio_bytestream(data: numpy.ndarray, step: int, lookahead_steps: int,
             opcodes.candidate_opcodes(
                 frame_horizon(frame_offset, lookahead_steps),
                 lookahead_steps, opcode if frame_offset == 2047 else None)
-        all_positions = lookahead.evolve(
-            sp, y1, y2, voltage1, voltage2, voltage1 * voltages)
+        opcode_idx = lookahead.evolve_return_best(
+            sp, y1, y2, voltage1, voltage2, voltage1 * voltages,
+            data[i:i+lookahead_steps])
 
         # Pick the opcode sequence that minimizes the total squared error
         # relative to the data waveform.
-        errors = total_error(
-            all_positions * sp.scale, data[i:i + lookahead_steps])
-        opcode_idx = numpy.argmin(errors).item()
+        # errors = total_error(
+        #     all_positions * sp.scale, data[i:i + lookahead_steps])
+        # opcode_idx = numpy.argmin(errors).item()
         # if frame_offset == 2046:
         #     print("XXX", lookahead_steps)
         #     print(opcode_idx)
diff --git a/lookahead.pyx b/lookahead.pyx
new file mode 100644
index 0000000..655efc3
--- /dev/null
+++ b/lookahead.pyx
@@ -0,0 +1,74 @@
+# cython: infer_types=True
+# cython: profile=False
+# cython: boundscheck=False
+# cython: wraparound=False
+
+
+cimport cython
+from libc.stdlib cimport malloc, free
+import numpy as np
+
+
+def evolve_return_best(object speaker, float position1, float position2, float voltage1, float voltage2, float[:, ::1] voltages, float[::1] data):
+    cdef double c1 = speaker.c1
+    cdef double c2 = speaker.c2
+    cdef double b1 = speaker.b1
+    cdef double b2 = speaker.b2
+    cdef double scale = speaker.scale
+
+    cdef int i, j
+    cdef double y, y1, y2
+    cdef float x1, x2
+
+    cdef int lowest_idx
+    cdef double lowest_err = 1e9
+    cdef double error
+    for i in range(voltages.shape[0]):
+        x1 = voltage1
+        x2 = voltage2
+        y1 = position1
+        y2 = position2
+        error = 0
+        for j in range(voltages.shape[1]):
+            y = c1 * y1 - c2 * y2 + b1 * x1 + b2 * x2
+            error += (y * scale - data[j]) ** 2
+            if error > lowest_err:
+                break
+            # output[i, j] = y
+            x2 = x1
+            x1 = voltages[i, j]  # XXX does this really always lag?
+            y2 = y1
+            y1 = y
+        if error < lowest_err:
+            lowest_err = error
+            lowest_idx = i
+
+    return lowest_idx
+
+
+def evolve(object speaker, float position1, float position2, float voltage1, float voltage2, float[:, ::1] voltages):
+    cdef double[:,::1] output = np.empty_like(voltages, dtype=np.float64)
+
+    cdef double c1 = speaker.c1
+    cdef double c2 = speaker.c2
+    cdef double b1 = speaker.b1
+    cdef double b2 = speaker.b2
+
+    cdef int i, j
+    cdef double y, y1, y2
+    cdef float x1, x2
+
+    for i in range(voltages.shape[0]):
+        x1 = voltage1
+        x2 = voltage2
+        y1 = position1
+        y2 = position2
+        for j in range(voltages.shape[1]):
+            y = c1 * y1 - c2 * y2 + b1 * x1 + b2 * x2
+            output[i, j] = y
+            x2 = x1
+            x1 = voltages[i, j]  # XXX does this really always lag?
+            y2 = y1
+            y1 = y
+
+    return output