AnalogTapeModel

commit d5a4d1ef176377a10e1ab1802c11de25de936d0c
parent ece26bddf60f1029f237569947906baa09fa30c5
Author: jatinchowdhury18 <jatinchowdhury18@users.noreply.github.com>
Date:   Wed, 27 Feb 2019 09:39:23 -0800

Optimize hysteresis code

Diffstat:
M
Plugin/Source/Processors/Hysteresis/HysteresisProcessing.cpp
 | 
40
++++++++++++++++++++++++----------------
M
Plugin/Source/Processors/Hysteresis/HysteresisProcessing.h
 | 
8
+++++++-
M
Plugin/Source/Processors/Hysteresis/HysteresisProcessor.cpp
 | 
5
+++--

3 files changed, 34 insertions(+), 19 deletions(-)
diff --git a/Plugin/Source/Processors/Hysteresis/HysteresisProcessing.cpp b/Plugin/Source/Processors/Hysteresis/HysteresisProcessing.cpp
@@ -6,21 +6,25 @@ HysteresisProcessing::HysteresisProcessing()
 
 }
 
+float HysteresisProcessing::tanhApprox (float x)
+{
+    const auto xSquare = x * x;
+
+    return x / (1.0f + (xSquare / (3.0f + (xSquare / (5.0f + (xSquare / (7.0f)))))));
+}
+
 float HysteresisProcessing::langevin (float x)
 {
-    if (std::abs (x) > (float) (10e-4))
-        return (1.0f / (float) tanh (x)) - (1.0f / x);
+    if (nearZero)
+        return (tanhRecip) - (1.0f / x);
     else
         return (x / 3.0f);
 }
 
 float HysteresisProcessing::langevinD (float x)
 {
-    if (std::abs (x) > (float) (10e-4))
-    {
-        float tanhRecip = 1.0f / tanh (x);
+    if (nearZero)
         return (1.0f / (x * x)) - (tanhRecip * tanhRecip) + 1.0f;
-    }
     else
         return (1.0f / 3.0f);
 }
@@ -33,6 +37,10 @@ float HysteresisProcessing::deriv (float x_n, float x_n1, float x_d_n1)
 float HysteresisProcessing::hysteresisFunc (float M, float H, float H_d)
 {
     const float Q = (H + alpha * M) / a;
+    tanHyp = (float) tanhApprox (Q);
+    tanhRecip = 1.0f / tanHyp;
+    nearZero = std::abs (Q) > (float) (10e-4);
+
     const float M_diff = M_s * langevin (Q) - M;
 
     const float delta = H_d > 0 ? 1.0f : -1.0f;
@@ -40,15 +48,16 @@ float HysteresisProcessing::hysteresisFunc (float M, float H, float H_d)
 
     const float 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;
+    //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 * (((1 - c) * delta_M * M_diff) / (((1 - c) * delta * k) - (alpha * M_diff))) + (c * (M_s / a) * H_d * L_prime)) / (1 - (c * alpha * (M_s / a) * L_prime));
 }
 
 float HysteresisProcessing::M_n (float prevM, float k1, float k2, float k3, float k4)
@@ -60,7 +69,6 @@ float HysteresisProcessing::process (float H)
 {
     const float H_d = deriv (H, H_n1, H_d_n1);
 
-    const float T = (1.0f / fs);
     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);
     const float k3 = T * hysteresisFunc (M_n1 + (k2 / 2.0f), (H + H_n1) / 2.0f, (H_d + H_d_n1) / 2.0f);
diff --git a/Plugin/Source/Processors/Hysteresis/HysteresisProcessing.h b/Plugin/Source/Processors/Hysteresis/HysteresisProcessing.h
@@ -10,9 +10,10 @@ public:
 
     float process (float H);
 
-    void setSampleRate (float newSR) { fs = newSR; }
+    void setSampleRate (float newSR) { fs = newSR; T = 1.0f / fs; }
 
 private:
+    float tanhApprox (float x);
 
     float langevin (float x);
     float langevinD (float x);
@@ -22,6 +23,7 @@ private:
     float M_n (float M_n1, float k1, float k2, float k3, float k4);
 
     float fs = 48000.0f;
+    float T = 1.0f / fs;
     const float M_s = 350000.0f;
     const float a = (float) 2.2e4;
     const float alpha = (float) 1.6e-3;
@@ -32,6 +34,10 @@ private:
     float H_n1 = 0.0f;
     float H_d_n1 = 0.0f;
 
+    float tanHyp = 0.0f;
+    float tanhRecip = 0.0f;
+    bool nearZero = false;
+
     //JUCE_DECLARE_NONCOPYABLE_WITH_LEAK_DETECTOR (HysteresisProcessing)
 };
 
diff --git a/Plugin/Source/Processors/Hysteresis/HysteresisProcessor.cpp b/Plugin/Source/Processors/Hysteresis/HysteresisProcessor.cpp
@@ -67,8 +67,7 @@ void HysteresisProcessor::processBlock (AudioBuffer<float>& buffer, MidiBuffer& 
         for (int samp = 0; samp < osBuffer.getNumSamples(); samp++)
         {
             //x[samp] = hProcs[channel].process ((float) 5e4 * (x[samp]));
-            
-            x[samp] = 10.0f * hProcs[channel].process ((float) 5e3 * (x[samp] + biasGain * sinf (currentAngle[channel])));
+            x[samp] = hProcs[channel].process ((float) 5e3 * (x[samp] + biasGain * sinf (currentAngle[channel])));
             //x[samp] = biasGain * sinf (currentAngle[channel]);
 
             currentAngle[channel] += angleDelta;
@@ -85,6 +84,8 @@ void HysteresisProcessor::processBlock (AudioBuffer<float>& buffer, MidiBuffer& 
         overSample4->processSamplesDown (block);
     else
         overSample2->processSamplesDown (block);
+
+    buffer.applyGain (10.0f);
 }
 
 void HysteresisProcessor::updateAngleDelta()