From 8ae3bd52494864fb3e84d53782147279c81ae35f Mon Sep 17 00:00:00 2001 From: Ivan COHEN Date: Thu, 5 Jul 2018 11:34:07 +0200 Subject: [PATCH 1/8] Improved VCF HP output and added proposal for a BP output I have replaced the previous Fundamental VCF HP output formula with the correct one, and I have added a commented instruction which could be used to provide a BP output as well if wanted. Note : I have not tested yet the submitted change, but these formulas are the one I use all the time in my own Moog-like / transistor ladder filter algorithms, such as the one in Sonic Academy ANA 2. For more information about the formula I have used, see there : http://www.kvraudio.com/forum/viewtopic.php?t=207647 --- src/VCF.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/VCF.cpp b/src/VCF.cpp index 195d8b4..49619f6 100644 --- a/src/VCF.cpp +++ b/src/VCF.cpp @@ -140,7 +140,8 @@ void VCF::step() { // Set outputs outputs[LPF_OUTPUT].value = 5.0f * filter.state[3]; - outputs[HPF_OUTPUT].value = 5.0f * (input - filter.state[3]); + outputs[HPF_OUTPUT].value = 5.0f * (input - 4*filter.state[0] + 6*filter.state[1] - 4*filter.state[2] + filter.state[3]); + //outputs[BPF_OUTPUT].value = 5.0f * (filter.state[1] - 2*filter.state[2] + filter.state[3]); } From 7443734c23d6c0e79aeeedc54992babf2020f236 Mon Sep 17 00:00:00 2001 From: Ivan COHEN Date: Thu, 5 Jul 2018 13:48:57 +0200 Subject: [PATCH 2/8] Replaced original VCF algorithm with Teemu Voipio's aka mystran's one This new algorithm for a transistor ladder filter is stable for all the cutoff frequencies, uses the TPT structure which ensures by design that any kind of cutoff frequency modulation will not create any artefact. And it uses a trick which makes it 0df on the linear side only, so that it is less CPU consuming than standard 0df filters using implicit solvers, but with its frequency response being the one expected. --- src/VCF.cpp | 293 ++++++++++++++++++++++++++-------------------------- 1 file changed, 146 insertions(+), 147 deletions(-) diff --git a/src/VCF.cpp b/src/VCF.cpp index 49619f6..dbbb0ef 100644 --- a/src/VCF.cpp +++ b/src/VCF.cpp @@ -1,175 +1,174 @@ -/* -The filter DSP code has been derived from -Miller Puckette's code hosted at -https://github.com/ddiakopoulos/MoogLadders/blob/master/src/RKSimulationModel.h -which is licensed for use under the following terms (MIT license): - - -Copyright (c) 2015, Miller Puckette. All rights reserved. - -Redistribution and use in source and binary forms, with or without -modification, are permitted provided that the following conditions are met: - -* Redistributions of source code must retain the above copyright notice, this - list of conditions and the following disclaimer. -* Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimer in the documentation - and/or other materials provided with the distribution. - -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE -FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -*/ +// LICENSE TERMS: Copyright 2012 Teemu Voipio +// +// You can use this however you like for pretty much any purpose, +// as long as you don't claim you wrote it. There is no warranty. +// +// Distribution of substantial portions of this code in source form +// must include this copyright notice and list of conditions. +// +// Improved by Ivan COHEN from musical entropy +// More about this algorithm here : +// http://www.kvraudio.com/forum/viewtopic.php?t=349859 +// http://www.kvraudio.com/forum/viewtopic.php?t=207647 #include "Fundamental.hpp" // The clipping function of a transistor pair is approximately tanh(x) -// TODO: Put this in a lookup table. 5th order approx doesn't seem to cut it +// This one is a Pade-approx for tanh(sqrt(x))/sqrt(x) inline float clip(float x) { - return tanhf(x); + float a = x*x; + return ((a + 105)*a + 945) / ((15*a + 420)*a + 945); } struct LadderFilter { - float cutoff = 1000.0f; - float resonance = 1.0f; - float state[4] = {}; - - void calculateDerivatives(float input, float *dstate, const float *state) { - float cutoff2Pi = 2*M_PI * cutoff; - - float satstate0 = clip(state[0]); - float satstate1 = clip(state[1]); - float satstate2 = clip(state[2]); - - dstate[0] = cutoff2Pi * (clip(input - resonance * state[3]) - satstate0); - dstate[1] = cutoff2Pi * (satstate0 - satstate1); - dstate[2] = cutoff2Pi * (satstate1 - satstate2); - dstate[3] = cutoff2Pi * (satstate2 - clip(state[3])); - } - - void process(float input, float dt) { - float deriv1[4], deriv2[4], deriv3[4], deriv4[4], tempState[4]; - - calculateDerivatives(input, deriv1, state); - for (int i = 0; i < 4; i++) - tempState[i] = state[i] + 0.5f * dt * deriv1[i]; - - calculateDerivatives(input, deriv2, tempState); - for (int i = 0; i < 4; i++) - tempState[i] = state[i] + 0.5f * dt * deriv2[i]; - - calculateDerivatives(input, deriv3, tempState); - for (int i = 0; i < 4; i++) - tempState[i] = state[i] + dt * deriv3[i]; - - calculateDerivatives(input, deriv4, tempState); - for (int i = 0; i < 4; i++) - state[i] += (1.0f / 6.0f) * dt * (deriv1[i] + 2.0f * deriv2[i] + 2.0f * deriv3[i] + deriv4[i]); - } - void reset() { - for (int i = 0; i < 4; i++) { - state[i] = 0.0f; - } - } + float g = 0.1f; + float resonance = 0.5f; + float state[4] = {}; + float zi = 0.f; + float output[3] = {}; + + void process(float input) { + // input with half delay, for non-linearities + const float ih = 0.5f * (input + zi); + + // evaluate the non-linear gains + const float t0 = g * clip(ih - resonance * state[3]); + const float t1 = g * clip(state[0]); + const float t2 = g * clip(state[1]); + const float t3 = g * clip(state[2]); + const float t4 = g * clip(state[3]); + + // update last LP1 output + const float t2t3 = t2*t3; + + float y3 = (s[3]*(1+t3) + s[2]*t3)*(1+t2); + y3 = (y3 + t2t3*s[1])*(1+t1); + y3 = (y3 + t1*t2t3*(s[0]+t0*input)); + y3 = y3 / ((1+t1)*(1+t2)*(1+t3)*(1+t4) + resonance*t0*t1*t2t3); + + // update other LP1 outputs + const float xx = t0 * (input - resonance * y3); + const float y0 = t1 * (s[0] + xx) / (1+t1); + const float y1 = t2 * (s[1] + y0) / (1+t2); + const float y2 = t3 * (s[2] + y1) / (1+t3); + + // update states + s[0] += 2 * (xx - y0); + s[1] += 2 * (y0 - y1); + s[2] += 2 * (y1 - y2); + s[3] += 2 * (y2 - t4*y3); + + // returns LP, HP and BP outputs + float y1t2 = y1/t2; + float y2t3 = y2/t3; + + output[0] = y3; + output[1] = xx/t0 - 4*y0/t1 + 6*y1t2 - 4*y2t3 + y3; + output[2] = y1t2 - 2*y2t3 + y3; + + // update delay input state + zi = input; + } + + void reset() { + for (int i = 0; i < 4; i++) { + state[i] = 0.0f; + } + zi = 0.f; + } }; struct VCF : Module { - enum ParamIds { - FREQ_PARAM, - FINE_PARAM, - RES_PARAM, - FREQ_CV_PARAM, - DRIVE_PARAM, - NUM_PARAMS - }; - enum InputIds { - FREQ_INPUT, - RES_INPUT, - DRIVE_INPUT, - IN_INPUT, - NUM_INPUTS - }; - enum OutputIds { - LPF_OUTPUT, - HPF_OUTPUT, - NUM_OUTPUTS - }; - - LadderFilter filter; - - VCF() : Module(NUM_PARAMS, NUM_INPUTS, NUM_OUTPUTS) {} - void step() override; - void onReset() override { - filter.reset(); - } + enum ParamIds { + FREQ_PARAM, + FINE_PARAM, + RES_PARAM, + FREQ_CV_PARAM, + DRIVE_PARAM, + NUM_PARAMS + }; + enum InputIds { + FREQ_INPUT, + RES_INPUT, + DRIVE_INPUT, + IN_INPUT, + NUM_INPUTS + }; + enum OutputIds { + LPF_OUTPUT, + HPF_OUTPUT, + NUM_OUTPUTS + }; + + LadderFilter filter; + + VCF() : Module(NUM_PARAMS, NUM_INPUTS, NUM_OUTPUTS) {} + void step() override; + void onReset() override { + filter.reset(); + } }; void VCF::step() { - float input = inputs[IN_INPUT].value / 5.0f; - float drive = params[DRIVE_PARAM].value + inputs[DRIVE_INPUT].value / 10.0f; - float gain = powf(100.0f, drive); - input *= gain; - // Add -60dB noise to bootstrap self-oscillation - input += 1e-6f * (2.0f*randomUniform() - 1.0f); - - // Set resonance - float res = params[RES_PARAM].value + inputs[RES_INPUT].value / 5.0f; - res = 5.5f * clamp(res, 0.0f, 1.0f); - filter.resonance = res; - - // Set cutoff frequency - float cutoffExp = params[FREQ_PARAM].value + params[FREQ_CV_PARAM].value * inputs[FREQ_INPUT].value / 5.0f; - cutoffExp = clamp(cutoffExp, 0.0f, 1.0f); - const float minCutoff = 15.0f; - const float maxCutoff = 8400.0f; - filter.cutoff = minCutoff * powf(maxCutoff / minCutoff, cutoffExp); - - // Push a sample to the state filter - filter.process(input, 1.0f/engineGetSampleRate()); - - // Set outputs - outputs[LPF_OUTPUT].value = 5.0f * filter.state[3]; - outputs[HPF_OUTPUT].value = 5.0f * (input - 4*filter.state[0] + 6*filter.state[1] - 4*filter.state[2] + filter.state[3]); - //outputs[BPF_OUTPUT].value = 5.0f * (filter.state[1] - 2*filter.state[2] + filter.state[3]); + float input = inputs[IN_INPUT].value / 5.0f; + float drive = params[DRIVE_PARAM].value + inputs[DRIVE_INPUT].value / 10.0f; + float gain = powf(100.0f, drive); + input *= gain; + // Add -60dB noise to bootstrap self-oscillation + input += 1e-6f * (2.0f*randomUniform() - 1.0f); + + // Set resonance + float res = params[RES_PARAM].value + inputs[RES_INPUT].value / 5.0f; + res = clamp(res, 0.0f, 1.0f); // resonance must be between 0 and 1 + filter.resonance = res; + + // Set cutoff frequency + float cutoffExp = params[FREQ_PARAM].value + params[FREQ_CV_PARAM].value * inputs[FREQ_INPUT].value / 5.0f; + cutoffExp = clamp(cutoffExp, 0.0f, 1.0f); + const float minCutoff = 15.0f; + const float maxCutoff = 20000.0f; + float cutoff = minCutoff * powf(maxCutoff / minCutoff, cutoffExp); + filter.g = tanf(float_Pi * cutoff / engineGetSampleRate()); + + // Push a sample to the state filter + filter.process(input); + + // Set outputs + outputs[LPF_OUTPUT].value = 5.0f * filter.output[0]; + outputs[HPF_OUTPUT].value = 5.0f * filter.output[1]; + //outputs[BPF_OUTPUT].value = 5.0f * filter.output[2]; } struct VCFWidget : ModuleWidget { - VCFWidget(VCF *module); + VCFWidget(VCF *module); }; VCFWidget::VCFWidget(VCF *module) : ModuleWidget(module) { - setPanel(SVG::load(assetPlugin(plugin, "res/VCF.svg"))); - - addChild(Widget::create(Vec(15, 0))); - addChild(Widget::create(Vec(box.size.x-30, 0))); - addChild(Widget::create(Vec(15, 365))); - addChild(Widget::create(Vec(box.size.x-30, 365))); - - addParam(ParamWidget::create(Vec(33, 61), module, VCF::FREQ_PARAM, 0.0f, 1.0f, 0.5f)); - addParam(ParamWidget::create(Vec(12, 143), module, VCF::FINE_PARAM, 0.0f, 1.0f, 0.5f)); - addParam(ParamWidget::create(Vec(71, 143), module, VCF::RES_PARAM, 0.0f, 1.0f, 0.0f)); - addParam(ParamWidget::create(Vec(12, 208), module, VCF::FREQ_CV_PARAM, -1.0f, 1.0f, 0.0f)); - addParam(ParamWidget::create(Vec(71, 208), module, VCF::DRIVE_PARAM, 0.0f, 1.0f, 0.0f)); - - addInput(Port::create(Vec(10, 276), Port::INPUT, module, VCF::FREQ_INPUT)); - addInput(Port::create(Vec(48, 276), Port::INPUT, module, VCF::RES_INPUT)); - addInput(Port::create(Vec(85, 276), Port::INPUT, module, VCF::DRIVE_INPUT)); - addInput(Port::create(Vec(10, 320), Port::INPUT, module, VCF::IN_INPUT)); - - addOutput(Port::create(Vec(48, 320), Port::OUTPUT, module, VCF::LPF_OUTPUT)); - addOutput(Port::create(Vec(85, 320), Port::OUTPUT, module, VCF::HPF_OUTPUT)); + setPanel(SVG::load(assetPlugin(plugin, "res/VCF.svg"))); + + addChild(Widget::create(Vec(15, 0))); + addChild(Widget::create(Vec(box.size.x-30, 0))); + addChild(Widget::create(Vec(15, 365))); + addChild(Widget::create(Vec(box.size.x-30, 365))); + + addParam(ParamWidget::create(Vec(33, 61), module, VCF::FREQ_PARAM, 0.0f, 1.0f, 0.5f)); + addParam(ParamWidget::create(Vec(12, 143), module, VCF::FINE_PARAM, 0.0f, 1.0f, 0.5f)); + addParam(ParamWidget::create(Vec(71, 143), module, VCF::RES_PARAM, 0.0f, 1.0f, 0.0f)); + addParam(ParamWidget::create(Vec(12, 208), module, VCF::FREQ_CV_PARAM, -1.0f, 1.0f, 0.0f)); + addParam(ParamWidget::create(Vec(71, 208), module, VCF::DRIVE_PARAM, 0.0f, 1.0f, 0.0f)); + + addInput(Port::create(Vec(10, 276), Port::INPUT, module, VCF::FREQ_INPUT)); + addInput(Port::create(Vec(48, 276), Port::INPUT, module, VCF::RES_INPUT)); + addInput(Port::create(Vec(85, 276), Port::INPUT, module, VCF::DRIVE_INPUT)); + addInput(Port::create(Vec(10, 320), Port::INPUT, module, VCF::IN_INPUT)); + + addOutput(Port::create(Vec(48, 320), Port::OUTPUT, module, VCF::LPF_OUTPUT)); + addOutput(Port::create(Vec(85, 320), Port::OUTPUT, module, VCF::HPF_OUTPUT)); } From 7932e0cbb7152a0124e22343fcc401aebeff2367 Mon Sep 17 00:00:00 2001 From: Ivan COHEN Date: Thu, 5 Jul 2018 13:52:58 +0200 Subject: [PATCH 3/8] Update indentation --- src/VCF.cpp | 261 ++++++++++++++++++++++++++-------------------------- 1 file changed, 128 insertions(+), 133 deletions(-) diff --git a/src/VCF.cpp b/src/VCF.cpp index dbbb0ef..754aa0e 100644 --- a/src/VCF.cpp +++ b/src/VCF.cpp @@ -18,157 +18,152 @@ // The clipping function of a transistor pair is approximately tanh(x) // This one is a Pade-approx for tanh(sqrt(x))/sqrt(x) inline float clip(float x) { - float a = x*x; - return ((a + 105)*a + 945) / ((15*a + 420)*a + 945); + float a = x*x; + return ((a + 105)*a + 945) / ((15*a + 420)*a + 945); } struct LadderFilter { - float g = 0.1f; - float resonance = 0.5f; - float state[4] = {}; - float zi = 0.f; - float output[3] = {}; - - void process(float input) { - // input with half delay, for non-linearities - const float ih = 0.5f * (input + zi); - - // evaluate the non-linear gains - const float t0 = g * clip(ih - resonance * state[3]); - const float t1 = g * clip(state[0]); - const float t2 = g * clip(state[1]); - const float t3 = g * clip(state[2]); - const float t4 = g * clip(state[3]); - - // update last LP1 output - const float t2t3 = t2*t3; - - float y3 = (s[3]*(1+t3) + s[2]*t3)*(1+t2); - y3 = (y3 + t2t3*s[1])*(1+t1); - y3 = (y3 + t1*t2t3*(s[0]+t0*input)); - y3 = y3 / ((1+t1)*(1+t2)*(1+t3)*(1+t4) + resonance*t0*t1*t2t3); - - // update other LP1 outputs - const float xx = t0 * (input - resonance * y3); - const float y0 = t1 * (s[0] + xx) / (1+t1); - const float y1 = t2 * (s[1] + y0) / (1+t2); - const float y2 = t3 * (s[2] + y1) / (1+t3); - - // update states - s[0] += 2 * (xx - y0); - s[1] += 2 * (y0 - y1); - s[2] += 2 * (y1 - y2); - s[3] += 2 * (y2 - t4*y3); - - // returns LP, HP and BP outputs - float y1t2 = y1/t2; - float y2t3 = y2/t3; - - output[0] = y3; - output[1] = xx/t0 - 4*y0/t1 + 6*y1t2 - 4*y2t3 + y3; - output[2] = y1t2 - 2*y2t3 + y3; - - // update delay input state - zi = input; - } - - void reset() { - for (int i = 0; i < 4; i++) { - state[i] = 0.0f; - } - zi = 0.f; - } + float g = 0.1f; + float resonance = 0.5f; + float state[4] = {}; + float zi = 0.f; + float output[3] = {}; + + void process(float input) { + // input with half delay, for non-linearities + const float ih = 0.5f * (input + zi); + + // evaluate the non-linear gains + const float t0 = g * clip(ih - resonance * state[3]); + const float t1 = g * clip(state[0]); + const float t2 = g * clip(state[1]); + const float t3 = g * clip(state[2]); + const float t4 = g * clip(state[3]); + + // update last LP1 output + float y3 = (s[3]*(1+t3) + s[2]*t3)*(1+t2); + y3 = (y3 + t2*t3*s[1])*(1+t1); + y3 = (y3 + t1*t2*t3*(s[0]+t0*input)); + y3 = y3 / ((1+t1)*(1+t2)*(1+t3)*(1+t4) + resonance*t0*t1*t2*t3); + + // update other LP1 outputs + const float xx = t0 * (input - resonance * y3); + const float y0 = t1 * (s[0] + xx) / (1+t1); + const float y1 = t2 * (s[1] + y0) / (1+t2); + const float y2 = t3 * (s[2] + y1) / (1+t3); + + // update states + s[0] += 2 * (xx - y0); + s[1] += 2 * (y0 - y1); + s[2] += 2 * (y1 - y2); + s[3] += 2 * (y2 - t4*y3); + + // returns LP, HP and BP outputs + output[0] = y3; + output[1] = xx/t0 - 4*y0/t1 + 6*y1/t2 - 4*y2/t3 + y3; + output[2] = y1/t2 - 2*y2/t3 + y3; + + // update delay input state + zi = input; + } + + void reset() { + for (int i = 0; i < 4; i++) { + state[i] = 0.0f; + } + zi = 0.f; + } }; struct VCF : Module { - enum ParamIds { - FREQ_PARAM, - FINE_PARAM, - RES_PARAM, - FREQ_CV_PARAM, - DRIVE_PARAM, - NUM_PARAMS - }; - enum InputIds { - FREQ_INPUT, - RES_INPUT, - DRIVE_INPUT, - IN_INPUT, - NUM_INPUTS - }; - enum OutputIds { - LPF_OUTPUT, - HPF_OUTPUT, - NUM_OUTPUTS - }; - - LadderFilter filter; - - VCF() : Module(NUM_PARAMS, NUM_INPUTS, NUM_OUTPUTS) {} - void step() override; - void onReset() override { - filter.reset(); - } + enum ParamIds { + FREQ_PARAM, + FINE_PARAM, + RES_PARAM, + FREQ_CV_PARAM, + DRIVE_PARAM, + NUM_PARAMS + }; + enum InputIds { + FREQ_INPUT, + RES_INPUT, + DRIVE_INPUT, + IN_INPUT, + NUM_INPUTS + }; + enum OutputIds { + LPF_OUTPUT, + HPF_OUTPUT, + NUM_OUTPUTS + }; + + LadderFilter filter; + + VCF() : Module(NUM_PARAMS, NUM_INPUTS, NUM_OUTPUTS) {} + void step() override; + void onReset() override { + filter.reset(); + } }; void VCF::step() { - float input = inputs[IN_INPUT].value / 5.0f; - float drive = params[DRIVE_PARAM].value + inputs[DRIVE_INPUT].value / 10.0f; - float gain = powf(100.0f, drive); - input *= gain; - // Add -60dB noise to bootstrap self-oscillation - input += 1e-6f * (2.0f*randomUniform() - 1.0f); - - // Set resonance - float res = params[RES_PARAM].value + inputs[RES_INPUT].value / 5.0f; - res = clamp(res, 0.0f, 1.0f); // resonance must be between 0 and 1 - filter.resonance = res; - - // Set cutoff frequency - float cutoffExp = params[FREQ_PARAM].value + params[FREQ_CV_PARAM].value * inputs[FREQ_INPUT].value / 5.0f; - cutoffExp = clamp(cutoffExp, 0.0f, 1.0f); - const float minCutoff = 15.0f; - const float maxCutoff = 20000.0f; - float cutoff = minCutoff * powf(maxCutoff / minCutoff, cutoffExp); - filter.g = tanf(float_Pi * cutoff / engineGetSampleRate()); - - // Push a sample to the state filter - filter.process(input); - - // Set outputs - outputs[LPF_OUTPUT].value = 5.0f * filter.output[0]; - outputs[HPF_OUTPUT].value = 5.0f * filter.output[1]; - //outputs[BPF_OUTPUT].value = 5.0f * filter.output[2]; + float input = inputs[IN_INPUT].value / 5.0f; + float drive = params[DRIVE_PARAM].value + inputs[DRIVE_INPUT].value / 10.0f; + float gain = powf(100.0f, drive); + input *= gain; + // Add -60dB noise to bootstrap self-oscillation + input += 1e-6f * (2.0f*randomUniform() - 1.0f); + + // Set resonance + float res = params[RES_PARAM].value + inputs[RES_INPUT].value / 5.0f; + res = clamp(res, 0.0f, 1.0f); // resonance must be between 0 and 1 + filter.resonance = res; + + // Set cutoff frequency + float cutoffExp = params[FREQ_PARAM].value + params[FREQ_CV_PARAM].value * inputs[FREQ_INPUT].value / 5.0f; + cutoffExp = clamp(cutoffExp, 0.0f, 1.0f); + const float minCutoff = 15.0f; + const float maxCutoff = 20000.0f; + float cutoff = minCutoff * powf(maxCutoff / minCutoff, cutoffExp); + filter.g = tanf(float_Pi * cutoff / engineGetSampleRate()); + + // Push a sample to the state filter + filter.process(input); + + // Set outputs + outputs[LPF_OUTPUT].value = 5.0f * filter.output[0]; + outputs[HPF_OUTPUT].value = 5.0f * filter.output[1]; + //outputs[BPF_OUTPUT].value = 5.0f * filter.output[2]; } struct VCFWidget : ModuleWidget { - VCFWidget(VCF *module); + VCFWidget(VCF *module); }; VCFWidget::VCFWidget(VCF *module) : ModuleWidget(module) { - setPanel(SVG::load(assetPlugin(plugin, "res/VCF.svg"))); - - addChild(Widget::create(Vec(15, 0))); - addChild(Widget::create(Vec(box.size.x-30, 0))); - addChild(Widget::create(Vec(15, 365))); - addChild(Widget::create(Vec(box.size.x-30, 365))); - - addParam(ParamWidget::create(Vec(33, 61), module, VCF::FREQ_PARAM, 0.0f, 1.0f, 0.5f)); - addParam(ParamWidget::create(Vec(12, 143), module, VCF::FINE_PARAM, 0.0f, 1.0f, 0.5f)); - addParam(ParamWidget::create(Vec(71, 143), module, VCF::RES_PARAM, 0.0f, 1.0f, 0.0f)); - addParam(ParamWidget::create(Vec(12, 208), module, VCF::FREQ_CV_PARAM, -1.0f, 1.0f, 0.0f)); - addParam(ParamWidget::create(Vec(71, 208), module, VCF::DRIVE_PARAM, 0.0f, 1.0f, 0.0f)); - - addInput(Port::create(Vec(10, 276), Port::INPUT, module, VCF::FREQ_INPUT)); - addInput(Port::create(Vec(48, 276), Port::INPUT, module, VCF::RES_INPUT)); - addInput(Port::create(Vec(85, 276), Port::INPUT, module, VCF::DRIVE_INPUT)); - addInput(Port::create(Vec(10, 320), Port::INPUT, module, VCF::IN_INPUT)); - - addOutput(Port::create(Vec(48, 320), Port::OUTPUT, module, VCF::LPF_OUTPUT)); - addOutput(Port::create(Vec(85, 320), Port::OUTPUT, module, VCF::HPF_OUTPUT)); + setPanel(SVG::load(assetPlugin(plugin, "res/VCF.svg"))); + + addChild(Widget::create(Vec(15, 0))); + addChild(Widget::create(Vec(box.size.x-30, 0))); + addChild(Widget::create(Vec(15, 365))); + addChild(Widget::create(Vec(box.size.x-30, 365))); + + addParam(ParamWidget::create(Vec(33, 61), module, VCF::FREQ_PARAM, 0.0f, 1.0f, 0.5f)); + addParam(ParamWidget::create(Vec(12, 143), module, VCF::FINE_PARAM, 0.0f, 1.0f, 0.5f)); + addParam(ParamWidget::create(Vec(71, 143), module, VCF::RES_PARAM, 0.0f, 1.0f, 0.0f)); + addParam(ParamWidget::create(Vec(12, 208), module, VCF::FREQ_CV_PARAM, -1.0f, 1.0f, 0.0f)); + addParam(ParamWidget::create(Vec(71, 208), module, VCF::DRIVE_PARAM, 0.0f, 1.0f, 0.0f)); + + addInput(Port::create(Vec(10, 276), Port::INPUT, module, VCF::FREQ_INPUT)); + addInput(Port::create(Vec(48, 276), Port::INPUT, module, VCF::RES_INPUT)); + addInput(Port::create(Vec(85, 276), Port::INPUT, module, VCF::DRIVE_INPUT)); + addInput(Port::create(Vec(10, 320), Port::INPUT, module, VCF::IN_INPUT)); + + addOutput(Port::create(Vec(48, 320), Port::OUTPUT, module, VCF::LPF_OUTPUT)); + addOutput(Port::create(Vec(85, 320), Port::OUTPUT, module, VCF::HPF_OUTPUT)); } From 226880f4ae6b60aab3fd5a22f4b74af5aa6c1f05 Mon Sep 17 00:00:00 2001 From: Ivan COHEN Date: Thu, 5 Jul 2018 14:36:36 +0200 Subject: [PATCH 4/8] Optimized a little VCF algorithm --- src/VCF.cpp | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/VCF.cpp b/src/VCF.cpp index 754aa0e..d01d61f 100644 --- a/src/VCF.cpp +++ b/src/VCF.cpp @@ -41,10 +41,12 @@ struct LadderFilter { const float t4 = g * clip(state[3]); // update last LP1 output + const float t2t3 = t2*t3; + float y3 = (s[3]*(1+t3) + s[2]*t3)*(1+t2); - y3 = (y3 + t2*t3*s[1])*(1+t1); - y3 = (y3 + t1*t2*t3*(s[0]+t0*input)); - y3 = y3 / ((1+t1)*(1+t2)*(1+t3)*(1+t4) + resonance*t0*t1*t2*t3); + y3 = (y3 + t2t3*s[1])*(1+t1); + y3 = (y3 + t1*t2t3*(s[0]+t0*input)); + y3 = y3 / ((1+t1)*(1+t2)*(1+t3)*(1+t4) + resonance*t0*t1*t2t3); // update other LP1 outputs const float xx = t0 * (input - resonance * y3); @@ -59,9 +61,12 @@ struct LadderFilter { s[3] += 2 * (y2 - t4*y3); // returns LP, HP and BP outputs + const float y1t2 = y1 / t2; + const float y2t3 = y2 / t3; + output[0] = y3; - output[1] = xx/t0 - 4*y0/t1 + 6*y1/t2 - 4*y2/t3 + y3; - output[2] = y1/t2 - 2*y2/t3 + y3; + output[1] = xx/t0 - 4*y0/t1 + 6*y1t2 - 4*y2t3 + y3; + output[2] = y1t2 - 2*y2t3 + y3; // update delay input state zi = input; From a2a9fc48d0e8384c58a3e4b034785b4a4b0c46a7 Mon Sep 17 00:00:00 2001 From: Andrew Belt Date: Thu, 5 Jul 2018 13:22:44 -0400 Subject: [PATCH 5/8] Fix build, update parameter scaling --- src/VCF.cpp | 78 +++++++++++++++++++++++++++-------------------------- 1 file changed, 40 insertions(+), 38 deletions(-) diff --git a/src/VCF.cpp b/src/VCF.cpp index d01d61f..cef12a1 100644 --- a/src/VCF.cpp +++ b/src/VCF.cpp @@ -13,18 +13,20 @@ #include "Fundamental.hpp" +#include "dsp/functions.hpp" // The clipping function of a transistor pair is approximately tanh(x) // This one is a Pade-approx for tanh(sqrt(x))/sqrt(x) inline float clip(float x) { - float a = x*x; + float a = x*x; return ((a + 105)*a + 945) / ((15*a + 420)*a + 945); } + struct LadderFilter { - float g = 0.1f; - float resonance = 0.5f; + float g = 0.f; + float resonance = 0.f; float state[4] = {}; float zi = 0.f; float output[3] = {}; @@ -42,28 +44,28 @@ struct LadderFilter { // update last LP1 output const float t2t3 = t2*t3; - - float y3 = (s[3]*(1+t3) + s[2]*t3)*(1+t2); - y3 = (y3 + t2t3*s[1])*(1+t1); - y3 = (y3 + t1*t2t3*(s[0]+t0*input)); + + float y3 = (state[3]*(1+t3) + state[2]*t3)*(1+t2); + y3 = (y3 + t2t3*state[1])*(1+t1); + y3 = (y3 + t1*t2t3*(state[0]+t0*input)); y3 = y3 / ((1+t1)*(1+t2)*(1+t3)*(1+t4) + resonance*t0*t1*t2t3); - // update other LP1 outputs + // update other LP1 outputs const float xx = t0 * (input - resonance * y3); - const float y0 = t1 * (s[0] + xx) / (1+t1); - const float y1 = t2 * (s[1] + y0) / (1+t2); - const float y2 = t3 * (s[2] + y1) / (1+t3); + const float y0 = t1 * (state[0] + xx) / (1+t1); + const float y1 = t2 * (state[1] + y0) / (1+t2); + const float y2 = t3 * (state[2] + y1) / (1+t3); // update states - s[0] += 2 * (xx - y0); - s[1] += 2 * (y0 - y1); - s[2] += 2 * (y1 - y2); - s[3] += 2 * (y2 - t4*y3); + state[0] += 2 * (xx - y0); + state[1] += 2 * (y0 - y1); + state[2] += 2 * (y1 - y2); + state[3] += 2 * (y2 - t4*y3); // returns LP, HP and BP outputs const float y1t2 = y1 / t2; const float y2t3 = y2 / t3; - + output[0] = y3; output[1] = xx/t0 - 4*y0/t1 + 6*y1t2 - 4*y2t3 + y3; output[2] = y1t2 - 2*y2t3 + y3; @@ -74,7 +76,7 @@ struct LadderFilter { void reset() { for (int i = 0; i < 4; i++) { - state[i] = 0.0f; + state[i] = 0.f; } zi = 0.f; } @@ -114,33 +116,33 @@ struct VCF : Module { void VCF::step() { - float input = inputs[IN_INPUT].value / 5.0f; - float drive = params[DRIVE_PARAM].value + inputs[DRIVE_INPUT].value / 10.0f; - float gain = powf(100.0f, drive); + float input = inputs[IN_INPUT].value / 5.f; + float gain = powf(1.f + params[DRIVE_PARAM].value, 5); + if (inputs[DRIVE_INPUT].active) + gain *= inputs[DRIVE_INPUT].value / 10.f; input *= gain; // Add -60dB noise to bootstrap self-oscillation - input += 1e-6f * (2.0f*randomUniform() - 1.0f); + input += 1e-6f * (2.f*randomUniform() - 1.f); // Set resonance - float res = params[RES_PARAM].value + inputs[RES_INPUT].value / 5.0f; - res = clamp(res, 0.0f, 1.0f); // resonance must be between 0 and 1 - filter.resonance = res; + float res = clamp(params[RES_PARAM].value + inputs[RES_INPUT].value / 10.f, 0.f, 1.f); + filter.resonance = powf(res, 2) * 10.f; // Set cutoff frequency - float cutoffExp = params[FREQ_PARAM].value + params[FREQ_CV_PARAM].value * inputs[FREQ_INPUT].value / 5.0f; - cutoffExp = clamp(cutoffExp, 0.0f, 1.0f); - const float minCutoff = 15.0f; - const float maxCutoff = 20000.0f; - float cutoff = minCutoff * powf(maxCutoff / minCutoff, cutoffExp); - filter.g = tanf(float_Pi * cutoff / engineGetSampleRate()); + float pitch = inputs[FREQ_INPUT].value * quadraticBipolar(params[FREQ_CV_PARAM].value); + pitch += params[FREQ_PARAM].value * 10.f - 3.f; + pitch += quadraticBipolar(params[FINE_PARAM].value * 2.f - 1.f) * 7.f/12.f; + float cutoff = 261.626f * powf(2.f, pitch); + cutoff = clamp(cutoff, 1.f, 20000.f); + filter.g = tanf(M_PI * cutoff * engineGetSampleTime()); // Push a sample to the state filter filter.process(input); // Set outputs - outputs[LPF_OUTPUT].value = 5.0f * filter.output[0]; - outputs[HPF_OUTPUT].value = 5.0f * filter.output[1]; - //outputs[BPF_OUTPUT].value = 5.0f * filter.output[2]; + outputs[LPF_OUTPUT].value = 5.f * filter.output[0]; + outputs[HPF_OUTPUT].value = 5.f * filter.output[1]; + //outputs[BPF_OUTPUT].value = 5.f * filter.output[2]; } @@ -156,11 +158,11 @@ VCFWidget::VCFWidget(VCF *module) : ModuleWidget(module) { addChild(Widget::create(Vec(15, 365))); addChild(Widget::create(Vec(box.size.x-30, 365))); - addParam(ParamWidget::create(Vec(33, 61), module, VCF::FREQ_PARAM, 0.0f, 1.0f, 0.5f)); - addParam(ParamWidget::create(Vec(12, 143), module, VCF::FINE_PARAM, 0.0f, 1.0f, 0.5f)); - addParam(ParamWidget::create(Vec(71, 143), module, VCF::RES_PARAM, 0.0f, 1.0f, 0.0f)); - addParam(ParamWidget::create(Vec(12, 208), module, VCF::FREQ_CV_PARAM, -1.0f, 1.0f, 0.0f)); - addParam(ParamWidget::create(Vec(71, 208), module, VCF::DRIVE_PARAM, 0.0f, 1.0f, 0.0f)); + addParam(ParamWidget::create(Vec(33, 61), module, VCF::FREQ_PARAM, 0.f, 1.f, 0.5f)); + addParam(ParamWidget::create(Vec(12, 143), module, VCF::FINE_PARAM, 0.f, 1.f, 0.5f)); + addParam(ParamWidget::create(Vec(71, 143), module, VCF::RES_PARAM, 0.f, 1.f, 0.f)); + addParam(ParamWidget::create(Vec(12, 208), module, VCF::FREQ_CV_PARAM, -1.f, 1.f, 0.f)); + addParam(ParamWidget::create(Vec(71, 208), module, VCF::DRIVE_PARAM, 0.f, 1.f, 0.f)); addInput(Port::create(Vec(10, 276), Port::INPUT, module, VCF::FREQ_INPUT)); addInput(Port::create(Vec(48, 276), Port::INPUT, module, VCF::RES_INPUT)); From 37160973e6ef2d8ba1cfcdc8d2c9684957acd0dd Mon Sep 17 00:00:00 2001 From: Ivan COHEN Date: Mon, 9 Jul 2018 10:17:36 +0200 Subject: [PATCH 6/8] Updated Fundamental VCF Added the correct extra HP output to original RK4 algorithm, and provided an oversampled 4 times version + a 0df version --- src/VCF.cpp | 504 +++++++++++++++++++++++++++++++++++++++------------- 1 file changed, 376 insertions(+), 128 deletions(-) diff --git a/src/VCF.cpp b/src/VCF.cpp index cef12a1..5f46068 100644 --- a/src/VCF.cpp +++ b/src/VCF.cpp @@ -1,148 +1,396 @@ -// LICENSE TERMS: Copyright 2012 Teemu Voipio -// -// You can use this however you like for pretty much any purpose, -// as long as you don't claim you wrote it. There is no warranty. -// -// Distribution of substantial portions of this code in source form -// must include this copyright notice and list of conditions. -// -// Improved by Ivan COHEN from musical entropy -// More about this algorithm here : -// http://www.kvraudio.com/forum/viewtopic.php?t=349859 -// http://www.kvraudio.com/forum/viewtopic.php?t=207647 - - #include "Fundamental.hpp" #include "dsp/functions.hpp" +#include "dsp/resampler.hpp" + + +// =================================================================================== +/** + The filter DSP code has been derived from + Miller Puckette's code hosted at + https://github.com/ddiakopoulos/MoogLadders/blob/master/src/RKSimulationModel.h + which is licensed for use under the following terms (MIT license): + + Copyright (c) 2015, Miller Puckette. All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ +struct LadderFilterRK4 +{ +public: + // =================================================================================== + void process(float input, float dt, float &outLP, float &outHP) + { + float deriv1[4], deriv2[4], deriv3[4], deriv4[4], tempState[4]; + + calculateDerivatives(input, deriv1, state); + + for (int i = 0; i < 4; i++) + tempState[i] = state[i] + 0.5f * dt * deriv1[i]; + calculateDerivatives(input, deriv2, tempState); + + for (int i = 0; i < 4; i++) + tempState[i] = state[i] + 0.5f * dt * deriv2[i]; + calculateDerivatives(input, deriv3, tempState); + + for (int i = 0; i < 4; i++) + tempState[i] = state[i] + dt * deriv3[i]; + calculateDerivatives(input, deriv4, tempState); + + for (int i = 0; i < 4; i++) + state[i] += (1.0f / 6.0f) * dt * (deriv1[i] + 2.0f * deriv2[i] + 2.0f * deriv3[i] + deriv4[i]); + + outLP = state[3]; + outHP = clip((input - resonance*state[3]) - 4 * state[0] + 6*state[1] - 4*state[2] + state[3]); + } + + void reset() + { + for (int i = 0; i < 4; i++) + state[i] = 0.0f; + } + // =================================================================================== + // parameters + float cutoff = 1000.0f; + float resonance = 1.0f; + +private: + // =================================================================================== + // The clipping function of a transistor pair is approximately tanh(x) + inline float clip(float x) + { + return std::tanh(x); + } -// The clipping function of a transistor pair is approximately tanh(x) -// This one is a Pade-approx for tanh(sqrt(x))/sqrt(x) -inline float clip(float x) { - float a = x*x; - return ((a + 105)*a + 945) / ((15*a + 420)*a + 945); -} + void calculateDerivatives(float input, float *dstate, const float *_state) + { + float cutoff2Pi = 2 * M_PI * cutoff; + float satstatein = clip(input - resonance * _state[3]); + float satstate0 = clip(_state[0]); + float satstate1 = clip(_state[1]); + float satstate2 = clip(_state[2]); + float satstate3 = clip(_state[3]); -struct LadderFilter { - float g = 0.f; - float resonance = 0.f; - float state[4] = {}; - float zi = 0.f; - float output[3] = {}; - - void process(float input) { - // input with half delay, for non-linearities - const float ih = 0.5f * (input + zi); - - // evaluate the non-linear gains - const float t0 = g * clip(ih - resonance * state[3]); - const float t1 = g * clip(state[0]); - const float t2 = g * clip(state[1]); - const float t3 = g * clip(state[2]); - const float t4 = g * clip(state[3]); - - // update last LP1 output - const float t2t3 = t2*t3; - - float y3 = (state[3]*(1+t3) + state[2]*t3)*(1+t2); - y3 = (y3 + t2t3*state[1])*(1+t1); - y3 = (y3 + t1*t2t3*(state[0]+t0*input)); - y3 = y3 / ((1+t1)*(1+t2)*(1+t3)*(1+t4) + resonance*t0*t1*t2t3); - - // update other LP1 outputs - const float xx = t0 * (input - resonance * y3); - const float y0 = t1 * (state[0] + xx) / (1+t1); - const float y1 = t2 * (state[1] + y0) / (1+t2); - const float y2 = t3 * (state[2] + y1) / (1+t3); - - // update states - state[0] += 2 * (xx - y0); - state[1] += 2 * (y0 - y1); - state[2] += 2 * (y1 - y2); - state[3] += 2 * (y2 - t4*y3); - - // returns LP, HP and BP outputs - const float y1t2 = y1 / t2; - const float y2t3 = y2 / t3; - - output[0] = y3; - output[1] = xx/t0 - 4*y0/t1 + 6*y1t2 - 4*y2t3 + y3; - output[2] = y1t2 - 2*y2t3 + y3; - - // update delay input state - zi = input; + dstate[0] = cutoff2Pi * (satstatein - satstate0); + dstate[1] = cutoff2Pi * (satstate0 - satstate1); + dstate[2] = cutoff2Pi * (satstate1 - satstate2); + dstate[3] = cutoff2Pi * (satstate2 - satstate3); } - void reset() { - for (int i = 0; i < 4; i++) { - state[i] = 0.f; + // =================================================================================== + // state variables + float state[4] = {}; +}; + + +// =================================================================================== +/** + Same than before but with oversampling x4 +*/ +struct LadderFilterRK4Oversampled +{ +public: + // =================================================================================== + void process(float input, float dt, float &outLP, float &outHP) + { + // upsampling + float inup[4], outLPup[4], outHPup[4]; + upsampler.process(input, inup); + + // processing + float deriv1[4], deriv2[4], deriv3[4], deriv4[4], tempState[4]; + + for (int n = 0; n < 4; n++) + { + for (int i = 0; i < 4; i++) + tempState[i] = state[i]; + calculateDerivatives(inup[n], deriv1, tempState); + + for (int i = 0; i < 4; i++) + tempState[i] = state[i] + 0.5f * dt * 0.25f * deriv1[i]; + calculateDerivatives(inup[n], deriv2, tempState); + + for (int i = 0; i < 4; i++) + tempState[i] = state[i] + 0.5f * dt * 0.25f * deriv2[i]; + calculateDerivatives(inup[n], deriv3, tempState); + + for (int i = 0; i < 4; i++) + tempState[i] = state[i] + dt * 0.25f * deriv3[i]; + calculateDerivatives(inup[n], deriv4, tempState); + + for (int i = 0; i < 4; i++) + state[i] += (1.0f / 6.0f) * dt * 0.25f * (deriv1[i] + 2.0f * deriv2[i] + 2.0f * deriv3[i] + deriv4[i]); + + outLPup[n] = state[3]; + outHPup[n] = clip((inup[n] - resonance*state[3]) - 4 * state[0] + 6*state[1] - 4*state[2] + state[3]); } - zi = 0.f; + + // downsampling + outLP = decimatorLP.process(outLPup); + outHP = decimatorHP.process(outHPup); + } + + void reset() + { + for (int i = 0; i < 4; i++) + state[i] = 0.0f; + + upsampler.reset(); + decimatorLP.reset(); + decimatorHP.reset(); } + + // =================================================================================== + // parameters + float cutoff = 1000.0f; + float resonance = 1.0f; + +private: + // =================================================================================== + // The clipping function of a transistor pair is approximately tanh(x) + inline float clip(float x) + { + return std::tanh(x); + } + + void calculateDerivatives(float input, float *dstate, const float *_state) + { + float cutoff2Pi = 2 * M_PI * cutoff; + + float satstatein = clip(input - resonance * _state[3]); + float satstate0 = clip(_state[0]); + float satstate1 = clip(_state[1]); + float satstate2 = clip(_state[2]); + float satstate3 = clip(_state[3]); + + dstate[0] = cutoff2Pi * (satstatein - satstate0); + dstate[1] = cutoff2Pi * (satstate0 - satstate1); + dstate[2] = cutoff2Pi * (satstate1 - satstate2); + dstate[3] = cutoff2Pi * (satstate2 - satstate3); + } + + // =================================================================================== + // oversampling classes + Upsampler<4, 32> upsampler; + Decimator<4, 24> decimatorLP, decimatorHP; + + // state variables + float state[4] = {}; }; +// =================================================================================== +/** + This code is derived from the algorithm described in "Modeling and Measuring a + Moog Voltage-Controlled Filter, by Effrosyni Paschou, Fabián Esqueda, Vesa Välimäki + and John Mourjopoulos, for APSIPA Annual Summit and Conference 2017. + + It is a 0df algorithm using Newton-Raphson's method (4 iterations) for solving + implicit equations, and midpoint integration method. + + Derived from the article and adapted to VCV Rack by Ivan COHEN. +*/ +struct LadderFilter0dfMidPoint +{ +public: + // =================================================================================== + // parameters + float cutoff = 1000.f; + float resonance = 1.0f; + + LadderFilter0dfMidPoint() + { + A = gainfactor / (2.f * VT); + } + + void process(float input, float dt, float &outLP, float &outHP) + { + float F[4], V[4]; + const float wc = 2*M_PI*cutoff / A; + + input *= gainfactor; + + for(int i = 0; i < 4; i++) + V[i] = lastV[i]; + + for(int i = 0; i < 4; i++) + { + const float S0 = A * (input + lastinput + resonance * (V[3] + lastV[3])) * 0.5f; + const float S1 = A * (V[0] + lastV[0]) * 0.5f; + const float S2 = A * (V[1] + lastV[1]) * 0.5f; + const float S3 = A * (V[2] + lastV[2]) * 0.5f; + const float S4 = A * (V[3] + lastV[3]) * 0.5f; + + const float t0 = clip(S0); + const float t1 = clip(S1); + const float t2 = clip(S2); + const float t3 = clip(S3); + const float t4 = clip(S4); + + // F function (the one for which we need to find the roots with NR) + F[0] = V[0] - lastV[0] + wc*dt*(t0 + t1); + F[1] = V[1] - lastV[1] - wc*dt*(t1 - t2); + F[2] = V[2] - lastV[2] - wc*dt*(t2 - t3); + F[3] = V[3] - lastV[3] - wc*dt*(t3 - t4); + + const float g = wc*dt*A*0.5f; + + // derivatives of ti + const float J0 = 1 - t0 * t0; + const float J1 = 1 - t1 * t1; + const float J2 = 1 - t2 * t2; + const float J3 = 1 - t3 * t3; + const float J4 = 1 - t4 * t4; + + // Jacobian matrix elements + const float a11 = 1 + g*J1; + const float a14 = resonance*g*J0; + const float a21 = -g*J1; + const float a22 = 1 + g*J2; + const float a32 = -g*J2; + const float a33 = 1 + g*J3; + const float a43 = -g*J3; + const float a44 = 1 + g*J4; + + // Newton-Raphson algorithm with Jacobian inverting + const float deninv = 1.f / (a11*a22*a33*a44 - a14*a21*a32*a43); + + float delta0 = ( F[0]*a22*a33*a44 - F[1]*a14*a32*a43 + F[2]*a14*a22*a43 - F[3]*a14*a22*a33) * deninv; + float delta1 = (-F[0]*a21*a33*a44 + F[1]*a11*a33*a44 - F[2]*a14*a21*a43 + F[3]*a14*a21*a33) * deninv; + float delta2 = ( F[0]*a21*a32*a44 - F[1]*a11*a32*a44 + F[2]*a11*a22*a44 - F[3]*a14*a21*a32) * deninv; + float delta3 = (-F[0]*a21*a32*a43 + F[1]*a11*a32*a43 - F[2]*a11*a22*a43 + F[3]*a11*a22*a33) * deninv; + + delta0 = isfinite(delta0) ? delta0 : 0.f; + delta1 = isfinite(delta1) ? delta1 : 0.f; + delta2 = isfinite(delta2) ? delta2 : 0.f; + delta3 = isfinite(delta3) ? delta3 : 0.f; + + V[0] -= delta0; + V[1] -= delta1; + V[2] -= delta2; + V[3] -= delta3; + } + + lastinput = input; + for(int i = 0; i < 4; i++) + lastV[i] = V[i]; + + // outputs are inverted + outLP = -clip(V[3] / gainfactor); + outHP = -clip(((input + resonance*V[3]) + 4 * V[0] - 6*V[1] + 4*V[2] - V[3]) / gainfactor); + } + + void reset() + { + for (int i = 0; i < 4; i++) + lastV[i] = 0.0f; + + lastinput = 0.f; + } + +private: + // =================================================================================== + inline float clip(float x) + { + return std::tanh(x); + } + + // =================================================================================== + // constants + float VT = 26e-3f; + float gainfactor = 0.2f; + float A; + // state variables + float lastinput = 0.f; + float lastV[4] = {}; +}; + + +// =================================================================================== struct VCF : Module { - enum ParamIds { - FREQ_PARAM, - FINE_PARAM, - RES_PARAM, - FREQ_CV_PARAM, - DRIVE_PARAM, - NUM_PARAMS - }; - enum InputIds { - FREQ_INPUT, - RES_INPUT, - DRIVE_INPUT, - IN_INPUT, - NUM_INPUTS - }; - enum OutputIds { - LPF_OUTPUT, - HPF_OUTPUT, - NUM_OUTPUTS - }; - - LadderFilter filter; - - VCF() : Module(NUM_PARAMS, NUM_INPUTS, NUM_OUTPUTS) {} - void step() override; - void onReset() override { - filter.reset(); - } + enum ParamIds { + FREQ_PARAM, + FINE_PARAM, + RES_PARAM, + FREQ_CV_PARAM, + DRIVE_PARAM, + NUM_PARAMS + }; + enum InputIds { + FREQ_INPUT, + RES_INPUT, + DRIVE_INPUT, + IN_INPUT, + NUM_INPUTS + }; + enum OutputIds { + LPF_OUTPUT, + HPF_OUTPUT, + NUM_OUTPUTS + }; + + //LadderFilter0dfMidPoint filter; + LadderFilterRK4Oversampled filter; + //LadderFilterRK4 filter; + + VCF() : Module(NUM_PARAMS, NUM_INPUTS, NUM_OUTPUTS) {} + void step() override; + void onReset() override { + filter.reset(); + } }; -void VCF::step() { - float input = inputs[IN_INPUT].value / 5.f; - float gain = powf(1.f + params[DRIVE_PARAM].value, 5); - if (inputs[DRIVE_INPUT].active) - gain *= inputs[DRIVE_INPUT].value / 10.f; - input *= gain; - // Add -60dB noise to bootstrap self-oscillation - input += 1e-6f * (2.f*randomUniform() - 1.f); - - // Set resonance - float res = clamp(params[RES_PARAM].value + inputs[RES_INPUT].value / 10.f, 0.f, 1.f); - filter.resonance = powf(res, 2) * 10.f; - - // Set cutoff frequency - float pitch = inputs[FREQ_INPUT].value * quadraticBipolar(params[FREQ_CV_PARAM].value); - pitch += params[FREQ_PARAM].value * 10.f - 3.f; - pitch += quadraticBipolar(params[FINE_PARAM].value * 2.f - 1.f) * 7.f/12.f; - float cutoff = 261.626f * powf(2.f, pitch); - cutoff = clamp(cutoff, 1.f, 20000.f); - filter.g = tanf(M_PI * cutoff * engineGetSampleTime()); - - // Push a sample to the state filter - filter.process(input); - - // Set outputs - outputs[LPF_OUTPUT].value = 5.f * filter.output[0]; - outputs[HPF_OUTPUT].value = 5.f * filter.output[1]; - //outputs[BPF_OUTPUT].value = 5.f * filter.output[2]; +void VCF::step() +{ + // input stage + float input = inputs[IN_INPUT].value / 5.f; + float gain = std::pow(1.f + params[DRIVE_PARAM].value, 5); + + if (inputs[DRIVE_INPUT].active) + gain *= inputs[DRIVE_INPUT].value / 10.f; + + input *= gain; + // Add -60dB noise to bootstrap self-oscillation + input += 1e-6f * (2.f*randomUniform() - 1.f); + + // Set resonance + float res = clamp(params[RES_PARAM].value + inputs[RES_INPUT].value / 10.f, 0.f, 1.f); + filter.resonance = std::pow(res, 2) * 10.f; + + // Set cutoff frequency + float pitch = inputs[FREQ_INPUT].value * quadraticBipolar(params[FREQ_CV_PARAM].value); + pitch += params[FREQ_PARAM].value * 10.f - 3.f; + pitch += quadraticBipolar(params[FINE_PARAM].value * 2.f - 1.f) * 7.f/12.f; + float cutoff = 261.626f * std::pow(2.f, pitch); + filter.cutoff = clamp(cutoff, 1.f, 20000.f); + + // Push a sample to the state filter + float outLP, outHP; + filter.process(input, engineGetSampleTime(), outLP, outHP); + + // Set outputs + outputs[LPF_OUTPUT].value = 5.0f * outLP; + outputs[HPF_OUTPUT].value = 5.0f * outHP; } From 4d5e9de05cf4f00f2c67bbe9489046e6a88fac28 Mon Sep 17 00:00:00 2001 From: Andrew Belt Date: Mon, 9 Jul 2018 09:10:59 -0400 Subject: [PATCH 7/8] Fix code style, keep 0df implementation, remove others --- src/VCF.cpp | 556 ++++++++++++++++------------------------------------ 1 file changed, 174 insertions(+), 382 deletions(-) diff --git a/src/VCF.cpp b/src/VCF.cpp index 5f46068..3c12086 100644 --- a/src/VCF.cpp +++ b/src/VCF.cpp @@ -3,423 +3,215 @@ #include "dsp/resampler.hpp" -// =================================================================================== /** - The filter DSP code has been derived from - Miller Puckette's code hosted at - https://github.com/ddiakopoulos/MoogLadders/blob/master/src/RKSimulationModel.h - which is licensed for use under the following terms (MIT license): - - Copyright (c) 2015, Miller Puckette. All rights reserved. - - Redistribution and use in source and binary forms, with or without - modification, are permitted provided that the following conditions are met: - - * Redistributions of source code must retain the above copyright notice, this - list of conditions and the following disclaimer. - * Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimer in the documentation - and/or other materials provided with the distribution. - - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" - AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE - DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE - FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR - SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER - CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, - OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE - OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -*/ -struct LadderFilterRK4 -{ -public: - // =================================================================================== - void process(float input, float dt, float &outLP, float &outHP) - { - float deriv1[4], deriv2[4], deriv3[4], deriv4[4], tempState[4]; - - calculateDerivatives(input, deriv1, state); - - for (int i = 0; i < 4; i++) - tempState[i] = state[i] + 0.5f * dt * deriv1[i]; - calculateDerivatives(input, deriv2, tempState); - - for (int i = 0; i < 4; i++) - tempState[i] = state[i] + 0.5f * dt * deriv2[i]; - calculateDerivatives(input, deriv3, tempState); - - for (int i = 0; i < 4; i++) - tempState[i] = state[i] + dt * deriv3[i]; - calculateDerivatives(input, deriv4, tempState); - - for (int i = 0; i < 4; i++) - state[i] += (1.0f / 6.0f) * dt * (deriv1[i] + 2.0f * deriv2[i] + 2.0f * deriv3[i] + deriv4[i]); +This code is derived from the algorithm described in "Modeling and Measuring a +Moog Voltage-Controlled Filter, by Effrosyni Paschou, Fabián Esqueda, Vesa Välimäki +and John Mourjopoulos, for APSIPA Annual Summit and Conference 2017. - outLP = state[3]; - outHP = clip((input - resonance*state[3]) - 4 * state[0] + 6*state[1] - 4*state[2] + state[3]); - } - - void reset() - { - for (int i = 0; i < 4; i++) - state[i] = 0.0f; - } +It is a 0df algorithm using Newton-Raphson's method (4 iterations) for solving +implicit equations, and midpoint integration method. - // =================================================================================== - // parameters - float cutoff = 1000.0f; +Derived from the article and adapted to VCV Rack by Ivan COHEN. +*/ +struct LadderFilter0dfMidPoint { + float cutoff = 1000.f; float resonance = 1.0f; + float gainFactor = 0.2f; + float A; + /** State variables */ + float lastInput = 0.f; + float lastV[4] = {}; + /** Outputs */ + float lowpass; + float highpass; -private: - // =================================================================================== - // The clipping function of a transistor pair is approximately tanh(x) - inline float clip(float x) - { - return std::tanh(x); - } - - void calculateDerivatives(float input, float *dstate, const float *_state) - { - float cutoff2Pi = 2 * M_PI * cutoff; - float satstatein = clip(input - resonance * _state[3]); - float satstate0 = clip(_state[0]); - float satstate1 = clip(_state[1]); - float satstate2 = clip(_state[2]); - float satstate3 = clip(_state[3]); + LadderFilter0dfMidPoint() { + float VT = 26e-3f; + A = gainFactor / (2.f * VT); + } - dstate[0] = cutoff2Pi * (satstatein - satstate0); - dstate[1] = cutoff2Pi * (satstate0 - satstate1); - dstate[2] = cutoff2Pi * (satstate1 - satstate2); - dstate[3] = cutoff2Pi * (satstate2 - satstate3); + inline float clip(float x) { + return tanhf(x); } - // =================================================================================== - // state variables - float state[4] = {}; -}; + void process(float input, float dt) { + float F[4]; + float V[4]; + float wc = 2 * M_PI * cutoff / A; + input *= gainFactor; -// =================================================================================== -/** - Same than before but with oversampling x4 -*/ -struct LadderFilterRK4Oversampled -{ -public: - // =================================================================================== - void process(float input, float dt, float &outLP, float &outHP) - { - // upsampling - float inup[4], outLPup[4], outHPup[4]; - upsampler.process(input, inup); - - // processing - float deriv1[4], deriv2[4], deriv3[4], deriv4[4], tempState[4]; - - for (int n = 0; n < 4; n++) - { - for (int i = 0; i < 4; i++) - tempState[i] = state[i]; - calculateDerivatives(inup[n], deriv1, tempState); - - for (int i = 0; i < 4; i++) - tempState[i] = state[i] + 0.5f * dt * 0.25f * deriv1[i]; - calculateDerivatives(inup[n], deriv2, tempState); - - for (int i = 0; i < 4; i++) - tempState[i] = state[i] + 0.5f * dt * 0.25f * deriv2[i]; - calculateDerivatives(inup[n], deriv3, tempState); - - for (int i = 0; i < 4; i++) - tempState[i] = state[i] + dt * 0.25f * deriv3[i]; - calculateDerivatives(inup[n], deriv4, tempState); - - for (int i = 0; i < 4; i++) - state[i] += (1.0f / 6.0f) * dt * 0.25f * (deriv1[i] + 2.0f * deriv2[i] + 2.0f * deriv3[i] + deriv4[i]); - - outLPup[n] = state[3]; - outHPup[n] = clip((inup[n] - resonance*state[3]) - 4 * state[0] + 6*state[1] - 4*state[2] + state[3]); + for (int i = 0; i < 4; i++) { + V[i] = lastV[i]; } - // downsampling - outLP = decimatorLP.process(outLPup); - outHP = decimatorHP.process(outHPup); - } + for (int i = 0; i < 4; i++) { + float S0 = A * (input + lastInput + resonance * (V[3] + lastV[3])) * 0.5f; + float S1 = A * (V[0] + lastV[0]) * 0.5f; + float S2 = A * (V[1] + lastV[1]) * 0.5f; + float S3 = A * (V[2] + lastV[2]) * 0.5f; + float S4 = A * (V[3] + lastV[3]) * 0.5f; + + float t0 = clip(S0); + float t1 = clip(S1); + float t2 = clip(S2); + float t3 = clip(S3); + float t4 = clip(S4); + + // F function (the one for which we need to find the roots with NR) + F[0] = V[0] - lastV[0] + wc * dt * (t0 + t1); + F[1] = V[1] - lastV[1] - wc * dt * (t1 - t2); + F[2] = V[2] - lastV[2] - wc * dt * (t2 - t3); + F[3] = V[3] - lastV[3] - wc * dt * (t3 - t4); + + float g = wc * dt * A * 0.5f; + + // derivatives of ti + float J0 = 1 - t0 * t0; + float J1 = 1 - t1 * t1; + float J2 = 1 - t2 * t2; + float J3 = 1 - t3 * t3; + float J4 = 1 - t4 * t4; + + // Jacobian matrix elements + float a11 = 1 + g * J1; + float a14 = resonance * g * J0; + float a21 = -g * J1; + float a22 = 1 + g * J2; + float a32 = -g * J2; + float a33 = 1 + g * J3; + float a43 = -g * J3; + float a44 = 1 + g * J4; + + // Newton-Raphson algorithm with Jacobian inverting + float deninv = 1.f / (a11 * a22 * a33 * a44 - a14 * a21 * a32 * a43); + + float delta0 = ( F[0] * a22 * a33 * a44 - F[1] * a14 * a32 * a43 + F[2] * a14 * a22 * a43 - F[3] * a14 * a22 * a33) * deninv; + float delta1 = (-F[0] * a21 * a33 * a44 + F[1] * a11 * a33 * a44 - F[2] * a14 * a21 * a43 + F[3] * a14 * a21 * a33) * deninv; + float delta2 = ( F[0] * a21 * a32 * a44 - F[1] * a11 * a32 * a44 + F[2] * a11 * a22 * a44 - F[3] * a14 * a21 * a32) * deninv; + float delta3 = (-F[0] * a21 * a32 * a43 + F[1] * a11 * a32 * a43 - F[2] * a11 * a22 * a43 + F[3] * a11 * a22 * a33) * deninv; + + delta0 = isfinite(delta0) ? delta0 : 0.f; + delta1 = isfinite(delta1) ? delta1 : 0.f; + delta2 = isfinite(delta2) ? delta2 : 0.f; + delta3 = isfinite(delta3) ? delta3 : 0.f; + + V[0] -= delta0; + V[1] -= delta1; + V[2] -= delta2; + V[3] -= delta3; + } - void reset() - { + lastInput = input; for (int i = 0; i < 4; i++) - state[i] = 0.0f; + lastV[i] = V[i]; - upsampler.reset(); - decimatorLP.reset(); - decimatorHP.reset(); + // outputs are inverted + lowpass = -clip(V[3] / gainFactor); + highpass = -clip(((input + resonance * V[3]) + 4 * V[0] - 6 * V[1] + 4 * V[2] - V[3]) / gainFactor); } - // =================================================================================== - // parameters - float cutoff = 1000.0f; - float resonance = 1.0f; - -private: - // =================================================================================== - // The clipping function of a transistor pair is approximately tanh(x) - inline float clip(float x) - { - return std::tanh(x); + void reset() { + for (int i = 0; i < 4; i++) { + lastV[i] = 0.f; + } + lastInput = 0.f; } +}; - void calculateDerivatives(float input, float *dstate, const float *_state) - { - float cutoff2Pi = 2 * M_PI * cutoff; - - float satstatein = clip(input - resonance * _state[3]); - float satstate0 = clip(_state[0]); - float satstate1 = clip(_state[1]); - float satstate2 = clip(_state[2]); - float satstate3 = clip(_state[3]); - dstate[0] = cutoff2Pi * (satstatein - satstate0); - dstate[1] = cutoff2Pi * (satstate0 - satstate1); - dstate[2] = cutoff2Pi * (satstate1 - satstate2); - dstate[3] = cutoff2Pi * (satstate2 - satstate3); +struct VCF : Module { + enum ParamIds { + FREQ_PARAM, + FINE_PARAM, + RES_PARAM, + FREQ_CV_PARAM, + DRIVE_PARAM, + NUM_PARAMS + }; + enum InputIds { + FREQ_INPUT, + RES_INPUT, + DRIVE_INPUT, + IN_INPUT, + NUM_INPUTS + }; + enum OutputIds { + LPF_OUTPUT, + HPF_OUTPUT, + NUM_OUTPUTS + }; + + LadderFilter0dfMidPoint filter; + + VCF() : Module(NUM_PARAMS, NUM_INPUTS, NUM_OUTPUTS) {} + + void onReset() override { + filter.reset(); } - // =================================================================================== - // oversampling classes - Upsampler<4, 32> upsampler; - Decimator<4, 24> decimatorLP, decimatorHP; - - // state variables - float state[4] = {}; -}; - -// =================================================================================== -/** - This code is derived from the algorithm described in "Modeling and Measuring a - Moog Voltage-Controlled Filter, by Effrosyni Paschou, Fabián Esqueda, Vesa Välimäki - and John Mourjopoulos, for APSIPA Annual Summit and Conference 2017. + void step() override { + if (!outputs[LPF_OUTPUT].active && !outputs[HPF_OUTPUT].active) { + outputs[LPF_OUTPUT].value = 0.f; + outputs[HPF_OUTPUT].value = 0.f; + return; + } - It is a 0df algorithm using Newton-Raphson's method (4 iterations) for solving - implicit equations, and midpoint integration method. + float input = inputs[IN_INPUT].value / 5.f; + float drive = clamp(params[DRIVE_PARAM].value + inputs[DRIVE_INPUT].value / 10.f, 0.f, 1.f); + float gain = powf(1.f + drive, 5); + input *= gain; - Derived from the article and adapted to VCV Rack by Ivan COHEN. -*/ -struct LadderFilter0dfMidPoint -{ -public: - // =================================================================================== - // parameters - float cutoff = 1000.f; - float resonance = 1.0f; - - LadderFilter0dfMidPoint() - { - A = gainfactor / (2.f * VT); - } - - void process(float input, float dt, float &outLP, float &outHP) - { - float F[4], V[4]; - const float wc = 2*M_PI*cutoff / A; - - input *= gainfactor; - - for(int i = 0; i < 4; i++) - V[i] = lastV[i]; - - for(int i = 0; i < 4; i++) - { - const float S0 = A * (input + lastinput + resonance * (V[3] + lastV[3])) * 0.5f; - const float S1 = A * (V[0] + lastV[0]) * 0.5f; - const float S2 = A * (V[1] + lastV[1]) * 0.5f; - const float S3 = A * (V[2] + lastV[2]) * 0.5f; - const float S4 = A * (V[3] + lastV[3]) * 0.5f; - - const float t0 = clip(S0); - const float t1 = clip(S1); - const float t2 = clip(S2); - const float t3 = clip(S3); - const float t4 = clip(S4); - - // F function (the one for which we need to find the roots with NR) - F[0] = V[0] - lastV[0] + wc*dt*(t0 + t1); - F[1] = V[1] - lastV[1] - wc*dt*(t1 - t2); - F[2] = V[2] - lastV[2] - wc*dt*(t2 - t3); - F[3] = V[3] - lastV[3] - wc*dt*(t3 - t4); - - const float g = wc*dt*A*0.5f; - - // derivatives of ti - const float J0 = 1 - t0 * t0; - const float J1 = 1 - t1 * t1; - const float J2 = 1 - t2 * t2; - const float J3 = 1 - t3 * t3; - const float J4 = 1 - t4 * t4; - - // Jacobian matrix elements - const float a11 = 1 + g*J1; - const float a14 = resonance*g*J0; - const float a21 = -g*J1; - const float a22 = 1 + g*J2; - const float a32 = -g*J2; - const float a33 = 1 + g*J3; - const float a43 = -g*J3; - const float a44 = 1 + g*J4; - - // Newton-Raphson algorithm with Jacobian inverting - const float deninv = 1.f / (a11*a22*a33*a44 - a14*a21*a32*a43); - - float delta0 = ( F[0]*a22*a33*a44 - F[1]*a14*a32*a43 + F[2]*a14*a22*a43 - F[3]*a14*a22*a33) * deninv; - float delta1 = (-F[0]*a21*a33*a44 + F[1]*a11*a33*a44 - F[2]*a14*a21*a43 + F[3]*a14*a21*a33) * deninv; - float delta2 = ( F[0]*a21*a32*a44 - F[1]*a11*a32*a44 + F[2]*a11*a22*a44 - F[3]*a14*a21*a32) * deninv; - float delta3 = (-F[0]*a21*a32*a43 + F[1]*a11*a32*a43 - F[2]*a11*a22*a43 + F[3]*a11*a22*a33) * deninv; - - delta0 = isfinite(delta0) ? delta0 : 0.f; - delta1 = isfinite(delta1) ? delta1 : 0.f; - delta2 = isfinite(delta2) ? delta2 : 0.f; - delta3 = isfinite(delta3) ? delta3 : 0.f; - - V[0] -= delta0; - V[1] -= delta1; - V[2] -= delta2; - V[3] -= delta3; - } - - lastinput = input; - for(int i = 0; i < 4; i++) - lastV[i] = V[i]; - - // outputs are inverted - outLP = -clip(V[3] / gainfactor); - outHP = -clip(((input + resonance*V[3]) + 4 * V[0] - 6*V[1] + 4*V[2] - V[3]) / gainfactor); - } - - void reset() - { - for (int i = 0; i < 4; i++) - lastV[i] = 0.0f; - - lastinput = 0.f; - } - -private: - // =================================================================================== - inline float clip(float x) - { - return std::tanh(x); - } + // Add -60dB noise to bootstrap self-oscillation + input += 1e-6f * (2.f * randomUniform() - 1.f); - // =================================================================================== - // constants - float VT = 26e-3f; - float gainfactor = 0.2f; - float A; + // Set resonance + float res = clamp(params[RES_PARAM].value + inputs[RES_INPUT].value / 10.f, 0.f, 1.f); + filter.resonance = powf(res, 2) * 10.f; - // state variables - float lastinput = 0.f; - float lastV[4] = {}; -}; + // Set cutoff frequency + float pitch = inputs[FREQ_INPUT].value * quadraticBipolar(params[FREQ_CV_PARAM].value); + pitch += params[FREQ_PARAM].value * 10.f - 3.f; + pitch += quadraticBipolar(params[FINE_PARAM].value * 2.f - 1.f) * 7.f / 12.f; + float cutoff = 261.626f * powf(2.f, pitch); + filter.cutoff = clamp(cutoff, 1.f, 20000.f); + // Step the filter + filter.process(input, engineGetSampleTime()); -// =================================================================================== -struct VCF : Module { - enum ParamIds { - FREQ_PARAM, - FINE_PARAM, - RES_PARAM, - FREQ_CV_PARAM, - DRIVE_PARAM, - NUM_PARAMS - }; - enum InputIds { - FREQ_INPUT, - RES_INPUT, - DRIVE_INPUT, - IN_INPUT, - NUM_INPUTS - }; - enum OutputIds { - LPF_OUTPUT, - HPF_OUTPUT, - NUM_OUTPUTS - }; - - //LadderFilter0dfMidPoint filter; - LadderFilterRK4Oversampled filter; - //LadderFilterRK4 filter; - - VCF() : Module(NUM_PARAMS, NUM_INPUTS, NUM_OUTPUTS) {} - void step() override; - void onReset() override { - filter.reset(); - } + // Set outputs + outputs[LPF_OUTPUT].value = 5.f * filter.lowpass; + outputs[HPF_OUTPUT].value = 5.f * filter.highpass; + } }; -void VCF::step() -{ - // input stage - float input = inputs[IN_INPUT].value / 5.f; - float gain = std::pow(1.f + params[DRIVE_PARAM].value, 5); - - if (inputs[DRIVE_INPUT].active) - gain *= inputs[DRIVE_INPUT].value / 10.f; - - input *= gain; - // Add -60dB noise to bootstrap self-oscillation - input += 1e-6f * (2.f*randomUniform() - 1.f); - - // Set resonance - float res = clamp(params[RES_PARAM].value + inputs[RES_INPUT].value / 10.f, 0.f, 1.f); - filter.resonance = std::pow(res, 2) * 10.f; - - // Set cutoff frequency - float pitch = inputs[FREQ_INPUT].value * quadraticBipolar(params[FREQ_CV_PARAM].value); - pitch += params[FREQ_PARAM].value * 10.f - 3.f; - pitch += quadraticBipolar(params[FINE_PARAM].value * 2.f - 1.f) * 7.f/12.f; - float cutoff = 261.626f * std::pow(2.f, pitch); - filter.cutoff = clamp(cutoff, 1.f, 20000.f); - - // Push a sample to the state filter - float outLP, outHP; - filter.process(input, engineGetSampleTime(), outLP, outHP); - - // Set outputs - outputs[LPF_OUTPUT].value = 5.0f * outLP; - outputs[HPF_OUTPUT].value = 5.0f * outHP; -} - - struct VCFWidget : ModuleWidget { - VCFWidget(VCF *module); + VCFWidget(VCF *module) : ModuleWidget(module) { + setPanel(SVG::load(assetPlugin(plugin, "res/VCF.svg"))); + + addChild(Widget::create(Vec(15, 0))); + addChild(Widget::create(Vec(box.size.x - 30, 0))); + addChild(Widget::create(Vec(15, 365))); + addChild(Widget::create(Vec(box.size.x - 30, 365))); + + addParam(ParamWidget::create(Vec(33, 61), module, VCF::FREQ_PARAM, 0.f, 1.f, 0.5f)); + addParam(ParamWidget::create(Vec(12, 143), module, VCF::FINE_PARAM, 0.f, 1.f, 0.5f)); + addParam(ParamWidget::create(Vec(71, 143), module, VCF::RES_PARAM, 0.f, 1.f, 0.f)); + addParam(ParamWidget::create(Vec(12, 208), module, VCF::FREQ_CV_PARAM, -1.f, 1.f, 0.f)); + addParam(ParamWidget::create(Vec(71, 208), module, VCF::DRIVE_PARAM, 0.f, 1.f, 0.f)); + + addInput(Port::create(Vec(10, 276), Port::INPUT, module, VCF::FREQ_INPUT)); + addInput(Port::create(Vec(48, 276), Port::INPUT, module, VCF::RES_INPUT)); + addInput(Port::create(Vec(85, 276), Port::INPUT, module, VCF::DRIVE_INPUT)); + addInput(Port::create(Vec(10, 320), Port::INPUT, module, VCF::IN_INPUT)); + + addOutput(Port::create(Vec(48, 320), Port::OUTPUT, module, VCF::LPF_OUTPUT)); + addOutput(Port::create(Vec(85, 320), Port::OUTPUT, module, VCF::HPF_OUTPUT)); + } }; -VCFWidget::VCFWidget(VCF *module) : ModuleWidget(module) { - setPanel(SVG::load(assetPlugin(plugin, "res/VCF.svg"))); - - addChild(Widget::create(Vec(15, 0))); - addChild(Widget::create(Vec(box.size.x-30, 0))); - addChild(Widget::create(Vec(15, 365))); - addChild(Widget::create(Vec(box.size.x-30, 365))); - - addParam(ParamWidget::create(Vec(33, 61), module, VCF::FREQ_PARAM, 0.f, 1.f, 0.5f)); - addParam(ParamWidget::create(Vec(12, 143), module, VCF::FINE_PARAM, 0.f, 1.f, 0.5f)); - addParam(ParamWidget::create(Vec(71, 143), module, VCF::RES_PARAM, 0.f, 1.f, 0.f)); - addParam(ParamWidget::create(Vec(12, 208), module, VCF::FREQ_CV_PARAM, -1.f, 1.f, 0.f)); - addParam(ParamWidget::create(Vec(71, 208), module, VCF::DRIVE_PARAM, 0.f, 1.f, 0.f)); - - addInput(Port::create(Vec(10, 276), Port::INPUT, module, VCF::FREQ_INPUT)); - addInput(Port::create(Vec(48, 276), Port::INPUT, module, VCF::RES_INPUT)); - addInput(Port::create(Vec(85, 276), Port::INPUT, module, VCF::DRIVE_INPUT)); - addInput(Port::create(Vec(10, 320), Port::INPUT, module, VCF::IN_INPUT)); - - addOutput(Port::create(Vec(48, 320), Port::OUTPUT, module, VCF::LPF_OUTPUT)); - addOutput(Port::create(Vec(85, 320), Port::OUTPUT, module, VCF::HPF_OUTPUT)); -} Model *modelVCF = Model::create("Fundamental", "VCF", "VCF", FILTER_TAG); From 8bad806f7927c72df2c40105ded1a324913536fe Mon Sep 17 00:00:00 2001 From: Andrew Belt Date: Mon, 9 Jul 2018 10:52:48 -0400 Subject: [PATCH 8/8] Rewrite LadderFilter --- src/VCF.cpp | 175 ++++++++++++++++++++-------------------------------- 1 file changed, 66 insertions(+), 109 deletions(-) diff --git a/src/VCF.cpp b/src/VCF.cpp index 3c12086..7b51679 100644 --- a/src/VCF.cpp +++ b/src/VCF.cpp @@ -1,126 +1,58 @@ #include "Fundamental.hpp" #include "dsp/functions.hpp" #include "dsp/resampler.hpp" +#include "dsp/ode.hpp" -/** -This code is derived from the algorithm described in "Modeling and Measuring a -Moog Voltage-Controlled Filter, by Effrosyni Paschou, Fabián Esqueda, Vesa Välimäki -and John Mourjopoulos, for APSIPA Annual Summit and Conference 2017. +inline float clip(float x) { + return tanhf(x); +} -It is a 0df algorithm using Newton-Raphson's method (4 iterations) for solving -implicit equations, and midpoint integration method. - -Derived from the article and adapted to VCV Rack by Ivan COHEN. -*/ -struct LadderFilter0dfMidPoint { - float cutoff = 1000.f; +struct LadderFilter { + float omega0; float resonance = 1.0f; - float gainFactor = 0.2f; - float A; - /** State variables */ - float lastInput = 0.f; - float lastV[4] = {}; - /** Outputs */ + float state[4]; + float input; float lowpass; float highpass; - - LadderFilter0dfMidPoint() { - float VT = 26e-3f; - A = gainFactor / (2.f * VT); - } - - inline float clip(float x) { - return tanhf(x); + LadderFilter() { + reset(); + setCutoff(0.f); } - void process(float input, float dt) { - float F[4]; - float V[4]; - float wc = 2 * M_PI * cutoff / A; - - input *= gainFactor; - - for (int i = 0; i < 4; i++) { - V[i] = lastV[i]; - } - + void reset() { for (int i = 0; i < 4; i++) { - float S0 = A * (input + lastInput + resonance * (V[3] + lastV[3])) * 0.5f; - float S1 = A * (V[0] + lastV[0]) * 0.5f; - float S2 = A * (V[1] + lastV[1]) * 0.5f; - float S3 = A * (V[2] + lastV[2]) * 0.5f; - float S4 = A * (V[3] + lastV[3]) * 0.5f; - - float t0 = clip(S0); - float t1 = clip(S1); - float t2 = clip(S2); - float t3 = clip(S3); - float t4 = clip(S4); - - // F function (the one for which we need to find the roots with NR) - F[0] = V[0] - lastV[0] + wc * dt * (t0 + t1); - F[1] = V[1] - lastV[1] - wc * dt * (t1 - t2); - F[2] = V[2] - lastV[2] - wc * dt * (t2 - t3); - F[3] = V[3] - lastV[3] - wc * dt * (t3 - t4); - - float g = wc * dt * A * 0.5f; - - // derivatives of ti - float J0 = 1 - t0 * t0; - float J1 = 1 - t1 * t1; - float J2 = 1 - t2 * t2; - float J3 = 1 - t3 * t3; - float J4 = 1 - t4 * t4; - - // Jacobian matrix elements - float a11 = 1 + g * J1; - float a14 = resonance * g * J0; - float a21 = -g * J1; - float a22 = 1 + g * J2; - float a32 = -g * J2; - float a33 = 1 + g * J3; - float a43 = -g * J3; - float a44 = 1 + g * J4; - - // Newton-Raphson algorithm with Jacobian inverting - float deninv = 1.f / (a11 * a22 * a33 * a44 - a14 * a21 * a32 * a43); - - float delta0 = ( F[0] * a22 * a33 * a44 - F[1] * a14 * a32 * a43 + F[2] * a14 * a22 * a43 - F[3] * a14 * a22 * a33) * deninv; - float delta1 = (-F[0] * a21 * a33 * a44 + F[1] * a11 * a33 * a44 - F[2] * a14 * a21 * a43 + F[3] * a14 * a21 * a33) * deninv; - float delta2 = ( F[0] * a21 * a32 * a44 - F[1] * a11 * a32 * a44 + F[2] * a11 * a22 * a44 - F[3] * a14 * a21 * a32) * deninv; - float delta3 = (-F[0] * a21 * a32 * a43 + F[1] * a11 * a32 * a43 - F[2] * a11 * a22 * a43 + F[3] * a11 * a22 * a33) * deninv; - - delta0 = isfinite(delta0) ? delta0 : 0.f; - delta1 = isfinite(delta1) ? delta1 : 0.f; - delta2 = isfinite(delta2) ? delta2 : 0.f; - delta3 = isfinite(delta3) ? delta3 : 0.f; - - V[0] -= delta0; - V[1] -= delta1; - V[2] -= delta2; - V[3] -= delta3; + state[i] = 0.f; } + } - lastInput = input; - for (int i = 0; i < 4; i++) - lastV[i] = V[i]; - - // outputs are inverted - lowpass = -clip(V[3] / gainFactor); - highpass = -clip(((input + resonance * V[3]) + 4 * V[0] - 6 * V[1] + 4 * V[2] - V[3]) / gainFactor); + void setCutoff(float cutoff) { + omega0 = 2.f*M_PI * cutoff; } - void reset() { - for (int i = 0; i < 4; i++) { - lastV[i] = 0.f; - } - lastInput = 0.f; + void process(float input, float dt) { + stepRK4(0.f, dt, state, 4, [&](float t, const float y[], float dydt[]) { + float inputc = clip(input - resonance * y[3]); + float yc0 = clip(y[0]); + float yc1 = clip(y[1]); + float yc2 = clip(y[2]); + float yc3 = clip(y[3]); + + dydt[0] = omega0 * (inputc - yc0); + dydt[1] = omega0 * (yc0 - yc1); + dydt[2] = omega0 * (yc1 - yc2); + dydt[3] = omega0 * (yc2 - yc3); + }); + + lowpass = state[3]; + highpass = clip((input - resonance*state[3]) - 4 * state[0] + 6*state[1] - 4*state[2] + state[3]); } }; +static const int UPSAMPLE = 2; + struct VCF : Module { enum ParamIds { FREQ_PARAM, @@ -143,7 +75,10 @@ struct VCF : Module { NUM_OUTPUTS }; - LadderFilter0dfMidPoint filter; + LadderFilter filter; + // Upsampler inputUpsampler; + // Decimator lowpassDecimator; + // Decimator highpassDecimator; VCF() : Module(NUM_PARAMS, NUM_INPUTS, NUM_OUTPUTS) {} @@ -171,16 +106,38 @@ struct VCF : Module { filter.resonance = powf(res, 2) * 10.f; // Set cutoff frequency - float pitch = inputs[FREQ_INPUT].value * quadraticBipolar(params[FREQ_CV_PARAM].value); - pitch += params[FREQ_PARAM].value * 10.f - 3.f; + float pitch = 0.f; + if (inputs[FREQ_INPUT].active) + pitch += inputs[FREQ_INPUT].value * quadraticBipolar(params[FREQ_CV_PARAM].value); + pitch += params[FREQ_PARAM].value * 10.f - 5.f; pitch += quadraticBipolar(params[FINE_PARAM].value * 2.f - 1.f) * 7.f / 12.f; float cutoff = 261.626f * powf(2.f, pitch); - filter.cutoff = clamp(cutoff, 1.f, 20000.f); - - // Step the filter - filter.process(input, engineGetSampleTime()); + cutoff = clamp(cutoff, 1.f, 8000.f); + filter.setCutoff(cutoff); + + /* + // Process sample + float dt = engineGetSampleTime() / UPSAMPLE; + float inputBuf[UPSAMPLE]; + float lowpassBuf[UPSAMPLE]; + float highpassBuf[UPSAMPLE]; + inputUpsampler.process(input, inputBuf); + for (int i = 0; i < UPSAMPLE; i++) { + // Step the filter + filter.process(inputBuf[i], dt); + lowpassBuf[i] = filter.lowpass; + highpassBuf[i] = filter.highpass; + } // Set outputs + if (outputs[LPF_OUTPUT].active) { + outputs[LPF_OUTPUT].value = 5.f * lowpassDecimator.process(lowpassBuf); + } + if (outputs[HPF_OUTPUT].active) { + outputs[HPF_OUTPUT].value = 5.f * highpassDecimator.process(highpassBuf); + } + */ + filter.process(input, engineGetSampleTime()); outputs[LPF_OUTPUT].value = 5.f * filter.lowpass; outputs[HPF_OUTPUT].value = 5.f * filter.highpass; }