AnalogTapeModel

Physical modelling signal processing for analog tape recording
Log | Files | Refs | Submodules | README | LICENSE

commit cf0a1537ad20610d4eebd22df5a5ba88039d7f88
parent 83aafdfce9409d79bc755afd2c163b00a972210f
Author: jatinchowdhury18 <jatinchowdhury18@gmail.com>
Date:   Wed, 22 Jul 2020 22:18:57 -0700

Update hysteresis processing (#56)

* Update hysteresis processing to use different ODE solvers

* Add version 1 mode to hysteresis process

* Use double precision for hysteresis processing

* Fine-tune V1 parameters

* Minor hysteresis tweaks

Co-authored-by: jatinchowdhury18 <jatinchowdhury18@users.noreply.github.com>
Diffstat:
MPlugin/Source/GUI/gui.xml | 1+
MPlugin/Source/Processors/Hysteresis/DCFilters.h | 2+-
MPlugin/Source/Processors/Hysteresis/HysteresisProcessing.cpp | 171+++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------
MPlugin/Source/Processors/Hysteresis/HysteresisProcessing.h | 104+++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------
MPlugin/Source/Processors/Hysteresis/HysteresisProcessor.cpp | 215++++++++++++++++++++++++++++++++++++++++++++++---------------------------------
MPlugin/Source/Processors/Hysteresis/HysteresisProcessor.h | 26++++++++++++++++----------
6 files changed, 334 insertions(+), 185 deletions(-)

diff --git a/Plugin/Source/GUI/gui.xml b/Plugin/Source/GUI/gui.xml @@ -46,6 +46,7 @@ <Slider caption="Saturation" parameter="sat" slider-textbox="textbox-below" slider-type="rotary-horizontal-vertical"/> <Slider caption="Drive" parameter="drive" slider-type="rotary-horizontal-vertical" slider-textbox="textbox-below"/> + <ComboBox caption="Mode" parameter="mode" max-height="70"/> </View> <View display="tabbed" padding="0"> <View flex-direction="column" tab-caption="Loss" tab-color="" background-color="" diff --git a/Plugin/Source/Processors/Hysteresis/DCFilters.h b/Plugin/Source/Processors/Hysteresis/DCFilters.h @@ -38,7 +38,7 @@ public: buffer[n] = processSample (buffer[n]); } - inline float processSample (float x) + inline float processSample (float x) noexcept { // direct form II transposed float y = z[1] + x * b[0]; diff --git a/Plugin/Source/Processors/Hysteresis/HysteresisProcessing.cpp b/Plugin/Source/Processors/Hysteresis/HysteresisProcessing.cpp @@ -1,9 +1,15 @@ #include "HysteresisProcessing.h" #include <math.h> -inline int sign (float x) +namespace { - return (x > 0.0f) - (x < 0.0f); + constexpr double ONE_THIRD = 1.0 / 3.0; + constexpr double NEG_TWO_OVER_15 = -2.0 / 15.0; +} + +inline int sign (double x) +{ + return (x > 0.0) - (x < 0.0); } HysteresisProcessing::HysteresisProcessing() @@ -12,108 +18,163 @@ HysteresisProcessing::HysteresisProcessing() void HysteresisProcessing::reset() { - M_n1 = 0.0f; - H_n1 = 0.0f; - H_d_n1 = 0.0f; + M_n1 = 0.0; + H_n1 = 0.0; + H_d_n1 = 0.0; - coth = 0.0f; + coth = 0.0; nearZero = false; } -void HysteresisProcessing::cook (float drive, float width, float sat) +void HysteresisProcessing::setSampleRate (double newSR) +{ + fs = newSR; + T = 1.0 / fs; + Talpha = T / 1.9; +} + +void HysteresisProcessing::cook (float drive, float width, float sat, bool v1) { - M_s = 0.5f + 1.5f * (1.0f - sat); - a = M_s / (0.01f + 6.0f * drive); - c = sqrtf (1.0f - width) - 0.01f; + M_s = 0.5 + 1.5 * (1.0 - (double) sat); + a = M_s / (0.01 + 6.0 * (double) drive); + c = std::sqrt (1.0f - (double) width) - 0.01; + k = 0.47875; + upperLim = 20.0; + + if (v1) + { + k = 27.0e3; + c = 1.7e-1; + M_s *= 50000.0; + a = M_s / (0.01 + 40.0 * (double) drive); + upperLim = 100000.0; + } - nc = 1.0f - c; + nc = 1.0 - c; M_s_oa = M_s / a; + M_s_oa_talpha = alpha * M_s_oa; M_s_oa_tc = c * M_s_oa; M_s_oa_tc_talpha = alpha * M_s_oa_tc; + M_s_oaSq_tc_talpha = M_s_oa_tc_talpha / a; + M_s_oaSq_tc_talphaSq = alpha * M_s_oaSq_tc_talpha; } -inline float HysteresisProcessing::cothApprox (float x) +void HysteresisProcessing::setSolver (SolverType solverType) { - const auto exp_2x = expf (2.0f * x); // expf_approx (2.0f * x); + // defaults + numIter = 0; + solver = &HysteresisProcessing::NR; + + switch (solverType) + { + case SolverType::RK4: + solver = &HysteresisProcessing::RK4; + return; - return (exp_2x + 1.0f) / (exp_2x - 1.0f); + case SolverType::NR5: + numIter = 5; + return; - // const auto xSquare = x * x; + case SolverType::NR10: + numIter = 10; + return; - // return (1.0f + (xSquare / (3.0f + (xSquare / (5.0f + (xSquare / (7.0f))))))) / x; - // return (1.0f + (xSquare / (3.0f + (xSquare / (5.0f))))) / x; - // return (1.0f + (xSquare / 3.0f)) / x; + default: // RK2 + solver = &HysteresisProcessing::RK2; + }; } -inline float HysteresisProcessing::langevin (float x) +inline double HysteresisProcessing::langevin (double x) const noexcept { if (! nearZero) - return (coth) - (1.0f / x); + return (coth) - (1.0 / x); else - return (x / 3.0f); + return x / 3.0; } -inline float HysteresisProcessing::langevinD (float x) +inline double HysteresisProcessing::langevinD (double x) const noexcept { if (! nearZero) - return (1.0f / (x * x)) - (coth * coth) + 1.0f; + return (1.0 / (x * x)) - (coth * coth) + 1.0; else - return (1.0f / 3.0f); + return ONE_THIRD; } -inline float HysteresisProcessing::deriv (float x_n, float x_n1, float x_d_n1) +inline double HysteresisProcessing::langevinD2 (double x) const noexcept { - constexpr float dAlpha = 0.9f; - return (((1.0f + dAlpha) / T) * (x_n - x_n1)) - dAlpha * x_d_n1; + if (! nearZero) + return 2.0 * coth * (coth * coth - 1.0) - (2.0 / (x * x * x)); + else + return NEG_TWO_OVER_15 * x;; } -inline float HysteresisProcessing::hysteresisFunc (float M, float H, float H_d) +inline double HysteresisProcessing::hysteresisFunc (double M, double H, double H_d) noexcept { Q = (H + alpha * M) / a; - coth = 1.0f / std::tanh (Q); - nearZero = Q < 0.001f && Q > -0.001f; + coth = 1.0 / std::tanh (Q); + nearZero = Q < 0.001 && Q > -0.001; M_diff = M_s * langevin (Q) - M; - delta = (float) ((H_d >= 0.0f) - (H_d < 0.0f)); - delta_M = (float) (sign (delta) == sign (M_diff)); + delta = (double) ((H_d >= 0.0) - (H_d < 0.0)); + delta_M = (double) (sign (delta) == sign (M_diff)); L_prime = langevinD (Q); - // const float denominator = 1 - (c * alpha * (M_s / a) * L_prime); - // - // const float t1_num = (1 - c) * delta_M * M_diff; - // const float t1_den = ((1 - c) * delta * k) - (alpha * M_diff); - // const float t1 = (t1_num / t1_den) * H_d; - // - // const float t2 = c * (M_s / a) * H_d * L_prime; - // - // return (t1 + t2) / denominator; - return H_d * (((nc * delta_M * M_diff) / ((nc * delta * k) - (alpha * M_diff))) + (M_s_oa_tc * L_prime)) / (1.0f - (M_s_oa_tc_talpha * L_prime)); + kap1 = nc * delta_M; + f1Denom = nc * delta * k - alpha * M_diff; + f1 = kap1 * M_diff / f1Denom; + f2 = M_s_oa_tc * L_prime; + f3 = 1.0 - (M_s_oa_tc_talpha * L_prime); + + return H_d * (f1 + f2) / f3; } -float HysteresisProcessing::M_n (float prevM, float k1, float k2, float k3, float k4) +inline double HysteresisProcessing::hysteresisFuncPrime (double H_d, double dMdt) noexcept { - return prevM + (k1 / 6.0f) + (k2 / 3.0f) + (k3 / 3.0f) + (k4 / 6.0f); + const double L_prime2 = langevinD2 (Q); + const double M_diff2 = M_s_oa_talpha * L_prime - 1.0; + + const double f1_p = kap1 * ((M_diff2 / f1Denom) + M_diff * alpha * M_diff2 / (f1Denom * f1Denom)); + const double f2_p = M_s_oaSq_tc_talpha * L_prime2; + const double f3_p = -M_s_oaSq_tc_talphaSq * L_prime2; + + return H_d * (f1_p + f2_p) / f3 - dMdt * f3_p / f3; } -float HysteresisProcessing::process (float H) +inline double HysteresisProcessing::RK2 (double H, double H_d) noexcept { - float H_d = deriv (H, H_n1, H_d_n1); + const double k1 = T * hysteresisFunc (M_n1, H_n1, H_d_n1); + const double k2 = T * hysteresisFunc (M_n1 + (k1 / 2.0), (H + H_n1) / 2.0, (H_d + H_d_n1) / 2.0); - const float k1 = T * hysteresisFunc (M_n1, H_n1, H_d_n1); - const float k2 = T * hysteresisFunc (M_n1 + (k1 / 2.0f), (H + H_n1) / 2.0f, (H_d + H_d_n1) / 2.0f); + return M_n1 + k2; +} + +double HysteresisProcessing::RK4 (double H, double H_d) noexcept +{ + const double H_1_2 = (H + H_n1) / 2.0; + const double H_d_1_2 = (H_d + H_d_n1) / 2.0; - float M = M_n1 + k2; + const double k1 = T * hysteresisFunc (M_n1, H_n1, H_d_n1); + const double k2 = T * hysteresisFunc (M_n1 + (k1 / 2.0), H_1_2, H_d_1_2); + const double k3 = T * hysteresisFunc (M_n1 + (k2 / 2.0), H_1_2, H_d_1_2); + const double k4 = T * hysteresisFunc (M_n1 + k3, H, H_d); - if (std::isnan (M) || abs (M) > 20.0f) + return M_n1 + k1 / 6.0 + k2 / 3.0 + k3 / 3.0 + k4 / 6.0; +} + +inline double HysteresisProcessing::NR (double H, double H_d) noexcept +{ + double M = M_n1; + const double last_dMdt = hysteresisFunc (M_n1, H_n1, H_d_n1); + for (int n = 0; n < numIter; ++n) { - M = 0.0f; - } + const double dMdt = hysteresisFunc (M, H, H_d); + const double dMdtPrime = hysteresisFuncPrime (H_d, dMdt); - M_n1 = M; - H_n1 = H; - H_d_n1 = H_d; + const double deltaNR = (M - M_n1 - Talpha * (dMdt + last_dMdt)) / (1.0 - Talpha * dMdtPrime); + M -= deltaNR; + } return M; } diff --git a/Plugin/Source/Processors/Hysteresis/HysteresisProcessing.h b/Plugin/Source/Processors/Hysteresis/HysteresisProcessing.h @@ -3,6 +3,14 @@ #include "JuceHeader.h" +enum SolverType +{ + RK2 = 0, + RK4, + NR5, + NR10, +}; + /* Hysteresis processing for a model of an analog tape machine. For more information on the DSP happening here, see: @@ -13,50 +21,88 @@ class HysteresisProcessing public: HysteresisProcessing(); + void reset(); + void setSampleRate (double newSR); + + void cook (float drive, float width, float sat, bool v1); + void setSolver (SolverType solverType); + /* Process a single sample */ - float process (float H); + inline double process (double H) noexcept + { + double H_d = deriv (H, H_n1, H_d_n1); + double M = (*this.*solver) (H, H_d); - void reset(); - void setSampleRate (float newSR) { fs = newSR; T = 1.0f / fs; } + if (std::isnan (M) || M > upperLim) + { + M = 0.0; + } + + M_n1 = M; + H_n1 = H; + H_d_n1 = H_d; - void cook (float drive, float width, float sat); + return M; + } private: - /* DEPRECATED (Continued fraction approximation for hyperbolic cotangent) */ - inline float cothApprox (float x); + inline double langevin (double x) const noexcept; // Langevin function + inline double langevinD (double x) const noexcept; // Derivative of Langevin function + inline double langevinD2 (double x) const noexcept; // 2nd derivative of Langevin function + inline double deriv (double x_n, double x_n1, double x_d_n1) const noexcept // Derivative by alpha transform + { + constexpr double dAlpha = 0.9; + return (((1.0 + dAlpha) / T) * (x_n - x_n1)) - dAlpha * x_d_n1; + } + + // hysteresis function dM/dt + inline double hysteresisFunc (double M, double H, double H_d) noexcept; - inline float langevin (float x); // Langevin function - inline float langevinD (float x); // Derivative of Langevin function - inline float deriv (float x_n, float x_n1, float x_d_n1); // Derivative by alpha transform + // derivative of hysteresis func w.r.t M (depends on cached values from computing hysteresisFunc) + inline double hysteresisFuncPrime (double H_d, double dMdt) noexcept; - inline float hysteresisFunc (float M, float H, float H_d); - float M_n (float M_n1, float k1, float k2, float k3, float k4); // DEPRECATED (from RK4 version) + // runge-kutta solvers + inline double RK2 (double H, double H_d) noexcept; + inline double RK4 (double H, double H_d) noexcept; - float fs = 48000.0f; - float T = 1.0f / fs; - float M_s = (float) 1; - float a = M_s / 4.0f; - const float alpha = (float) 1.6e-3; - const float k = 0.47875f; - float c = (float) 1.7e-1; + // newton-raphson solvers + inline double NR (double H, double H_d) noexcept; + int numIter = 0; + + // solver function pointer + double (HysteresisProcessing::*solver) (double, double) = &HysteresisProcessing::NR; + + // parameter values + double fs = 48000.0; + double T = 1.0 / fs; + double Talpha = T / 1.9; + double M_s = 1.0; + double a = M_s / 4.0; + const double alpha = 1.6e-3; + double k = 0.47875; + double c = 1.7e-1; + double upperLim = 20.0; // Save calculations - float nc = 1-c; - float M_s_oa = M_s / a; - float M_s_oa_tc = c * M_s / a; - float M_s_oa_tc_talpha = alpha * c * M_s / a; + double nc = 1-c; + double M_s_oa = M_s / a; + double M_s_oa_talpha = alpha * M_s / a; + double M_s_oa_tc = c * M_s / a; + double M_s_oa_tc_talpha = alpha * c * M_s / a; + double M_s_oaSq_tc_talpha = alpha * c * M_s / (a * a); + double M_s_oaSq_tc_talphaSq = alpha * alpha * c * M_s / (a * a); - float M_n1 = 0.0f; - float H_n1 = 0.0f; - float H_d_n1 = 0.0f; + // state variables + double M_n1 = 0.0; + double H_n1 = 0.0; + double H_d_n1 = 0.0; // temp vars - float Q, M_diff, delta, delta_M, L_prime; - - float coth = 0.0f; + double Q, M_diff, delta, delta_M, L_prime, kap1, f1Denom, f1, f2, f3; + double coth = 0.0; bool nearZero = false; - // JUCE_DECLARE_NONCOPYABLE_WITH_LEAK_DETECTOR (HysteresisProcessing) + JUCE_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR (HysteresisProcessing) }; #endif diff --git a/Plugin/Source/Processors/Hysteresis/HysteresisProcessor.cpp b/Plugin/Source/Processors/Hysteresis/HysteresisProcessor.cpp @@ -2,7 +2,6 @@ enum { - reset = 1048576, numSteps = 500, }; @@ -12,8 +11,7 @@ HysteresisProcessor::HysteresisProcessor (AudioProcessorValueTreeState& vts) satParam = vts.getRawParameterValue ("sat"); widthParam = vts.getRawParameterValue ("width"); osParam = vts.getRawParameterValue ("os"); - - delayParam = vts.getRawParameterValue ("delay_factor"); + modeParam = vts.getRawParameterValue ("mode"); for (int i = 0; i < 5; ++i) overSample[i] = std::make_unique<dsp::Oversampling<float>> @@ -34,11 +32,26 @@ void HysteresisProcessor::createParameterLayout (std::vector<std::unique_ptr<Ran params.push_back (std::make_unique<AudioParameterFloat> ("sat", "Saturation", 0.0f, 1.0f, 0.5f)); params.push_back (std::make_unique<AudioParameterFloat> ("width", "Bias", 0.0f, 1.0f, 0.5f)); - params.push_back (std::make_unique<AudioParameterFloat> ("delay_factor", "Delay", 0.0f, 16.0f, 4.0f)); - + params.push_back (std::make_unique<AudioParameterChoice> ("mode", "Mode", StringArray ({"RK2", "RK4", "NR5", "NR10", "V1"}), 0)); params.push_back (std::make_unique<AudioParameterChoice> ("os", "Oversampling", StringArray ({"1x", "2x", "4x", "8x", "16x"}), 1)); } +void HysteresisProcessor::setSolver (int newSolver) +{ + if (newSolver == 4) // V1 + { + useV1 = true; + newSolver = 1; // RK4 + } + else + { + useV1 = false; + } + + for (int ch = 0; ch < 2; ++ch) + hProcs[ch].setSolver (static_cast<SolverType> (newSolver)); +} + float HysteresisProcessor::calcMakeup() { return (1.0f + 0.6f * width[0].getTargetValue()) / (0.5f + 1.5f * (1.0f - sat[0].getTargetValue())); @@ -70,18 +83,35 @@ void HysteresisProcessor::setSaturation (float newSaturation) } } -void HysteresisProcessor::toggleOnOff (bool shouldBeOn) +void HysteresisProcessor::setOversampling() { - if (shouldBeOn == isOn) - return; + if ((int) *osParam != prevOS) + { + overSamplingFactor = (int) powf(2.0f, *osParam); + prevOS = (int) *osParam; - isChanging = true; + for (int ch = 0; ch < 2; ++ch) + { + hProcs[ch].setSampleRate (fs * overSamplingFactor); + hProcs[ch].cook (drive[ch].getCurrentValue(), width[ch].getCurrentValue(), sat[ch].getCurrentValue(), wasV1); + hProcs[ch].reset(); + } + + calcBiasFreq(); + } +} + +void HysteresisProcessor::calcBiasFreq() +{ + biasFreq = fs * overSamplingFactor / 2.0f; } void HysteresisProcessor::prepareToPlay (double sampleRate, int samplesPerBlock) { fs = (float) sampleRate; overSamplingFactor = (int) powf(2.0f, *osParam); + wasV1 = useV1; + calcBiasFreq(); for (int ch = 0; ch < 2; ++ch) { @@ -90,9 +120,11 @@ void HysteresisProcessor::prepareToPlay (double sampleRate, int samplesPerBlock) sat[ch].skip (numSteps); makeup[ch].skip (numSteps); - hProcs[ch].setSampleRate ((float) (sampleRate * overSamplingFactor)); - hProcs[ch].cook (drive[ch].getCurrentValue(), width[ch].getCurrentValue(), sat[ch].getCurrentValue()); + hProcs[ch].setSampleRate (sampleRate * overSamplingFactor); + hProcs[ch].cook (drive[ch].getCurrentValue(), width[ch].getCurrentValue(), sat[ch].getCurrentValue(), wasV1); hProcs[ch].reset(); + + biasAngle[ch] = 0.0f; } for (int i = 0; i < 5; ++i) @@ -108,11 +140,6 @@ void HysteresisProcessor::prepareToPlay (double sampleRate, int samplesPerBlock) // dcLower[0].calcCoefs (dcShelfFreq, 0.707f, Decibels::decibelsToGain (-12.0f)); // dcLower[1].reset (sampleRate); // dcLower[1].calcCoefs (dcShelfFreq, 0.707f, Decibels::decibelsToGain (-12.0f)); - - fadeBuffer.setSize (2, samplesPerBlock); - resetCount = 0; - - toggleOnOff (true); } void HysteresisProcessor::releaseResources() @@ -129,31 +156,21 @@ float HysteresisProcessor::getLatencySamples() const noexcept void HysteresisProcessor::processBlock (AudioBuffer<float>& buffer, MidiBuffer& /*midi*/) { + setSolver ((int) *modeParam); setDrive (*driveParam); setSaturation (*satParam); setWidth (1.0f - *widthParam); + setOversampling(); - if ((int) *osParam != prevOS) - { - overSamplingFactor = (int) powf(2.0f, *osParam); - prevOS = (int) *osParam; + bool needsSmoothing = drive[0].isSmoothing() || width[0].isSmoothing() || sat[0].isSmoothing() || wasV1 != useV1; + if (useV1 != wasV1) + { for (int ch = 0; ch < 2; ++ch) - { - hProcs[ch].setSampleRate (fs * overSamplingFactor); - hProcs[ch].cook (drive[ch].getCurrentValue(), width[ch].getCurrentValue(), sat[ch].getCurrentValue()); hProcs[ch].reset(); - } } - if (! isOn && ! isChanging) //bypass - return; - - if (isChanging || resetCount + buffer.getNumSamples() >= reset) - { - for (int ch = 0; ch < buffer.getNumChannels(); ++ch) - fadeBuffer.copyFrom (ch, 0, buffer, ch, 0, buffer.getNumSamples()); - } + wasV1 = useV1; // clip input to +9 dB for (int ch = 0; ch < buffer.getNumChannels(); ++ch) @@ -161,87 +178,105 @@ void HysteresisProcessor::processBlock (AudioBuffer<float>& buffer, MidiBuffer& buffer.getWritePointer (ch), -8.0f, 8.0f, buffer.getNumSamples()); dsp::AudioBlock<float> block (buffer); - dsp::AudioBlock<float> osBlock; - - osBlock = overSample[(int) *osParam]->processSamplesUp (block); + dsp::AudioBlock<float> osBlock = overSample[(int) *osParam]->processSamplesUp (block); - float* ptrArray[] = { osBlock.getChannelPointer(0), osBlock.getChannelPointer(1) }; - AudioBuffer<float> osBuffer (ptrArray, 2, static_cast<int> (osBlock.getNumSamples())); - - for (int channel = 0; channel < osBuffer.getNumChannels(); ++channel) + if (needsSmoothing) { - auto* x = osBuffer.getWritePointer (channel); - for (int samp = 0; samp < osBuffer.getNumSamples(); samp++) - { - if (drive[channel].isSmoothing() || width[channel].isSmoothing() || sat[channel].isSmoothing()) - hProcs[channel].cook (drive[channel].getNextValue(), width[channel].getNextValue(), sat[channel].getNextValue()); - - x[samp] = hProcs[channel].process (x[samp]) * makeup[channel].getNextValue(); - } + if (useV1) + processSmoothV1 (osBlock); + else + processSmooth (osBlock); + } + else + { + if (useV1) + processV1 (osBlock); + else + process (osBlock); } overSample[(int) *osParam]->processSamplesDown (block); - resetCount += buffer.getNumSamples(); - if (resetCount >= reset && ! isChanging) + applyDCBlockers (buffer); +} + +void HysteresisProcessor::process (dsp::AudioBlock<float>& block) +{ + for (int channel = 0; channel < block.getNumChannels(); ++channel) { - resetCount = 0; + auto* x = block.getChannelPointer (channel); + for (int samp = 0; samp < block.getNumSamples(); samp++) + { + x[samp] = (float) hProcs[channel].process ((double) x[samp]) * makeup[channel].getNextValue(); + } + } +} - dsp::AudioBlock<float> fadeBlock (fadeBuffer); - dsp::AudioBlock<float> osFadeBlock; +void HysteresisProcessor::processSmooth (dsp::AudioBlock<float>& block) +{ + for (int channel = 0; channel < block.getNumChannels(); ++channel) + { + auto* x = block.getChannelPointer (channel); + for (int samp = 0; samp < block.getNumSamples(); samp++) + { + hProcs[channel].cook (drive[channel].getNextValue(), width[channel].getNextValue(), sat[channel].getNextValue(), false); - osFadeBlock = overSample[(int) *osParam]->processSamplesUp (fadeBlock); + x[samp] = (float) hProcs[channel].process ((double) x[samp]) * makeup[channel].getNextValue(); + } + } +} - float* fadePtrArray[] = { osFadeBlock.getChannelPointer(0), osFadeBlock.getChannelPointer(1) }; - AudioBuffer<float> osFadeBuffer (fadePtrArray, 2, static_cast<int> (osFadeBlock.getNumSamples())); +void HysteresisProcessor::processV1 (dsp::AudioBlock<float>& block) +{ + const auto angleDelta = MathConstants<float>::twoPi * biasFreq / (fs * overSamplingFactor); - for (int channel = 0; channel < osFadeBuffer.getNumChannels(); ++channel) + for (int channel = 0; channel < block.getNumChannels(); ++channel) + { + auto* x = block.getChannelPointer (channel); + for (int samp = 0; samp < block.getNumSamples(); samp++) { - hProcs[channel].reset(); + float bias = biasGain * (1.0f - width[channel].getCurrentValue()) * std::sin (biasAngle[channel]); + biasAngle[channel] += angleDelta; + biasAngle[channel] -= MathConstants<float>::twoPi * (biasAngle[channel] >= MathConstants<float>::twoPi); - auto* x = osFadeBuffer.getWritePointer (channel); - for (int samp = 0; samp < osFadeBuffer.getNumSamples(); samp++) - { - x[samp] = hProcs[channel].process (x[samp]) * makeup[channel].getCurrentValue(); - } + x[samp] = (float) hProcs[channel].process (10000.0 * (double) (x[samp] + bias)); } - overSample[(int) *osParam]->processSamplesDown (fadeBlock); - - buffer.applyGainRamp (0, buffer.getNumSamples(), 1.0f, 0.0f); - - for (int ch = 0; ch < buffer.getNumChannels(); ++ch) - buffer.addFromWithRamp (ch, 0, fadeBuffer.getReadPointer (ch), buffer.getNumSamples(), 0.0f, 1.0f); + FloatVectorOperations::multiply (x, 1.414f / 10000.0f, (int) block.getNumSamples()); } +} - for (int channel = 0; channel < buffer.getNumChannels(); ++channel) - { - // auto* x = buffer.getWritePointer (channel); - // for (int samp = 0; samp < buffer.getNumSamples(); samp++) - // { - // // x[samp] = dcBlocker[channel].processSample (x[samp]); - // // x[samp] = dcLower[channel].processSample (x[samp]); - // } - } +void HysteresisProcessor::processSmoothV1 (dsp::AudioBlock<float>& block) +{ + const auto angleDelta = MathConstants<float>::twoPi * biasFreq / (fs * overSamplingFactor); - if (isChanging) + for (int channel = 0; channel < block.getNumChannels(); ++channel) { - if (! isOn) //fading in + auto* x = block.getChannelPointer (channel); + for (int samp = 0; samp < block.getNumSamples(); samp++) { - buffer.applyGainRamp (0, buffer.getNumSamples(), 0.0f, 1.0f); + hProcs[channel].cook (drive[channel].getNextValue(), width[channel].getNextValue(), sat[channel].getNextValue(), true); - for (int ch = 0; ch < buffer.getNumChannels(); ++ch) - buffer.addFromWithRamp (ch, 0, fadeBuffer.getReadPointer (ch), buffer.getNumSamples(), 1.0f, 0.0f); - } - else if (isOn) //fading out - { - buffer.applyGainRamp (0, buffer.getNumSamples(), 1.0f, 0.0f); + float bias = biasGain * (1.0f - width[channel].getCurrentValue()) * std::sin (biasAngle[channel]); + biasAngle[channel] += angleDelta; + biasAngle[channel] -= MathConstants<float>::twoPi * (biasAngle[channel] >= MathConstants<float>::twoPi); - for (int ch = 0; ch < buffer.getNumChannels(); ++ch) - buffer.addFromWithRamp (ch, 0, fadeBuffer.getReadPointer (ch), buffer.getNumSamples(), 0.0f, 1.0f); + x[samp] = (float) hProcs[channel].process (10000.0 * (double) (x[samp] + bias)); } - isChanging = false; - isOn = ! isOn; + FloatVectorOperations::multiply (x, 1.414f / 10000.0f, (int) block.getNumSamples()); + } +} + +void HysteresisProcessor::applyDCBlockers (AudioBuffer<float>& buffer) +{ + for (int channel = 0; channel < buffer.getNumChannels(); ++channel) + { + auto* x = buffer.getWritePointer (channel); + for (int samp = 0; samp < buffer.getNumSamples(); samp++) + { + x[samp] = dcBlocker[channel].processSample (x[samp]); + // x[samp] = dcLower[channel].processSample (x[samp]); + } } } diff --git a/Plugin/Source/Processors/Hysteresis/HysteresisProcessor.h b/Plugin/Source/Processors/Hysteresis/HysteresisProcessor.h @@ -19,23 +19,30 @@ public: /* Proceess a buffer. */ void processBlock (AudioBuffer<float>& buffer, MidiBuffer& midiBuffer); - void toggleOnOff (bool shouldBeOn); - static void createParameterLayout (std::vector<std::unique_ptr<RangedAudioParameter>>& params); float getLatencySamples() const noexcept; +private: + void setSolver (int newSolver); void setDrive (float newDrive); void setSaturation (float newSat); void setWidth (float newWidth); + void setOversampling(); float calcMakeup(); + void calcBiasFreq(); + + void process (dsp::AudioBlock<float>& block); + void processSmooth (dsp::AudioBlock<float>& block); + void processV1 (dsp::AudioBlock<float>& block); + void processSmoothV1 (dsp::AudioBlock<float>& block); + void applyDCBlockers (AudioBuffer<float>& buffer); -private: std::atomic<float>* driveParam = nullptr; std::atomic<float>* satParam = nullptr; std::atomic<float>* widthParam = nullptr; std::atomic<float>* osParam = nullptr; - std::atomic<float>* delayParam = nullptr; + std::atomic<float>* modeParam = nullptr; SmoothedValue<float, ValueSmoothingTypes::Linear> drive[2]; SmoothedValue<float, ValueSmoothingTypes::Linear> width[2]; @@ -49,15 +56,14 @@ private: TransformerHPF dcBlocker[2]; // TransformerShelf dcLower[2]; - AudioBuffer<float> fadeBuffer; - bool isOn = false; - bool isChanging = true; - int overSamplingFactor = 2; - const float dcFreq = 10.0f; + const float dcFreq = 35.0f; const float dcShelfFreq = 60.0f; - int resetCount = 0; + float biasGain = 10.0f; + float biasFreq = 48000.0f; + float biasAngle[2]; + bool wasV1 = false, useV1 = false; JUCE_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR (HysteresisProcessor) };