iir: Start of testing code

Author: jonnor

benchmarks/iir/SimpleFilters.ipynb

@@ -11,12 +11,12 @@
     "%autoreload 2\n",
     "\n",
     "import sys\n",
-    "sys.path.insert(0, 'benchmarks/iir/')"
+    "#sys.path.insert(0, '.benchmarks/iir/')"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 42,
+   "execution_count": 2,
    "id": "10e5b8fc-6972-4502-a397-300a800df9ee",
    "metadata": {},
    "outputs": [],
@@ -27,12 +27,14 @@
     "from matplotlib import pyplot as plt\n",
     "\n",
     "import numpy\n",
-    "import pandas"
+    "import pandas\n",
+    "import librosa\n",
+    "import seaborn"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 43,
+   "execution_count": 30,
    "id": "fe122a40-66ed-4add-b59a-3254a35ece9a",
    "metadata": {},
    "outputs": [
@@ -54,7 +56,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/tmp/ipykernel_12537/4127207659.py:11: RuntimeWarning: divide by zero encountered in log10\n",
+      "/tmp/ipykernel_27591/4127207659.py:11: RuntimeWarning: divide by zero encountered in log10\n",
       "  db = 20 * numpy.log10(abs(h))\n"
      ]
     },
@@ -118,17 +120,17 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 53,
+   "execution_count": 5,
    "id": "0dff57e5-69cf-4ad7-80f0-d4b721a77910",
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "array('i', [29321, 0, -58643, 29321, -58280, 26239, 2210, 0, 4420, 2210, -37452, 13526])"
+       "array('i', [26233, 0, -52466, 26233, -51149, 21015])"
       ]
      },
-     "execution_count": 53,
+     "execution_count": 5,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -166,7 +168,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 48,
+   "execution_count": 6,
    "id": "7857c33a-6f7f-405b-9d01-16e9d86043f3",
    "metadata": {},
    "outputs": [
@@ -184,7 +186,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/tmp/ipykernel_12537/4127207659.py:11: RuntimeWarning: divide by zero encountered in log10\n",
+      "/tmp/ipykernel_27591/4127207659.py:11: RuntimeWarning: divide by zero encountered in log10\n",
       "  db = 20 * numpy.log10(abs(h))\n"
      ]
     },
@@ -222,13 +224,98 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 13,
    "id": "4cb895db-d746-4af0-8e12-e7805865586d",
    "metadata": {},
    "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 43,
+   "id": "51a151c4-cf39-430d-bd49-c0f8fa0baab3",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "run:  micropython /home/jon/projects/emlearn-micropython/benchmarks/iir/iir_run.py /tmp/tmp3kr9y4bm/input.npy /tmp/tmp3kr9y4bm/output.npy /tmp/tmp3kr9y4bm/filter.npy\n",
+      "b''\n"
+     ]
+    },
+    {
+     "ename": "AssertionError",
+     "evalue": "/tmp/tmp3kr9y4bm/output.npy",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
+      "Cell \u001b[0;32mIn[43], line 32\u001b[0m\n\u001b[1;32m     30\u001b[0m noise \u001b[38;5;241m=\u001b[39m numpy\u001b[38;5;241m.\u001b[39mrandom\u001b[38;5;241m.\u001b[39mrandom(\u001b[38;5;241m1000\u001b[39m)\n\u001b[1;32m     31\u001b[0m out_scipy \u001b[38;5;241m=\u001b[39m scipy\u001b[38;5;241m.\u001b[39msignal\u001b[38;5;241m.\u001b[39msosfilt(coeff, noise)\n\u001b[0;32m---> 32\u001b[0m out_emlearn \u001b[38;5;241m=\u001b[39m \u001b[43msosfilt_emlearn\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnumpy\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43marray\u001b[49m\u001b[43m(\u001b[49m\u001b[43mcoeff\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnoise\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     34\u001b[0m fig, ax \u001b[38;5;241m=\u001b[39m plt\u001b[38;5;241m.\u001b[39msubplots(\u001b[38;5;241m1\u001b[39m, figsize\u001b[38;5;241m=\u001b[39m(\u001b[38;5;241m10\u001b[39m, \u001b[38;5;241m5\u001b[39m))\n\u001b[1;32m     35\u001b[0m plot_spectrum_difference(ax, noise, out_scipy, fs\u001b[38;5;241m=\u001b[39mfs, label\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mscipy\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
+      "Cell \u001b[0;32mIn[43], line 24\u001b[0m, in \u001b[0;36msosfilt_emlearn\u001b[0;34m(coeff, input)\u001b[0m\n\u001b[1;32m     19\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m     20\u001b[0m \u001b[38;5;124;03mCalls MicroPython in a subprocess, using eml_iir modules from emlearn-micropython\u001b[39;00m\n\u001b[1;32m     21\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m     22\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01miir_run_subprocess\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m run_iir\n\u001b[0;32m---> 24\u001b[0m out \u001b[38;5;241m=\u001b[39m \u001b[43mrun_iir\u001b[49m\u001b[43m(\u001b[49m\u001b[43marray\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43marray\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mf\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43minput\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43marray\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43marray\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mf\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcoeff\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mflatten\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     26\u001b[0m out \u001b[38;5;241m=\u001b[39m numpy\u001b[38;5;241m.\u001b[39marray(out)\n\u001b[1;32m     28\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m out\n",
+      "File \u001b[0;32m~/projects/emlearn-micropython/benchmarks/iir/iir_run_subprocess.py:48\u001b[0m, in \u001b[0;36mrun_iir\u001b[0;34m(data, coefficients, micropython_bin)\u001b[0m\n\u001b[1;32m     45\u001b[0m \u001b[38;5;28mprint\u001b[39m(out)\n\u001b[1;32m     47\u001b[0m \u001b[38;5;66;03m# load the output\u001b[39;00m\n\u001b[0;32m---> 48\u001b[0m \u001b[38;5;28;01massert\u001b[39;00m os\u001b[38;5;241m.\u001b[39mpath\u001b[38;5;241m.\u001b[39mexists(output_path), output_path\n\u001b[1;32m     49\u001b[0m output_shape, output \u001b[38;5;241m=\u001b[39m npyfile\u001b[38;5;241m.\u001b[39mload(output_path)\n\u001b[1;32m     51\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m output\n",
+      "\u001b[0;31mAssertionError\u001b[0m: /tmp/tmp3kr9y4bm/output.npy"
+     ]
+    }
+   ],
    "source": [
-    "# TODO: run tests with eml_iir and eml_iir_q15, using scipy sosfilt as reference"
+    "\n",
+    "def spectrum_welch(signal, sr, n_fft, window='hann'):\n",
+    "    freqs, power = scipy.signal.welch(signal, fs=sr, nperseg=n_fft, window=window, scaling=\"spectrum\", average='mean')\n",
+    "    db = librosa.power_to_db(power, ref=0.0, top_db=240)\n",
+    "    return pandas.Series(db, index=freqs)\n",
+    "\n",
+    "\n",
+    "def plot_spectrum_difference(ax, input, output, fs, label=None, alpha=0.5, n_fft = 512):\n",
+    "    \n",
+    "    S_in = spectrum_welch(input, sr=fs, n_fft=n_fft)\n",
+    "    S_out = spectrum_welch(output, sr=fs, n_fft=n_fft)\n",
+    "    \n",
+    "    FR = S_out - S_in\n",
+    "    FR += numpy.random.random(1)\n",
+    "    ax.plot(FR.index, FR.values, label=label, alpha=alpha)\n",
+    "    ax.set_xscale('log')\n",
+    "    ax.set_xlim(1, fs)\n",
+    "\n",
+    "def sosfilt_emlearn(coeff : numpy.array, input : numpy.array):\n",
+    "    \"\"\"\n",
+    "    Calls MicroPython in a subprocess, using eml_iir modules from emlearn-micropython\n",
+    "    \"\"\"\n",
+    "    from iir_run_subprocess import run_iir\n",
+    "    \n",
+    "    out = run_iir(array.array('f', input), array.array('f', coeff.flatten()))\n",
+    "\n",
+    "    out = numpy.array(out)\n",
+    "    \n",
+    "    return out\n",
+    "\n",
+    "noise = numpy.random.random(1000)\n",
+    "out_scipy = scipy.signal.sosfilt(coeff, noise)\n",
+    "out_emlearn = sosfilt_emlearn(numpy.array(coeff), noise)\n",
+    "\n",
+    "fig, ax = plt.subplots(1, figsize=(10, 5))\n",
+    "plot_spectrum_difference(ax, noise, out_scipy, fs=fs, label='scipy')\n",
+    "plot_spectrum_difference(ax, noise, out_emlearn, fs=fs, label='emlearn')\n",
+    "ax.legend()\n",
+    "\n",
+    "#numpy.save('noise.npy', noise)\n"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e1597230-78a9-45e2-a1af-e8bdc0e5388e",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7546d0f2-29fa-420c-8b0a-13faf8919fc7",
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {

benchmarks/iir/iir_run.py

@@ -0,0 +1,55 @@
+
+import sys
+import npyfile
+
+def iir_process_file(inp, out, filters, chunksize):
+
+    coefficients_shape, coefficients = npyfile.load(filters)
+    print('c', coefficients_shape, coefficients)
+
+    # NOTE: filter shape and type checked by the C modules
+
+    # Processing is done by streaming, to keep RAM usage low
+    with npyfile.Reader(inp) as reader:
+
+        # Check inputs and setup filters
+        print('r', reader.shape, reader.typecode, reader.itemsize)
+
+        if len(reader.shape) != 1:
+            raise ValueError("Input must be 1d")
+        if reader.typecode == 'f':
+            import emliir
+            filter = emliir.new(coefficients)
+        elif reader.typecode == 'h':
+            import eml_iir_q15
+            filter = eml_iir_q15.new(coefficients)
+        else:
+            raise ValueError("Input must either be float32/f or int16/h")
+
+
+        with npyfile.Writer(inp, reader.shape, reader.typecode) as writer:
+
+            # Apply filters
+            for chunk in reader.read_data_chunks(chunksize):
+                print(len(chunk))
+
+                filter.run(chunk) # operates in-place
+
+                writer.write_values(chunk)
+
+
+def main():
+
+    args = sys.argv
+    if len(args) < 3:
+        print('Usage: micropython iir_run.py INPUT.npy OUTPUT.npy FILTERS.npy')
+
+    _, inp, out, filters = args
+
+    chunksize = 100
+
+    iir_process_file(inp, out, filters, chunksize=chunksize)
+
+if __name__ == '__main__':
+    main()
+    

benchmarks/iir/iir_run_subprocess.py

@@ -0,0 +1,51 @@
+
+import os
+import array
+import tempfile
+import subprocess
+
+import npyfile
+
+
+here = os.path.dirname(__file__)
+
+
+def run_iir(data : array.array,
+        coefficients : array.array,
+        micropython_bin='micropython',
+    ):
+
+    script_path = os.path.join(here, 'iir_run.py')
+
+    with tempfile.TemporaryDirectory() as temp_dir:
+
+        filter_path = os.path.join(temp_dir, 'filter.npy')
+        input_path = os.path.join(temp_dir, 'input.npy')
+        output_path = os.path.join(temp_dir, 'output.npy')
+
+        # write input data to files
+        npyfile.save(filter_path, coefficients)
+        npyfile.save(input_path, data)
+
+        assert not os.path.exists(output_path), output_path
+
+        # run the processing function
+        args = [
+            micropython_bin,
+            script_path,
+            input_path,
+            output_path,
+            filter_path,            
+        ]
+        cmd = ' '.join(args)
+
+        print('run: ', cmd)
+
+        out = subprocess.check_output(args)
+        print(out)
+
+        # load the output
+        assert os.path.exists(output_path), output_path
+        output_shape, output = npyfile.load(output_path)
+
+        return output