diff --git a/.gitignore b/.gitignore index 25f8c4e..a8bf03b 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,5 @@ docs/.jekyll-cache .cache .DS_Store + +__pycache__ diff --git a/CMakeLists.txt b/CMakeLists.txt index 31faa17..51a8612 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/scripts/aliasing_analysis/README.md b/scripts/aliasing_analysis/README.md new file mode 100644 index 0000000..798b6b6 --- /dev/null +++ b/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. diff --git a/scripts/aliasing_analysis/antiderivatives.ipynb b/scripts/aliasing_analysis/antiderivatives.ipynb new file mode 100644 index 0000000..8ad6812 --- /dev/null +++ b/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 +} diff --git a/scripts/aliasing_analysis/dsplib.py b/scripts/aliasing_analysis/dsplib.py new file mode 100644 index 0000000..429c5b8 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/scripts/aliasing_analysis/environment.yml b/scripts/aliasing_analysis/environment.yml new file mode 100644 index 0000000..71a2797 --- /dev/null +++ b/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 diff --git a/scripts/aliasing_analysis/learning_antialiasing.ipynb b/scripts/aliasing_analysis/learning_antialiasing.ipynb new file mode 100644 index 0000000..d0cc184 --- /dev/null +++ b/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 +} diff --git a/scripts/aliasing_analysis/models.py b/scripts/aliasing_analysis/models.py new file mode 100644 index 0000000..3380da4 --- /dev/null +++ b/scripts/aliasing_analysis/models.py @@ -0,0 +1,30 @@ +import torch +import torch.nn as nn + +class LutWithMemory(nn.Module): + def __init__(self, input_range, bins_num=10, memory_depth=3): + super(LutWithMemory, self).__init__() + self.embeddings = nn.ModuleList([nn.Embedding(num_embeddings=bins_num, embedding_dim=1) for _ in range(1+memory_depth)]) + self.bins_num = bins_num + self.memory_depth = memory_depth + self.min_x = input_range[0] + self.max_x = input_range[1] + + def forward(self, X): + X_normalized = (X - self.min_x) / (self.max_x - self.min_x) * (self.bins_num - 1) + lower_indices = X_normalized.floor().long() + upper_indices = torch.clamp(lower_indices + 1, max=self.bins_num-1) + + y_pred = torch.zeros((X_normalized.shape[1], 1), device=X_normalized.device) + + for mem_ind in range(1+self.memory_depth): + lower_embeddings = self.embeddings[mem_ind](lower_indices[mem_ind,:]) + upper_embeddings = self.embeddings[mem_ind](upper_indices[mem_ind,:]) + + weights = (X_normalized[mem_ind,:] - lower_indices[mem_ind,:].float()).unsqueeze(-1) + lut_out = lower_embeddings + weights * (upper_embeddings - lower_embeddings) + + y_pred += lut_out * X[mem_ind,:].unsqueeze(-1) + + return y_pred + diff --git a/scripts/aliasing_analysis/models_test.py b/scripts/aliasing_analysis/models_test.py new file mode 100644 index 0000000..3e1ab79 --- /dev/null +++ b/scripts/aliasing_analysis/models_test.py @@ -0,0 +1,48 @@ +import torch +import torch.nn as nn +import torch.optim as optim +import numpy as np + +from models import LutWithMemory +from training import Trainer + +import unittest + +class TestLutWithMemoryFit(unittest.TestCase): + def test_tanh_fit(self): + """Test that LutWithMemory can fit tanh(x) on np.arange data""" + + min_x = -5 + max_x = 5 + memory_depth = 3 + + model = LutWithMemory(input_range=(min_x, max_x), bins_num=64, memory_depth=memory_depth) + criterion = nn.MSELoss() + optimizer = optim.Adam(model.parameters(), lr=0.01) + + + x = torch.arange(min_x, max_x, 0.01) + y = np.tanh(x) + + X = x.repeat(memory_depth+1, 1) + + # Train + trainer = Trainer(model, optimizer=optimizer, criterion=criterion) + trainer.train(X, y, num_epochs=3000) + + # Generate predictions from the model for plotting + x_eval = torch.arange(min_x, max_x, 0.1) + X = x_eval.repeat(memory_depth+1, 1) + y_pred_eval = model(X).detach().squeeze() + + # Convert to numpy for plotting + x_eval_np = x_eval.numpy() + y_eval_np = np.tanh(x_eval_np) + y_pred_eval_np = y_pred_eval.numpy() + + tol = 1e-3 + self.assertTrue(np.allclose(y_eval_np, y_pred_eval_np, atol=tol), + msg=f"The prediction error is greater than the threshold {tol}.") + +if __name__ == '__main__': + unittest.main() diff --git a/scripts/aliasing_analysis/plotlib.py b/scripts/aliasing_analysis/plotlib.py new file mode 100644 index 0000000..b6342d7 --- /dev/null +++ b/scripts/aliasing_analysis/plotlib.py @@ -0,0 +1,29 @@ +def apply_sliders(fig, lines_num, labels, prefix="Prefix: ", suffix=" suffix"): + active_ind = 1 + + for k in range(lines_num): + fig.data[lines_num*active_ind+k].visible = True + + steps = [] + for i in range(len(fig.data)//lines_num): + step = dict( + method="update", + args=[{"visible": [False] * len(fig.data)},], + label=str(labels[i]) + ) + + for k in range(lines_num): + step["args"][0]["visible"][lines_num*i+k] = True + + steps.append(step) + + sliders = [dict( + active=active_ind, + currentvalue={"prefix": prefix, "suffix": suffix}, + pad={"t": 50}, + steps=steps + )] + + fig.update_layout(sliders=sliders, hovermode="x unified") + + return fig diff --git a/scripts/aliasing_analysis/tanh_antiderivative.py b/scripts/aliasing_analysis/tanh_antiderivative.py new file mode 100644 index 0000000..fb84ce9 --- /dev/null +++ b/scripts/aliasing_analysis/tanh_antiderivative.py @@ -0,0 +1,76 @@ +""" +Applying the antiderivative antialiasing idea to tanh(x) function. + +Stefan Bilbao, Fabián Esqueda, Julian D. Parker, Vesa Välimäki. +Antiderivative antialiasing for memoryless nonlinearities. +IEEE Signal Processing Letters, Vol. 24, No. 7, pp. 1049-1053, July 2017. +DOI: 10.1109/LSP.2017.2675541 +""" + +import numpy as np +import mpmath + + +def order1(X): + assert X.shape[0] >= 2 + + def D(z0, z1): + return (_F1(z0) - _F1(z1)) / (z0 - z1) + + # X[k, :] contains a signal delayed by k taps + return D(X[0,:], X[1,:]) + + +def order2(X): + assert X.shape[0] >= 3 + + def D(z0, z1): + return (_F2(z0) - _F2(z1)) / (z0 - z1) + + # X[k, :] contains a signal delayed by k taps + A = 2 / (X[0,:] - X[2,:]) + D_01 = D(X[0,:], X[1,:]) + D_12 = D(X[1,:], X[2,:]) + + return A * (D_01 - D_12) + + +def order3(X): + assert X.shape[0] >= 4 + + def D(z0, z1): + return (_F3(z0) - _F3(z1)) / (z0 - z1) + + # X[k, :] contains a signal delayed by k taps + A = 1 / (X[1,:] - X[2,:]) + B = 2 / (X[0,:] - X[2,:]) + C = 2 / (X[1,:] - X[3,:]) + + D_01 = D(X[0,:], X[1,:]) + D_12 = D(X[1,:], X[2,:]) + D_23 = D(X[2,:], X[3,:]) + + return A * (B * (D_01 - D_12) - C * (D_12 - D_23)) + +def _Li_n(x, n): + np_polylog = np.frompyfunc(lambda z: mpmath.polylog(n, z), 1, 1) + return np_polylog(x) + + +def _F1(x): + return np.log(np.cosh(x)) + + +def _F2(x): + """ Calculated with Wolfram Alpha """ + exp = np.exp(-2*x) + C1 = 0 # This won't matter as C1 coeffs cancel each other in 'order2' equation + return C1*x + 0.5*_Li_n(-exp, 2) - 0.5*x**2 - x*np.log(1 + exp) + x*np.log(np.cosh(x)) + + +def _F3(x): + """ Calculated with Wolfram Alpha """ + exp = np.exp(-2*x) + C1 = 0 + C2 = 0 + return (1/6)*x*(3*C1*x + 6*C2 - 2*x**2 - 3*x*np.log(exp+1) + 3*x*np.log(np.cosh(x))) - (1/4)*_Li_n(-exp, 3) diff --git a/scripts/aliasing_analysis/tanh_impact.ipynb b/scripts/aliasing_analysis/tanh_impact.ipynb new file mode 100644 index 0000000..160cfc3 --- /dev/null +++ b/scripts/aliasing_analysis/tanh_impact.ipynb @@ -0,0 +1,342 @@ +{ + "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 dsplib\n", + "import plotlib" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fs_hz = 48e3\n", + "signal_duration_sec = 0.2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### The tanh(x) function" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x = np.linspace(-3, 3)\n", + "fig = go.Figure()\n", + "fig.add_trace(go.Scatter(x=x, y=np.tanh(x), name='tanh'))\n", + "fig.update_layout(xaxis_title='Input magnitude', yaxis_title='Output magnitude')\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def tanh_approx(x):\n", + " numerator = x**7 + 378*x**5 + 17325*x**3 + 135135*x\n", + " denominator = 28*x**6 + 3150*x**4 + 62370*x**2 + 135135\n", + "\n", + " return numerator / denominator" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = go.Figure()\n", + "fig.add_trace(go.Scatter(x=x, y=np.tanh(x), name='tanh'))\n", + "fig.add_trace(go.Scatter(x=x, y=tanh_approx(x), name='tanh_approx'))\n", + "fig.add_trace(go.Scatter(x=x, y=np.tanh(x)-tanh_approx(x), name='error'))\n", + "fig.update_layout(xaxis_title='Input magnitude', yaxis_title='Output magnitude')\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import sympy as sp\n", + "from sympy import fps, sin, exp" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x = sp.symbols('x')\n", + "fps(sin(1 + x)).truncate()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define the variable\n", + "x = sp.symbols('x')\n", + "\n", + "# Define the numerator and denominator\n", + "numerator = x**7 + 378*x**5 + 17325*x**3 + 135135*x\n", + "denominator = 28*x**6 + 3150*x**4 + 62370*x**2 + 135135\n", + "\n", + "# Perform the polynomial long division\n", + "q, r = sp.div(numerator, denominator, domain='QQ')\n", + "\n", + "# Display the results\n", + "print(q, r, sep='\\n')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from sympy import symbols\n", + "from sympy import simplify, lambdify\n", + "from sympy import symbols, tanh, series\n", + "\n", + "\n", + "def calc_tanh_approximation(order=6):\n", + " x = symbols('x')\n", + "\n", + " add = 2*order - 1\n", + " denom = 2*order+1\n", + "\n", + " for _ in range(order):\n", + " denom = add + x**2 / denom\n", + " add -= 2\n", + "\n", + " F = x / denom\n", + " simplified_F = simplify(F)\n", + " return simplified_F\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def calc_taylor(order=6, x0=0):\n", + " # Define the variable x\n", + " x = symbols('x')\n", + " f = tanh(x)\n", + " taylor_series = series(f, x, x0, order).removeO()\n", + "\n", + " return taylor_series" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "a = np.linspace(-3, 3)\n", + "\n", + "fig = go.Figure()\n", + "labels = []\n", + "errors = []\n", + "\n", + "for order in range(11):\n", + " F = calc_tanh_approximation(order=order)\n", + " F_numeric = lambdify(x, F)\n", + " error = np.tanh(a)-F_numeric(a)\n", + " fig.add_trace(go.Scatter(x=a, y=error, visible=False, name='error'))\n", + " labels.append(order)\n", + " errors.append(np.max(np.abs(error)))\n", + "\n", + "fig = plotlib.apply_sliders(fig, lines_num=1, labels=labels, prefix=\"Order: \", suffix=\"\")\n", + "\n", + "fig.update_layout(\n", + " title=f\"Tanh vs Padé\",\n", + " xaxis_title='Input',\n", + " yaxis_title='Output',\n", + ")\n", + "fig.show()\n", + "\n", + "fig_e = go.Figure()\n", + "fig_e.add_trace(go.Scatter(x=labels, y=np.log10(errors)))\n", + "fig_e.update_layout(xaxis_title='Padé order', yaxis_title='Max error, log10')\n", + "fig_e.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "a = np.linspace(-3, 3)\n", + "\n", + "fig = go.Figure()\n", + "labels = []\n", + "for x0 in np.linspace(-3, 3, num=7, endpoint=True):\n", + " F = calc_taylor(order=4, x0=x0)\n", + " F_numeric = lambdify(x, F)\n", + "\n", + " fig.add_trace(go.Scatter(x=a, y=np.tanh(a), visible=False, name='tanh'))\n", + " fig.add_trace(go.Scatter(x=a, y=F_numeric(a), visible=False, name='Taylor[tanh]'))\n", + "\n", + " labels.append(x0)\n", + "\n", + "fig = plotlib.apply_sliders(fig, lines_num=2, labels=labels, prefix=\"Taylor center: \", suffix=\"\")\n", + "\n", + "fig.update_layout(\n", + " title=f\"Tanh vs Taylor\",\n", + " xaxis_title='Input',\n", + " yaxis_title='Output',\n", + " yaxis_range=(-1.5, 1.5)\n", + ")\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "F" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will do the analysis for a sin wave, as it makes it easy to interpret frequency domain data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Frequency impact" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mag = 2\n", + "t = dsplib.generate_time_signal(signal_duration_sec=signal_duration_sec, fs_hz=fs_hz)\n", + "\n", + "fig = go.Figure()\n", + "freq_points_hz = np.arange(300, fs_hz/2, 200)\n", + "labels = []\n", + "for tone_freq_hz in freq_points_hz:\n", + " sig = mag * np.sin(2*np.pi * tone_freq_hz * t)\n", + " f, S_in = dsplib.calc_spectrum(sig, fs=fs_hz)\n", + " _, S_out = dsplib.calc_spectrum(np.tanh(sig), fs=fs_hz)\n", + "\n", + " fig.add_trace(go.Scatter(x=f, y=S_in, visible=False, name=\"input\"))\n", + " fig.add_trace(go.Scatter(x=f, y=S_out, visible=False, line=dict(width=2, dash='dash'), name=\"tanh\"))\n", + " labels.append(int(tone_freq_hz))\n", + "\n", + "fig = plotlib.apply_sliders(fig, lines_num=2, labels=labels, prefix=\"Frequency: \", suffix=\" Hz\")\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", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can see that as the tone frequency increases, the generated by the `tanh` harmonics end up behind the `fs/2` frequency which leads to aliasing." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Magnitude impact" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "tone_freq_hz = 900\n", + "t = dsplib.generate_time_signal(signal_duration_sec=signal_duration_sec, fs_hz=fs_hz, history_smp_num=0)\n", + "sig = np.sin(2*np.pi * tone_freq_hz * t)\n", + "f, S_in_0 = dsplib.calc_spectrum(sig, fs=fs_hz)\n", + "\n", + "fig = go.Figure()\n", + "mag_points = np.arange(0.1, 5, 0.1)\n", + "labels = []\n", + "for mag in mag_points:\n", + " S_in = S_in_0 + 20*np.log10(mag)\n", + " _, S_out = dsplib.calc_spectrum(np.tanh(mag*sig), fs=fs_hz)\n", + "\n", + " fig.add_trace(go.Scatter(x=f, y=S_in, visible=False, name=\"input\"))\n", + " fig.add_trace(go.Scatter(x=f, y=S_out, visible=False, line=dict(width=2, dash='dash'), name=\"tanh\"))\n", + " labels.append(f\"{mag:.1f}\")\n", + "\n", + "fig = plotlib.apply_sliders(fig, lines_num=2, labels=labels, prefix=\"Magnitude: \", suffix=\"\")\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 +} diff --git a/scripts/aliasing_analysis/training.py b/scripts/aliasing_analysis/training.py new file mode 100644 index 0000000..848c10e --- /dev/null +++ b/scripts/aliasing_analysis/training.py @@ -0,0 +1,20 @@ +import torch.nn as nn + +class Trainer: + def __init__(self, model, *, optimizer, criterion=nn.MSELoss()): + self.model = model + self.criterion = criterion + self.optimizer = optimizer + + self.loss_history = [] + + def train(self, X, y, num_epochs=2000): + for epoch in range(num_epochs): + self.optimizer.zero_grad() + y_pred = self.model(X) + loss = self.criterion(y_pred, y.unsqueeze(-1)) + loss.backward() + self.optimizer.step() + + if epoch % 10 == 0: + self.loss_history.append(loss.item()) diff --git a/scripts/aliasing_analysis/upsampling.ipynb b/scripts/aliasing_analysis/upsampling.ipynb new file mode 100644 index 0000000..9adc2b9 --- /dev/null +++ b/scripts/aliasing_analysis/upsampling.ipynb @@ -0,0 +1,119 @@ +{ + "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 dsplib\n", + "import plotlib" + ] + }, + { + "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": [ + "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)\n", + "sig = S[0,:]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "@dsplib.resample(factor=32)\n", + "def tanh_resample(x):\n", + " return np.tanh(x)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mag = 3\n", + "t = dsplib.generate_time_signal(signal_duration_sec=signal_duration_sec, fs_hz=fs_hz)\n", + "\n", + "fig = go.Figure()\n", + "freq_points_hz = np.arange(300, fs_hz/2, 200)\n", + "\n", + "labels = []\n", + "for tone_freq_hz in freq_points_hz:\n", + " sig = mag * np.sin(2*np.pi * tone_freq_hz * t)\n", + " f, S_in = dsplib.calc_spectrum(sig, fs=fs_hz)\n", + " _, S_out = dsplib.calc_spectrum(np.tanh(sig), fs=fs_hz)\n", + "\n", + " tanh_ups_dwn = tanh_resample(sig)\n", + " _, S_out_ups = dsplib.calc_spectrum(tanh_ups_dwn, fs=fs_hz)\n", + "\n", + " fig.add_trace(go.Scatter(x=f, y=S_in, visible=False, name=\"input\"))\n", + " fig.add_trace(go.Scatter(x=f, y=S_out, visible=False, line=dict(width=2, dash='dash'), name=\"tanh\"))\n", + " fig.add_trace(go.Scatter(x=f, y=S_out_ups, visible=False, line=dict(width=2, dash='dot'), name=\"tanh ups\"))\n", + "\n", + " labels.append(int(tone_freq_hz))\n", + "\n", + "fig = plotlib.apply_sliders(fig, lines_num=3, labels=labels, prefix=\"Frequency: \", suffix=\" Hz\")\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 +} diff --git a/source/dsp/LadderFilter.cpp b/source/dsp/LadderFilter.cpp new file mode 100644 index 0000000..bbb3d0d --- /dev/null +++ b/source/dsp/LadderFilter.cpp @@ -0,0 +1,173 @@ +#include "LadderFilter.h" + +namespace blackbird::dsp { + +//============================================================================== +template +LadderFilter::LadderFilter() : state(2) { + setSampleRate( + SampleType(1000)); // intentionally setting unrealistic default + // sample rate to catch missing initialisation bugs + setResonance(SampleType(0)); + setDrive(SampleType(1.2)); + + mode = Mode::LPF24; + setMode(Mode::LPF12); +} + +//============================================================================== +template +void LadderFilter::setMode(Mode newMode) noexcept { + if (newMode == mode) + return; + + switch (newMode) { + case Mode::LPF12: + A = {{SampleType(0), SampleType(0), SampleType(1), SampleType(0), + SampleType(0)}}; + comp = SampleType(0.5); + break; + case Mode::HPF12: + A = {{SampleType(1), SampleType(-2), SampleType(1), SampleType(0), + SampleType(0)}}; + comp = SampleType(0); + break; + case Mode::BPF12: + A = {{SampleType(0), SampleType(0), SampleType(-1), SampleType(1), + SampleType(0)}}; + comp = SampleType(0.5); + break; + case Mode::LPF24: + A = {{SampleType(0), SampleType(0), SampleType(0), SampleType(0), + SampleType(1)}}; + comp = SampleType(0.5); + break; + case Mode::HPF24: + A = {{SampleType(1), SampleType(-4), SampleType(6), SampleType(-4), + SampleType(1)}}; + comp = SampleType(0); + break; + case Mode::BPF24: + A = {{SampleType(0), SampleType(0), SampleType(1), SampleType(-2), + SampleType(1)}}; + comp = SampleType(0.5); + break; + default: + jassertfalse; + break; + } + + static constexpr auto outputGain = SampleType(1.2); + + for (auto &a : A) + a *= outputGain; + + mode = newMode; + reset(); +} + +//============================================================================== +template +void LadderFilter::prepare(const ProcessSpec &spec) { + setSampleRate(SampleType(spec.sampleRate)); + setNumChannels(spec.numChannels); + reset(); +} + +//============================================================================== +template void LadderFilter::reset() noexcept { + for (auto &s : state) + s.fill(SampleType(0)); + + cutoffTransformSmoother.setCurrentAndTargetValue( + cutoffTransformSmoother.getTargetValue()); + scaledResonanceSmoother.setCurrentAndTargetValue( + scaledResonanceSmoother.getTargetValue()); +} + +//============================================================================== +template +void LadderFilter::setCutoffFrequencyHz( + SampleType newCutoff) noexcept { + jassert(newCutoff > SampleType(0)); + cutoffFreqHz = newCutoff; + updateCutoffFreq(); +} + +//============================================================================== +template +void LadderFilter::setResonance(SampleType newResonance) noexcept { + jassert(newResonance >= SampleType(0) && newResonance <= SampleType(1)); + resonance = newResonance; + updateResonance(); +} + +//============================================================================== +template +void LadderFilter::setDrive(SampleType newDrive) noexcept { + jassert(newDrive >= SampleType(1)); + + drive = newDrive; + gain = std::pow(drive, SampleType(-2.642)) * SampleType(0.6103) + + SampleType(0.3903); + drive2 = drive * SampleType(0.04) + SampleType(0.96); + gain2 = std::pow(drive2, SampleType(-2.642)) * SampleType(0.6103) + + SampleType(0.3903); +} + +//============================================================================== +template +SampleType +LadderFilter::processSample(SampleType inputValue, + size_t channelToUse) noexcept { + auto &s = state[channelToUse]; + + const auto a1 = cutoffTransformValue; + const auto g = a1 * SampleType(-1) + SampleType(1); + const auto b0 = g * SampleType(0.76923076923); + const auto b1 = g * SampleType(0.23076923076); + + const auto dx = gain * saturationLUT(drive * inputValue); + const auto a = dx + scaledResonanceValue * SampleType(-4) * + (gain2 * saturationLUT(drive2 * s[4]) - dx * comp); + + const auto b = b1 * s[0] + a1 * s[1] + b0 * a; + const auto c = b1 * s[1] + a1 * s[2] + b0 * b; + const auto d = b1 * s[2] + a1 * s[3] + b0 * c; + const auto e = b1 * s[3] + a1 * s[4] + b0 * d; + + s[0] = a; + s[1] = b; + s[2] = c; + s[3] = d; + s[4] = e; + + return a * A[0] + b * A[1] + c * A[2] + d * A[3] + e * A[4]; +} + +//============================================================================== +template +void LadderFilter::updateSmoothers() noexcept { + cutoffTransformValue = cutoffTransformSmoother.getNextValue(); + scaledResonanceValue = scaledResonanceSmoother.getNextValue(); +} + +//============================================================================== +template +void LadderFilter::setSampleRate(SampleType newValue) noexcept { + jassert(newValue > SampleType(0)); + cutoffFreqScaler = + SampleType(-2.0 * juce::MathConstants::pi) / newValue; + + static constexpr SampleType smootherRampTimeSec = SampleType(0.05); + cutoffTransformSmoother.reset(newValue, smootherRampTimeSec); + scaledResonanceSmoother.reset(newValue, smootherRampTimeSec); + + updateCutoffFreq(); +} + +//============================================================================== +template class LadderFilter; +template class LadderFilter; + +} // namespace blackbird::dsp diff --git a/source/dsp/LadderFilter.h b/source/dsp/LadderFilter.h new file mode 100644 index 0000000..e098b0a --- /dev/null +++ b/source/dsp/LadderFilter.h @@ -0,0 +1,127 @@ +#pragma once + +#include + +namespace blackbird::dsp { + +using namespace juce; +using namespace juce::dsp; + +/** + Multi-mode filter based on the Moog ladder filter. + + @tags{DSP} +*/ +template class LadderFilter { +public: + //============================================================================== + using Mode = LadderFilterMode; + + //============================================================================== + /** Creates an uninitialised filter. Call prepare() before first use. */ + LadderFilter(); + + /** Enables or disables the filter. If disabled it will simply pass through + * the input signal. */ + void setEnabled(bool isEnabled) noexcept { enabled = isEnabled; } + + /** Sets filter mode. */ + void setMode(Mode newMode) noexcept; + + /** Initialises the filter. */ + void prepare(const ProcessSpec &spec); + + /** Returns the current number of channels. */ + size_t getNumChannels() const noexcept { return state.size(); } + + /** Resets the internal state variables of the filter. */ + void reset() noexcept; + + /** Sets the cutoff frequency of the filter. + + @param newCutoff cutoff frequency in Hz + */ + void setCutoffFrequencyHz(SampleType newCutoff) noexcept; + + /** Sets the resonance of the filter. + + @param newResonance a value between 0 and 1; higher values increase the + resonance and can result in self oscillation! + */ + void setResonance(SampleType newResonance) noexcept; + + /** Sets the amount of saturation in the filter. + + @param newDrive saturation amount; it can be any number greater than or + equal to one. Higher values result in more distortion. + */ + void setDrive(SampleType newDrive) noexcept; + + //============================================================================== + template + void process(const ProcessContext &context) noexcept { + const auto &inputBlock = context.getInputBlock(); + auto &outputBlock = context.getOutputBlock(); + const auto numChannels = outputBlock.getNumChannels(); + const auto numSamples = outputBlock.getNumSamples(); + + jassert(inputBlock.getNumChannels() <= getNumChannels()); + jassert(inputBlock.getNumChannels() == numChannels); + jassert(inputBlock.getNumSamples() == numSamples); + + if (!enabled || context.isBypassed) { + outputBlock.copyFrom(inputBlock); + return; + } + + for (size_t n = 0; n < numSamples; ++n) { + updateSmoothers(); + + for (size_t ch = 0; ch < numChannels; ++ch) + outputBlock.getChannelPointer(ch)[n] = + processSample(inputBlock.getChannelPointer(ch)[n], ch); + } + } + +protected: + //============================================================================== + SampleType processSample(SampleType inputValue, size_t channelToUse) noexcept; + void updateSmoothers() noexcept; + +private: + //============================================================================== + void setSampleRate(SampleType newValue) noexcept; + void setNumChannels(size_t newValue) { state.resize(newValue); } + void updateCutoffFreq() noexcept { + cutoffTransformSmoother.setTargetValue( + std::exp(cutoffFreqHz * cutoffFreqScaler)); + } + void updateResonance() noexcept { + scaledResonanceSmoother.setTargetValue( + jmap(resonance, SampleType(0.1), SampleType(1.0))); + } + + //============================================================================== + SampleType drive, drive2, gain, gain2, comp; + + static constexpr size_t numStates = 5; + std::vector> state; + std::array A; + + SmoothedValue cutoffTransformSmoother, scaledResonanceSmoother; + SampleType cutoffTransformValue, scaledResonanceValue; + + LookupTableTransform saturationLUT{ + [](SampleType x) { return std::tanh(x); }, SampleType(-5), SampleType(5), + 128}; + + SampleType cutoffFreqHz{SampleType(200)}; + SampleType resonance; + + SampleType cutoffFreqScaler; + + Mode mode; + bool enabled = true; +}; + +} // namespace blackbird::dsp diff --git a/source/dsp/Voice.h b/source/dsp/Voice.h index e352bcd..ac0c984 100644 --- a/source/dsp/Voice.h +++ b/source/dsp/Voice.h @@ -11,6 +11,7 @@ #pragma once #include "DSPParameters.h" +#include "LadderFilter.h" #include "LookupTablesBank.h" #include "VCAOscillator.h" #include @@ -222,7 +223,7 @@ class Voice : public SynthesiserVoice, private Timer { }; dsp::ProcessorChain, VCAOscillator, - dsp::LadderFilter, dsp::Gain> + blackbird::dsp::LadderFilter, dsp::Gain> processorChain; dsp::Oscillator lfo; @@ -241,7 +242,7 @@ class Voice : public SynthesiserVoice, private Timer { dsp::Gain &gainProcessor() { return processorChain.get(); } - dsp::LadderFilter &filter() { + blackbird::dsp::LadderFilter &filter() { return processorChain.get(); }