From a067634c1b9cc71dd3dd8ee69ecdfc7bcdca1c67 Mon Sep 17 00:00:00 2001 From: Dmitry Khrykin Date: Mon, 26 Feb 2024 21:12:55 +0200 Subject: [PATCH 01/13] Copy over LadderFilter implementation --- CMakeLists.txt | 1 + source/dsp/LadderFilter.cpp | 173 ++++++++++++++++++++++++++++++++++++ source/dsp/LadderFilter.h | 127 ++++++++++++++++++++++++++ source/dsp/Voice.h | 5 +- 4 files changed, 304 insertions(+), 2 deletions(-) create mode 100644 source/dsp/LadderFilter.cpp create mode 100644 source/dsp/LadderFilter.h 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/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(); } From d87f78c7855d9a846f1289f7ad2689d44d523dc5 Mon Sep 17 00:00:00 2001 From: Slava Primenko Date: Tue, 27 Feb 2024 00:34:49 +0200 Subject: [PATCH 02/13] Do quick-n-dirty aliasing analysis --- scripts/aliasing_analysis/README.md | 18 +++ scripts/aliasing_analysis/environment.yml | 11 ++ scripts/aliasing_analysis/nonlinearity.ipynb | 160 +++++++++++++++++++ 3 files changed, 189 insertions(+) create mode 100644 scripts/aliasing_analysis/README.md create mode 100644 scripts/aliasing_analysis/environment.yml create mode 100644 scripts/aliasing_analysis/nonlinearity.ipynb diff --git a/scripts/aliasing_analysis/README.md b/scripts/aliasing_analysis/README.md new file mode 100644 index 0000000..7226c98 --- /dev/null +++ b/scripts/aliasing_analysis/README.md @@ -0,0 +1,18 @@ +# Exploring LadderFilter aliasing issue + + +## 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 a recommended editor/IDE. diff --git a/scripts/aliasing_analysis/environment.yml b/scripts/aliasing_analysis/environment.yml new file mode 100644 index 0000000..4997943 --- /dev/null +++ b/scripts/aliasing_analysis/environment.yml @@ -0,0 +1,11 @@ +name: black-bird-test +channels: + - conda-forge + - defaults +dependencies: + - python=3.11 + - jupyter=1.0.0 + - numpy=1.26.4 + - plotly=5.19.0 + - scipy=1.12.0 + diff --git a/scripts/aliasing_analysis/nonlinearity.ipynb b/scripts/aliasing_analysis/nonlinearity.ipynb new file mode 100644 index 0000000..29b1849 --- /dev/null +++ b/scripts/aliasing_analysis/nonlinearity.ipynb @@ -0,0 +1,160 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from plotly import graph_objs as go\n", + "from scipy import fft" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fs_hz = 48e3\n", + "signal_duration_sec = 0.1\n", + "t = np.arange(0, signal_duration_sec, 1/fs_hz)\n", + "\n", + "print(f\"Number of samples: {len(t)}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def calc_spectrum(sig, fs=fs_hz):\n", + " S = 20*np.log10(np.abs(fft.rfft(sig)))\n", + " f = fft.rfftfreq(len(sig), d=1/fs)\n", + "\n", + " return f, S" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def saturationFunc(sig):\n", + " return np.tanh(sig)\n", + "\n", + "x = np.linspace(-3, 3)\n", + "fig = go.Figure()\n", + "fig.add_trace(go.Scatter(x=x, y=saturationFunc(x), name='saturationFunc'))\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "tone_freq_hz = 4.3e3\n", + "mag = 0.1\n", + "sig = mag * np.sin(2*np.pi * tone_freq_hz * t)\n", + "\n", + "f, S = calc_spectrum(sig)\n", + "_, nlS = calc_spectrum(saturationFunc(sig))\n", + "\n", + "fig = go.Figure()\n", + "fig.add_trace(go.Scatter(x=f, y=S, name='Before'))\n", + "fig.add_trace(go.Scatter(x=f, y=nlS, line=dict(width=2, dash='dash'), name='After'))\n", + "fig.update_layout(title='Nonlinearity impact',\n", + " xaxis_title='Frequency, Hz',\n", + " yaxis_title='Spectral density')\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mag = 0.1\n", + "\n", + "fig = go.Figure()\n", + "freq_points_hz = np.arange(300, fs_hz/2, 200)\n", + "for tone_freq_hz in freq_points_hz:\n", + " sig = mag * np.sin(2*np.pi * tone_freq_hz * t)\n", + " f, S = calc_spectrum(sig)\n", + " _, nlS = calc_spectrum(saturationFunc(sig))\n", + " \n", + " fig.add_trace(go.Scatter(x=f, y=S,\n", + " visible=False, line=dict(color='blue'),\n", + " name=\"Before\"))\n", + " \n", + " fig.add_trace(go.Scatter(x=f, y=nlS,\n", + " visible=False, line=dict(color='red', width=2, dash='dash'),\n", + " name=\"After\"))\n", + "\n", + "active_ind = 1\n", + "fig.data[2*active_ind].visible = True\n", + "fig.data[2*active_ind+1].visible = True\n", + "\n", + "steps = []\n", + "for i in range(len(freq_points_hz)):\n", + " step = dict(\n", + " method=\"update\",\n", + " args=[{\"visible\": [False] * len(fig.data)},\n", + " {\"title\": f\"Frequency: {freq_points_hz[i]} Hz\"}],\n", + " label=str(int(freq_points_hz[i]))\n", + " )\n", + " step[\"args\"][0][\"visible\"][2*i] = True\n", + " step[\"args\"][0][\"visible\"][2*i+1] = True\n", + " steps.append(step)\n", + "\n", + "sliders = [dict(\n", + " active=active_ind,\n", + " currentvalue={\"prefix\": \"Frequency: \", \"suffix\": \" Hz\"},\n", + " pad={\"t\": 50},\n", + " steps=steps\n", + ")]\n", + "\n", + "fig.update_layout(\n", + " sliders=sliders,\n", + " title=f\"Signal frequency: {freq_points_hz[active_ind]} Hz\"\n", + ")\n", + "\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From aca0dc80eb7e2a60e97f0d95a71d839a7200a96a Mon Sep 17 00:00:00 2001 From: Slava Primenko Date: Tue, 27 Feb 2024 01:17:50 +0200 Subject: [PATCH 03/13] Update aliasing_analysis/README --- scripts/aliasing_analysis/README.md | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/scripts/aliasing_analysis/README.md b/scripts/aliasing_analysis/README.md index 7226c98..c3eca6f 100644 --- a/scripts/aliasing_analysis/README.md +++ b/scripts/aliasing_analysis/README.md @@ -1,5 +1,8 @@ # Exploring LadderFilter aliasing issue +`nonlinearity.ipynb` demonstrates how `tanh(x)` function produces new +frequency components and how aliasing phenomenon looks like at various +frequencies and magnitudes. ## Configuring Python environment @@ -15,4 +18,4 @@ conda env create -f environment.yml conda activate black-bird ``` -VSCode is a recommended editor/IDE. +VSCode is the recommended editor/IDE. From 5729b0a83772c5b8fd248b28741be3357adb99de Mon Sep 17 00:00:00 2001 From: Slava Primenko Date: Mon, 11 Mar 2024 00:39:03 +0100 Subject: [PATCH 04/13] Try 'antiderivative antialiasing' method for sin wave MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The idea described in 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 was simulated for the case of tanh(x) nonlinearity and sin wave. --- .gitignore | 2 + scripts/aliasing_analysis/README.md | 2 +- scripts/aliasing_analysis/dsplib.py | 13 ++ scripts/aliasing_analysis/environment.yml | 1 + scripts/aliasing_analysis/nonlinearity.ipynb | 194 +++++++++++++----- .../aliasing_analysis/tanh_antiderivative.py | 76 +++++++ 6 files changed, 238 insertions(+), 50 deletions(-) create mode 100644 scripts/aliasing_analysis/dsplib.py create mode 100644 scripts/aliasing_analysis/tanh_antiderivative.py 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/scripts/aliasing_analysis/README.md b/scripts/aliasing_analysis/README.md index c3eca6f..d665fd5 100644 --- a/scripts/aliasing_analysis/README.md +++ b/scripts/aliasing_analysis/README.md @@ -1,7 +1,7 @@ # Exploring LadderFilter aliasing issue `nonlinearity.ipynb` demonstrates how `tanh(x)` function produces new -frequency components and how aliasing phenomenon looks like at various +frequency components and how the aliasing phenomenon looks like at various frequencies and magnitudes. ## Configuring Python environment diff --git a/scripts/aliasing_analysis/dsplib.py b/scripts/aliasing_analysis/dsplib.py new file mode 100644 index 0000000..0800b25 --- /dev/null +++ b/scripts/aliasing_analysis/dsplib.py @@ -0,0 +1,13 @@ +import numpy as np +from scipy import fft + + +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 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) diff --git a/scripts/aliasing_analysis/environment.yml b/scripts/aliasing_analysis/environment.yml index 4997943..543a184 100644 --- a/scripts/aliasing_analysis/environment.yml +++ b/scripts/aliasing_analysis/environment.yml @@ -8,4 +8,5 @@ dependencies: - numpy=1.26.4 - plotly=5.19.0 - scipy=1.12.0 + - mpmath=1.3.0 diff --git a/scripts/aliasing_analysis/nonlinearity.ipynb b/scripts/aliasing_analysis/nonlinearity.ipynb index 29b1849..9f0a4c0 100644 --- a/scripts/aliasing_analysis/nonlinearity.ipynb +++ b/scripts/aliasing_analysis/nonlinearity.ipynb @@ -8,7 +8,9 @@ "source": [ "import numpy as np\n", "from plotly import graph_objs as go\n", - "from scipy import fft" + "\n", + "import tanh_antiderivative\n", + "import dsplib" ] }, { @@ -18,23 +20,14 @@ "outputs": [], "source": [ "fs_hz = 48e3\n", - "signal_duration_sec = 0.1\n", - "t = np.arange(0, signal_duration_sec, 1/fs_hz)\n", - "\n", - "print(f\"Number of samples: {len(t)}\")" + "signal_duration_sec = 0.2" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "def calc_spectrum(sig, fs=fs_hz):\n", - " S = 20*np.log10(np.abs(fft.rfft(sig)))\n", - " f = fft.rfftfreq(len(sig), d=1/fs)\n", - "\n", - " return f, S" + "## The tanh(x) function" ] }, { @@ -43,71 +36,123 @@ "metadata": {}, "outputs": [], "source": [ - "def saturationFunc(sig):\n", - " return np.tanh(sig)\n", - "\n", "x = np.linspace(-3, 3)\n", "fig = go.Figure()\n", - "fig.add_trace(go.Scatter(x=x, y=saturationFunc(x), name='saturationFunc'))\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": "markdown", + "metadata": {}, + "source": [ + "## Sine wave" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Frequency impact" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "tone_freq_hz = 4.3e3\n", - "mag = 0.1\n", - "sig = mag * np.sin(2*np.pi * tone_freq_hz * t)\n", - "\n", - "f, S = calc_spectrum(sig)\n", - "_, nlS = calc_spectrum(saturationFunc(sig))\n", + "mag = 2\n", + "t = dsplib.generate_time_signal(signal_duration_sec=signal_duration_sec, fs_hz=fs_hz)\n", "\n", "fig = go.Figure()\n", - "fig.add_trace(go.Scatter(x=f, y=S, name='Before'))\n", - "fig.add_trace(go.Scatter(x=f, y=nlS, line=dict(width=2, dash='dash'), name='After'))\n", - "fig.update_layout(title='Nonlinearity impact',\n", - " xaxis_title='Frequency, Hz',\n", - " yaxis_title='Spectral density')\n", + "freq_points_hz = np.arange(300, fs_hz/2, 200)\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, line=dict(color='blue'), name=\"Before\"))\n", + " fig.add_trace(go.Scatter(x=f, y=S_out, visible=False, line=dict(color='red', width=2, dash='dash'), name=\"After\"))\n", + "\n", + "active_ind = 1\n", + "lines_num = 2\n", + "fig.data[lines_num*active_ind].visible = True\n", + "fig.data[lines_num*active_ind+1].visible = True\n", + "\n", + "steps = []\n", + "for i in range(len(freq_points_hz)):\n", + " step = dict(\n", + " method=\"update\",\n", + " args=[{\"visible\": [False] * len(fig.data)},],\n", + " label=str(int(freq_points_hz[i]))\n", + " )\n", + " step[\"args\"][0][\"visible\"][lines_num*i] = True\n", + " step[\"args\"][0][\"visible\"][lines_num*i+1] = True\n", + " steps.append(step)\n", + "\n", + "sliders = [dict(\n", + " active=active_ind,\n", + " currentvalue={\"prefix\": \"Frequency: \", \"suffix\": \" Hz\"},\n", + " pad={\"t\": 50},\n", + " steps=steps\n", + ")]\n", + "\n", + "fig.update_layout(\n", + " sliders=sliders,\n", + " title=f\"Signal spectrum before and after tanh\",\n", + " xaxis_title='Frequency, Hz',\n", + " yaxis_title='Spectral density'\n", + ")\n", + "fig.update_layout(hovermode=\"x unified\")\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": [ - "mag = 0.1\n", + "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", - "freq_points_hz = np.arange(300, fs_hz/2, 200)\n", - "for tone_freq_hz in freq_points_hz:\n", - " sig = mag * np.sin(2*np.pi * tone_freq_hz * t)\n", - " f, S = calc_spectrum(sig)\n", - " _, nlS = calc_spectrum(saturationFunc(sig))\n", - " \n", - " fig.add_trace(go.Scatter(x=f, y=S,\n", - " visible=False, line=dict(color='blue'),\n", - " name=\"Before\"))\n", - " \n", - " fig.add_trace(go.Scatter(x=f, y=nlS,\n", - " visible=False, line=dict(color='red', width=2, dash='dash'),\n", - " name=\"After\"))\n", + "mag_points = np.arange(0.1, 5, 0.1)\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, line=dict(color='blue'), name=\"Before\"))\n", + " fig.add_trace(go.Scatter(x=f, y=S_out, visible=False, line=dict(color='red', width=2, dash='dash'), name=\"After\"))\n", "\n", "active_ind = 1\n", "fig.data[2*active_ind].visible = True\n", "fig.data[2*active_ind+1].visible = True\n", "\n", "steps = []\n", - "for i in range(len(freq_points_hz)):\n", + "for i in range(len(mag_points)):\n", " step = dict(\n", " method=\"update\",\n", - " args=[{\"visible\": [False] * len(fig.data)},\n", - " {\"title\": f\"Frequency: {freq_points_hz[i]} Hz\"}],\n", - " label=str(int(freq_points_hz[i]))\n", + " args=[{\"visible\": [False] * len(fig.data)}],\n", + " label=f\"{mag_points[i]:.1f}\"\n", " )\n", " step[\"args\"][0][\"visible\"][2*i] = True\n", " step[\"args\"][0][\"visible\"][2*i+1] = True\n", @@ -115,25 +160,76 @@ "\n", "sliders = [dict(\n", " active=active_ind,\n", - " currentvalue={\"prefix\": \"Frequency: \", \"suffix\": \" Hz\"},\n", + " currentvalue={\"prefix\": \"Magnitude: \", \"suffix\": \"\"},\n", " pad={\"t\": 50},\n", " steps=steps\n", ")]\n", "\n", "fig.update_layout(\n", " sliders=sliders,\n", - " title=f\"Signal frequency: {freq_points_hz[active_ind]} Hz\"\n", + " title=f\"Signal spectrum before and after tanh\",\n", + " xaxis_title='Frequency, Hz',\n", + " yaxis_title='Spectral density'\n", ")\n", "\n", "fig.show()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Trying to avoid the aliasing" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "tone_freq_hz = 9100\n", + "mag = 0.5\n", + "history_samples = 3\n", + "\n", + "t = dsplib.generate_time_signal(signal_duration_sec=signal_duration_sec, fs_hz=fs_hz, history_smp_num=history_samples)\n", + "sig = mag * np.sin(2*np.pi * tone_freq_hz * t)\n", + "n = 5e-4*np.random.randn(len(t)) # noise\n", + "sig += n\n", + "\n", + "sig_len = len(sig) - history_samples\n", + "\n", + "# S is a matrix containing the delayed versions of 'sig': S[k,:] is 'sig' delayed by k taps\n", + "S = np.lib.stride_tricks.sliding_window_view(sig, sig_len)\n", + "S = np.flipud(S)" + ] + }, + { + "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=\"Before\"))\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": { 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) From d9ffcb84abfbb38b3a1717ea442847dea3675d3a Mon Sep 17 00:00:00 2001 From: Slava Primenko Date: Mon, 11 Mar 2024 21:07:26 +0100 Subject: [PATCH 05/13] Fix conda environment name --- scripts/aliasing_analysis/environment.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/aliasing_analysis/environment.yml b/scripts/aliasing_analysis/environment.yml index 543a184..538eb7c 100644 --- a/scripts/aliasing_analysis/environment.yml +++ b/scripts/aliasing_analysis/environment.yml @@ -1,4 +1,4 @@ -name: black-bird-test +name: black-bird channels: - conda-forge - defaults From 59049ca33cef862a3a6280f37ab0b61dec559655 Mon Sep 17 00:00:00 2001 From: Slava Primenko Date: Mon, 11 Mar 2024 23:55:15 +0100 Subject: [PATCH 06/13] Do refactoring --- .../aliasing_analysis/antiderivatives.ipynb | 95 +++++++ scripts/aliasing_analysis/dsplib.py | 31 ++- scripts/aliasing_analysis/environment.yml | 1 - scripts/aliasing_analysis/nonlinearity.ipynb | 256 ------------------ scripts/aliasing_analysis/plotlib.py | 29 ++ scripts/aliasing_analysis/tanh_impact.ipynb | 160 +++++++++++ 6 files changed, 311 insertions(+), 261 deletions(-) create mode 100644 scripts/aliasing_analysis/antiderivatives.ipynb delete mode 100644 scripts/aliasing_analysis/nonlinearity.ipynb create mode 100644 scripts/aliasing_analysis/plotlib.py create mode 100644 scripts/aliasing_analysis/tanh_impact.ipynb 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 index 0800b25..6c08649 100644 --- a/scripts/aliasing_analysis/dsplib.py +++ b/scripts/aliasing_analysis/dsplib.py @@ -2,12 +2,35 @@ 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) + + +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) + sig = mag * np.sin(2*np.pi * tone_freq_n * t) + 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 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) diff --git a/scripts/aliasing_analysis/environment.yml b/scripts/aliasing_analysis/environment.yml index 538eb7c..b8a7e89 100644 --- a/scripts/aliasing_analysis/environment.yml +++ b/scripts/aliasing_analysis/environment.yml @@ -9,4 +9,3 @@ dependencies: - plotly=5.19.0 - scipy=1.12.0 - mpmath=1.3.0 - diff --git a/scripts/aliasing_analysis/nonlinearity.ipynb b/scripts/aliasing_analysis/nonlinearity.ipynb deleted file mode 100644 index 9f0a4c0..0000000 --- a/scripts/aliasing_analysis/nonlinearity.ipynb +++ /dev/null @@ -1,256 +0,0 @@ -{ - "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": [ - "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": "markdown", - "metadata": {}, - "source": [ - "## Sine wave" - ] - }, - { - "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", - "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, line=dict(color='blue'), name=\"Before\"))\n", - " fig.add_trace(go.Scatter(x=f, y=S_out, visible=False, line=dict(color='red', width=2, dash='dash'), name=\"After\"))\n", - "\n", - "active_ind = 1\n", - "lines_num = 2\n", - "fig.data[lines_num*active_ind].visible = True\n", - "fig.data[lines_num*active_ind+1].visible = True\n", - "\n", - "steps = []\n", - "for i in range(len(freq_points_hz)):\n", - " step = dict(\n", - " method=\"update\",\n", - " args=[{\"visible\": [False] * len(fig.data)},],\n", - " label=str(int(freq_points_hz[i]))\n", - " )\n", - " step[\"args\"][0][\"visible\"][lines_num*i] = True\n", - " step[\"args\"][0][\"visible\"][lines_num*i+1] = True\n", - " steps.append(step)\n", - "\n", - "sliders = [dict(\n", - " active=active_ind,\n", - " currentvalue={\"prefix\": \"Frequency: \", \"suffix\": \" Hz\"},\n", - " pad={\"t\": 50},\n", - " steps=steps\n", - ")]\n", - "\n", - "fig.update_layout(\n", - " sliders=sliders,\n", - " title=f\"Signal spectrum before and after tanh\",\n", - " xaxis_title='Frequency, Hz',\n", - " yaxis_title='Spectral density'\n", - ")\n", - "fig.update_layout(hovermode=\"x unified\")\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", - "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, line=dict(color='blue'), name=\"Before\"))\n", - " fig.add_trace(go.Scatter(x=f, y=S_out, visible=False, line=dict(color='red', width=2, dash='dash'), name=\"After\"))\n", - "\n", - "active_ind = 1\n", - "fig.data[2*active_ind].visible = True\n", - "fig.data[2*active_ind+1].visible = True\n", - "\n", - "steps = []\n", - "for i in range(len(mag_points)):\n", - " step = dict(\n", - " method=\"update\",\n", - " args=[{\"visible\": [False] * len(fig.data)}],\n", - " label=f\"{mag_points[i]:.1f}\"\n", - " )\n", - " step[\"args\"][0][\"visible\"][2*i] = True\n", - " step[\"args\"][0][\"visible\"][2*i+1] = True\n", - " steps.append(step)\n", - "\n", - "sliders = [dict(\n", - " active=active_ind,\n", - " currentvalue={\"prefix\": \"Magnitude: \", \"suffix\": \"\"},\n", - " pad={\"t\": 50},\n", - " steps=steps\n", - ")]\n", - "\n", - "fig.update_layout(\n", - " sliders=sliders,\n", - " title=f\"Signal spectrum before and after tanh\",\n", - " xaxis_title='Frequency, Hz',\n", - " yaxis_title='Spectral density'\n", - ")\n", - "\n", - "fig.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Trying to avoid the aliasing" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tone_freq_hz = 9100\n", - "mag = 0.5\n", - "history_samples = 3\n", - "\n", - "t = dsplib.generate_time_signal(signal_duration_sec=signal_duration_sec, fs_hz=fs_hz, history_smp_num=history_samples)\n", - "sig = mag * np.sin(2*np.pi * tone_freq_hz * t)\n", - "n = 5e-4*np.random.randn(len(t)) # noise\n", - "sig += n\n", - "\n", - "sig_len = len(sig) - history_samples\n", - "\n", - "# S is a matrix containing the delayed versions of 'sig': S[k,:] is 'sig' delayed by k taps\n", - "S = np.lib.stride_tricks.sliding_window_view(sig, sig_len)\n", - "S = np.flipud(S)" - ] - }, - { - "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=\"Before\"))\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.7" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} 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_impact.ipynb b/scripts/aliasing_analysis/tanh_impact.ipynb new file mode 100644 index 0000000..f1c0e1c --- /dev/null +++ b/scripts/aliasing_analysis/tanh_impact.ipynb @@ -0,0 +1,160 @@ +{ + "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": "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 +} From 6166052e29a7c44c4f16a5e1f0167d8c959c8aab Mon Sep 17 00:00:00 2001 From: Slava Primenko Date: Mon, 11 Mar 2024 23:56:14 +0100 Subject: [PATCH 07/13] Add upsampling demo --- scripts/aliasing_analysis/dsplib.py | 16 +++ scripts/aliasing_analysis/upsampling.ipynb | 118 +++++++++++++++++++++ 2 files changed, 134 insertions(+) create mode 100644 scripts/aliasing_analysis/upsampling.ipynb diff --git a/scripts/aliasing_analysis/dsplib.py b/scripts/aliasing_analysis/dsplib.py index 6c08649..ba563a6 100644 --- a/scripts/aliasing_analysis/dsplib.py +++ b/scripts/aliasing_analysis/dsplib.py @@ -34,3 +34,19 @@ def calc_spectrum(sig, fs=1): 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 diff --git a/scripts/aliasing_analysis/upsampling.ipynb b/scripts/aliasing_analysis/upsampling.ipynb new file mode 100644 index 0000000..55c603e --- /dev/null +++ b/scripts/aliasing_analysis/upsampling.ipynb @@ -0,0 +1,118 @@ +{ + "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": [ + "mag = 3\n", + "t = dsplib.generate_time_signal(signal_duration_sec=signal_duration_sec, fs_hz=fs_hz)\n", + "ups_factor = 32\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", + " ups_sig = dsplib.upsample_fft(sig, ups_factor)\n", + " tanh_ups = np.tanh(ups_sig)\n", + " tanh_ups_dwn = dsplib.downsample_fft(tanh_ups, ups_factor)\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()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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 +} From 3ed077054a1de7e34d0dbdce3a3ad5d660c8153a Mon Sep 17 00:00:00 2001 From: Slava Primenko Date: Tue, 12 Mar 2024 00:02:42 +0100 Subject: [PATCH 08/13] Fix outdated README.md --- scripts/aliasing_analysis/README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/aliasing_analysis/README.md b/scripts/aliasing_analysis/README.md index d665fd5..798b6b6 100644 --- a/scripts/aliasing_analysis/README.md +++ b/scripts/aliasing_analysis/README.md @@ -1,6 +1,6 @@ # Exploring LadderFilter aliasing issue -`nonlinearity.ipynb` demonstrates how `tanh(x)` function produces new +`tanh_impact.ipynb` demonstrates how `tanh(x)` function produces new frequency components and how the aliasing phenomenon looks like at various frequencies and magnitudes. @@ -18,4 +18,4 @@ conda env create -f environment.yml conda activate black-bird ``` -VSCode is the recommended editor/IDE. +VSCode is the recommended editor/IDE. From a3657d0ff1466ea1f7b0e970b911f10971ea32ea Mon Sep 17 00:00:00 2001 From: Slava Primenko Date: Sun, 17 Mar 2024 00:33:03 +0100 Subject: [PATCH 09/13] Decorate resampling It is now possible to run 'upsample -> func -> downsample' chain using `@dsplib.resample(factor)` syntax. --- scripts/aliasing_analysis/dsplib.py | 19 +++++++++++++++--- scripts/aliasing_analysis/upsampling.ipynb | 23 +++++++++++----------- 2 files changed, 28 insertions(+), 14 deletions(-) diff --git a/scripts/aliasing_analysis/dsplib.py b/scripts/aliasing_analysis/dsplib.py index ba563a6..429c5b8 100644 --- a/scripts/aliasing_analysis/dsplib.py +++ b/scripts/aliasing_analysis/dsplib.py @@ -7,14 +7,14 @@ def calc_smp_num(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) + 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) - sig = mag * np.sin(2*np.pi * tone_freq_n * t) + 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 @@ -50,3 +50,16 @@ def downsample_fft(signal, downsample_factor): 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/upsampling.ipynb b/scripts/aliasing_analysis/upsampling.ipynb index 55c603e..9adc2b9 100644 --- a/scripts/aliasing_analysis/upsampling.ipynb +++ b/scripts/aliasing_analysis/upsampling.ipynb @@ -45,6 +45,17 @@ "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, @@ -53,7 +64,6 @@ "source": [ "mag = 3\n", "t = dsplib.generate_time_signal(signal_duration_sec=signal_duration_sec, fs_hz=fs_hz)\n", - "ups_factor = 32\n", "\n", "fig = go.Figure()\n", "freq_points_hz = np.arange(300, fs_hz/2, 200)\n", @@ -64,9 +74,7 @@ " f, S_in = dsplib.calc_spectrum(sig, fs=fs_hz)\n", " _, S_out = dsplib.calc_spectrum(np.tanh(sig), fs=fs_hz)\n", "\n", - " ups_sig = dsplib.upsample_fft(sig, ups_factor)\n", - " tanh_ups = np.tanh(ups_sig)\n", - " tanh_ups_dwn = dsplib.downsample_fft(tanh_ups, ups_factor)\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", @@ -85,13 +93,6 @@ "\n", "fig.show()" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { From d28428ff2ff7731549374305feaf2ddbb321226b Mon Sep 17 00:00:00 2001 From: Slava Primenko Date: Sun, 17 Mar 2024 00:46:52 +0100 Subject: [PATCH 10/13] Learn aliasing-free tanh for single frequency and magnitude --- .../learning_antialiasing.ipynb | 140 ++++++++++++++++++ scripts/aliasing_analysis/models.py | 29 ++++ scripts/aliasing_analysis/models_test.py | 48 ++++++ scripts/aliasing_analysis/training.py | 20 +++ 4 files changed, 237 insertions(+) create mode 100644 scripts/aliasing_analysis/learning_antialiasing.ipynb create mode 100644 scripts/aliasing_analysis/models.py create mode 100644 scripts/aliasing_analysis/models_test.py create mode 100644 scripts/aliasing_analysis/training.py diff --git a/scripts/aliasing_analysis/learning_antialiasing.ipynb b/scripts/aliasing_analysis/learning_antialiasing.ipynb new file mode 100644 index 0000000..f088045 --- /dev/null +++ b/scripts/aliasing_analysis/learning_antialiasing.ipynb @@ -0,0 +1,140 @@ +{ + "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": [ + "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=9100/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)\n", + "optimizer = optim.Adam(model.parameters(), lr=0.01)\n", + "\n", + "# Prepare data\n", + "X = torch.from_numpy(S.copy().astype(np.float32))\n", + "y = torch.from_numpy(ref_sig.astype(np.float32))\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..fdb1dab --- /dev/null +++ b/scripts/aliasing_analysis/models.py @@ -0,0 +1,29 @@ +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)) + + 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 \ No newline at end of file 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/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()) From 514cb25215c6620702f4f647c02b51909a6a8d25 Mon Sep 17 00:00:00 2001 From: Slava Primenko Date: Sun, 17 Mar 2024 01:30:51 +0100 Subject: [PATCH 11/13] Send data and model to(device) --- .../learning_antialiasing.ipynb | 26 ++++++++++++++++--- scripts/aliasing_analysis/models.py | 5 ++-- 2 files changed, 25 insertions(+), 6 deletions(-) diff --git a/scripts/aliasing_analysis/learning_antialiasing.ipynb b/scripts/aliasing_analysis/learning_antialiasing.ipynb index f088045..d0cc184 100644 --- a/scripts/aliasing_analysis/learning_antialiasing.ipynb +++ b/scripts/aliasing_analysis/learning_antialiasing.ipynb @@ -31,6 +31,24 @@ "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, @@ -62,7 +80,7 @@ "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", + "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,:])" ] @@ -75,12 +93,12 @@ "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)\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))\n", - "y = torch.from_numpy(ref_sig.astype(np.float32))\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", diff --git a/scripts/aliasing_analysis/models.py b/scripts/aliasing_analysis/models.py index fdb1dab..3380da4 100644 --- a/scripts/aliasing_analysis/models.py +++ b/scripts/aliasing_analysis/models.py @@ -15,7 +15,7 @@ def forward(self, X): 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)) + 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,:]) @@ -26,4 +26,5 @@ def forward(self, X): y_pred += lut_out * X[mem_ind,:].unsqueeze(-1) - return y_pred \ No newline at end of file + return y_pred + From 4ee11ad308b04d98ef4d8da79b676fed45ca2b86 Mon Sep 17 00:00:00 2001 From: Slava Primenko Date: Sun, 17 Mar 2024 01:45:05 +0100 Subject: [PATCH 12/13] Add pytorch to conda environment.yml --- scripts/aliasing_analysis/environment.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/scripts/aliasing_analysis/environment.yml b/scripts/aliasing_analysis/environment.yml index b8a7e89..71a2797 100644 --- a/scripts/aliasing_analysis/environment.yml +++ b/scripts/aliasing_analysis/environment.yml @@ -9,3 +9,4 @@ dependencies: - plotly=5.19.0 - scipy=1.12.0 - mpmath=1.3.0 + - pytorch::pytorch=2.2.1 From d58a4ce04433ac57b86cd20a6c529e84ed1a7514 Mon Sep 17 00:00:00 2001 From: Slava Primenko Date: Wed, 3 Jul 2024 22:48:48 +0200 Subject: [PATCH 13/13] Look into tanh approximations --- scripts/aliasing_analysis/tanh_impact.ipynb | 182 ++++++++++++++++++++ 1 file changed, 182 insertions(+) diff --git a/scripts/aliasing_analysis/tanh_impact.ipynb b/scripts/aliasing_analysis/tanh_impact.ipynb index f1c0e1c..160cfc3 100644 --- a/scripts/aliasing_analysis/tanh_impact.ipynb +++ b/scripts/aliasing_analysis/tanh_impact.ipynb @@ -43,6 +43,188 @@ "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": {},