Copy over LadderFilter implementation

.gitignore

@@ -6,3 +6,5 @@ docs/.jekyll-cache
 .cache
 
 .DS_Store
+
+__pycache__

CMakeLists.txt

@@ -45,6 +45,7 @@ target_include_directories(BlackBird
 target_sources(BlackBird
     PRIVATE
     source/dsp/DSPParameters.cpp
+    source/dsp/LadderFilter.cpp
     source/ui/EditorHeader.cpp
     source/ui/Knob.cpp
     source/ui/PluginEditor.cpp

scripts/aliasing_analysis/README.md

@@ -0,0 +1,21 @@
+# Exploring LadderFilter aliasing issue
+
+`tanh_impact.ipynb` demonstrates how `tanh(x)` function produces new
+frequency components and how the aliasing phenomenon looks like at various
+frequencies and magnitudes.
+
+## Configuring Python environment
+
+1. [Install miniconda](https://docs.anaconda.com/free/miniconda/miniconda-install/)
+
+2. Create `black-bird` conda environment:
+```
+conda env create -f environment.yml
+```
+
+3. Activate the environment:
+```
+conda activate black-bird
+```
+
+VSCode is the recommended editor/IDE.

scripts/aliasing_analysis/antiderivatives.ipynb

@@ -0,0 +1,95 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from plotly import graph_objs as go\n",
+    "\n",
+    "import tanh_antiderivative\n",
+    "import dsplib"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import importlib\n",
+    "importlib.reload(dsplib)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fs_hz = 48e3\n",
+    "signal_duration_sec = 0.2\n",
+    "smp_num = dsplib.calc_smp_num(signal_duration_sec, fs_hz)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "S = dsplib.generate_delayed_sin_matrix(smp_num=smp_num, tone_freq_n=9100/fs_hz, mag=0.5, history_smp_num=3, noise_level=1e-4)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ad_1_sig = tanh_antiderivative.order1(S)\n",
+    "ad_2_sig = tanh_antiderivative.order2(S)\n",
+    "ad_3_sig = tanh_antiderivative.order3(S)\n",
+    "\n",
+    "f, S_in = dsplib.calc_spectrum(S[0,:], fs=fs_hz)\n",
+    "_, S_tanh = dsplib.calc_spectrum(np.tanh(S[0,:]), fs=fs_hz)\n",
+    "_, S_ad_1 = dsplib.calc_spectrum(ad_1_sig, fs=fs_hz)\n",
+    "_, S_ad_2 = dsplib.calc_spectrum(ad_2_sig, fs=fs_hz)\n",
+    "_, S_ad_3 = dsplib.calc_spectrum(ad_3_sig, fs=fs_hz)\n",
+    "\n",
+    "\n",
+    "fig = go.Figure()\n",
+    "fig.add_trace(go.Scatter(x=f, y=S_in, name=\"input\"))\n",
+    "fig.add_trace(go.Scatter(x=f, y=S_tanh, line=dict(width=2, dash='dash'), name=\"tanh\"))\n",
+    "fig.add_trace(go.Scatter(x=f, y=S_ad_1, line=dict(width=2, dash='dot'), name=\"ad_1\"))\n",
+    "fig.add_trace(go.Scatter(x=f, y=S_ad_2, line=dict(width=1.5, dash='dot'), name=\"ad_2\"))\n",
+    "fig.add_trace(go.Scatter(x=f, y=S_ad_3, line=dict(width=1, dash='dot'), name=\"ad_3\"))\n",
+    "fig.update_layout(hovermode=\"x unified\")\n",
+    "fig.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "black-bird",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.8"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

scripts/aliasing_analysis/dsplib.py

@@ -0,0 +1,65 @@
+import numpy as np
+from scipy import fft
+
+
+def calc_smp_num(signal_duration_sec, fs_hz):
+    return int(signal_duration_sec*fs_hz)
+
+
+def generate_time_signal(signal_duration_sec, fs_hz, history_smp_num=0):
+    return np.arange(-history_smp_num/fs_hz, signal_duration_sec, 1/fs_hz, dtype=np.float64)
+
+
+def generate_delayed_sin_matrix(*, smp_num, tone_freq_n, mag, history_smp_num=0, noise_level=0):
+    assert np.abs(tone_freq_n) <= 1.0
+
+    t = np.arange(-history_smp_num, smp_num, 1, dtype=np.float64)
+    sig = mag * np.sin(2*np.pi * tone_freq_n * t, dtype=np.float64)
+    sig += noise_level*np.random.randn(len(sig))
+
+    # S is a matrix containing the delayed versions of 'sig': S[k,:] is 'sig' delayed by k taps
+    S = np.lib.stride_tricks.sliding_window_view(sig, smp_num)
+    S = np.flipud(S)
+
+    return S
+
+
+def generate_sin_signal(*, smp_num, tone_freq_n, mag, noise_level=0):
+    S = generate_delayed_sin_matrix(smp_num=smp_num, tone_freq_n=tone_freq_n, mag=mag, noise_level=noise_level)
+    return S[0,:]
+
+
+def calc_spectrum(sig, fs=1):
+    S = 20*np.log10(np.abs(fft.rfft(sig)))
+    f = fft.rfftfreq(len(sig), d=1/fs)
+
+    return f, S
+
+
+def upsample_fft(signal, upsample_factor):
+    fft_signal = fft.rfft(signal)
+
+    zeros_to_insert = np.zeros((upsample_factor - 1) * len(fft_signal))
+    fft_upsampled = np.concatenate((fft_signal, zeros_to_insert))
+    upsampled_signal = fft.irfft(fft_upsampled, n=upsample_factor*len(signal)) * upsample_factor
+    return upsampled_signal
+
+
+def downsample_fft(signal, downsample_factor):
+    fft_signal = fft.rfft(signal)
+    fft_downsampled = fft_signal[:(len(fft_signal) // downsample_factor)]
+    downsampled_signal = fft.irfft(fft_downsampled, n=len(signal)//downsample_factor) / downsample_factor
+    return downsampled_signal
+
+
+def resample(factor):
+    """Decorator to upsample before and downsample after the function."""
+    def decorator(func):
+        def wrapper(signal, *args, **kwargs):
+            upsampled_signal = upsample_fft(signal, factor)
+            result = func(upsampled_signal, *args, **kwargs)
+            downsampled_result = downsample_fft(result, factor)
+
+            return downsampled_result
+        return wrapper
+    return decorator

scripts/aliasing_analysis/environment.yml

@@ -0,0 +1,12 @@
+name: black-bird
+channels:
+  - conda-forge
+  - defaults
+dependencies:
+  - python=3.11
+  - jupyter=1.0.0
+  - numpy=1.26.4
+  - plotly=5.19.0
+  - scipy=1.12.0
+  - mpmath=1.3.0
+  - pytorch::pytorch=2.2.1

scripts/aliasing_analysis/learning_antialiasing.ipynb

@@ -0,0 +1,158 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from plotly import graph_objs as go\n",
+    "\n",
+    "import torch\n",
+    "import torch.nn as nn\n",
+    "import torch.optim as optim\n",
+    "\n",
+    "import dsplib\n",
+    "import plotlib\n",
+    "\n",
+    "from models import LutWithMemory\n",
+    "from training import Trainer"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import importlib\n",
+    "importlib.reload(dsplib)\n",
+    "importlib.reload(plotlib)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "enforce_cpu_on_mac = True\n",
+    "\n",
+    "if torch.cuda.is_available():\n",
+    "    device = torch.device(\"cuda\")\n",
+    "elif (torch.backends.mps.is_available()) and (not enforce_cpu_on_mac):\n",
+    "    device = torch.device(\"mps\")\n",
+    "else:\n",
+    "    device = torch.device(\"cpu\")\n",
+    "\n",
+    "device"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fs_hz = 48e3\n",
+    "signal_duration_sec = 1\n",
+    "smp_num = dsplib.calc_smp_num(signal_duration_sec, fs_hz)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "min_x = -5\n",
+    "max_x = 5\n",
+    "\n",
+    "@dsplib.resample(factor=32)\n",
+    "def tanh_resample(x):\n",
+    "    return np.tanh(x)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "S = dsplib.generate_delayed_sin_matrix(smp_num=smp_num, tone_freq_n=15100/fs_hz, mag=0.5, history_smp_num=3, noise_level=1e-4)\n",
+    "\n",
+    "ref_sig = tanh_resample(S[0,:])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Initialize model, loss, and optimizer\n",
+    "memory_depth = 3\n",
+    "model = LutWithMemory(input_range=(min_x, max_x), bins_num=64, memory_depth=memory_depth).to(device)\n",
+    "optimizer = optim.Adam(model.parameters(), lr=0.01)\n",
+    "\n",
+    "# Prepare data\n",
+    "X = torch.from_numpy(S.copy().astype(np.float32)).to(device)\n",
+    "y = torch.from_numpy(ref_sig.astype(np.float32)).to(device)\n",
+    "# Train\n",
+    "trainer = Trainer(model, optimizer=optimizer, criterion=nn.MSELoss())\n",
+    "trainer.train(X, y, num_epochs=int(2e4))\n",
+    "\n",
+    "fig = go.Figure()\n",
+    "fig.add_trace(go.Scatter(y=np.log10(trainer.loss_history)))\n",
+    "fig.show()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sig = S[0,:]\n",
+    "f, S_ref = dsplib.calc_spectrum(ref_sig, fs=fs_hz)\n",
+    "_, S_out = dsplib.calc_spectrum(np.tanh(sig), fs=fs_hz)\n",
+    "_, S_model = dsplib.calc_spectrum(model(X).detach().squeeze().numpy(), fs=fs_hz)\n",
+    "\n",
+    "fig = go.Figure()\n",
+    "fig.add_trace(go.Scatter(x=f, y=S_ref[:-1], name=\"ref\"))\n",
+    "fig.add_trace(go.Scatter(x=f, y=S_out, line=dict(width=2, dash='dash'), name=\"vanilla\"))\n",
+    "fig.add_trace(go.Scatter(x=f, y=S_model, line=dict(width=2, dash='dashdot'), name=\"model\"))\n",
+    "\n",
+    "fig.update_layout(\n",
+    "    title=f\"Signal spectrum before and after tanh\",\n",
+    "    xaxis_title='Frequency, Hz',\n",
+    "    yaxis_title='Spectral density'\n",
+    ")\n",
+    "\n",
+    "fig.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "black-bird",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.8"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}