HysteresisOps.h (5031B)
1 #pragma once 2 3 #include <JuceHeader.h> 4 #include <cmath> 5 6 #define HYSTERESIS_USE_SIMD 1 7 8 namespace HysteresisOps 9 { 10 using namespace chowdsp::SIMDUtils; 11 12 struct HysteresisState 13 { 14 // parameter values 15 double M_s = 1.0; 16 double a = M_s / 4.0; 17 static constexpr double alpha = 1.6e-3; 18 double k = 0.47875; 19 double c = 1.7e-1; 20 21 // Save calculations 22 double nc = 1 - c; 23 double M_s_oa = M_s / a; 24 double M_s_oa_talpha = alpha * M_s / a; 25 double M_s_oa_tc = c * M_s / a; 26 double M_s_oa_tc_talpha = alpha * c * M_s / a; 27 double M_s_oaSq_tc_talpha = alpha * c * M_s / (a * a); 28 double M_s_oaSq_tc_talphaSq = alpha * alpha * c * M_s / (a * a); 29 30 // temp vars 31 #if HYSTERESIS_USE_SIMD 32 xsimd::batch<double> Q, M_diff, L_prime, kap1, f1Denom, f1, f2, f3; 33 xsimd::batch<double> coth = 0.0; 34 xsimd::batch_bool<double> nearZero; 35 xsimd::batch<double> oneOverQ, oneOverQSq, oneOverQCubed, cothSq, oneOverF3, oneOverF1Denom; 36 #else 37 double Q, M_diff, L_prime, kap1, f1Denom, f1, f2, f3; 38 double coth = 0.0; 39 bool nearZero = false; 40 double oneOverQ, oneOverQSq, oneOverQCubed, cothSq, oneOverF3, oneOverF1Denom; 41 #endif 42 }; 43 44 constexpr double ONE_THIRD = 1.0 / 3.0; 45 constexpr double NEG_TWO_OVER_15 = -2.0 / 15.0; 46 47 constexpr inline int sign (double x) 48 { 49 return int (x > 0.0) - int (x < 0.0); 50 } 51 52 /** Langevin function */ 53 template <typename Float> 54 static inline Float langevin (const HysteresisState& hp) noexcept 55 { 56 #if HYSTERESIS_USE_SIMD 57 return xsimd::select (hp.nearZero, hp.Q * ONE_THIRD, hp.coth - hp.oneOverQ); 58 #else 59 return ! hp.nearZero ? (hp.coth) - (hp.oneOverQ) : hp.Q * ONE_THIRD; 60 #endif 61 } 62 63 /** Derivative of Langevin function */ 64 template <typename Float> 65 static inline Float langevinD (const HysteresisState& hp) noexcept 66 { 67 #if HYSTERESIS_USE_SIMD 68 return xsimd::select (hp.nearZero, (Float) ONE_THIRD, hp.oneOverQSq - hp.cothSq + (Float) 1.0); 69 #else 70 return ! hp.nearZero ? hp.oneOverQSq - hp.cothSq + 1.0 : ONE_THIRD; 71 #endif 72 } 73 74 /** 2nd derivative of Langevin function */ 75 template <typename Float> 76 static inline Float langevinD2 (const HysteresisState& hp) noexcept 77 { 78 #if HYSTERESIS_USE_SIMD 79 return xsimd::select (hp.nearZero, 80 NEG_TWO_OVER_15 * hp.Q, 81 (Float) 2.0 * hp.coth * (hp.cothSq - (Float) 1.0) - ((Float) 2.0 * hp.oneOverQCubed)); 82 #else 83 return ! hp.nearZero 84 ? 2.0 * hp.coth * (hp.cothSq - 1.0) - (2.0 * hp.oneOverQCubed) 85 : NEG_TWO_OVER_15 * hp.Q; 86 #endif 87 } 88 89 /** Derivative by alpha transform */ 90 template <typename Float> 91 static inline Float deriv (Float x_n, Float x_n1, Float x_d_n1, Float T) noexcept 92 { 93 const Float dAlpha = 0.75; 94 return ((((Float) 1.0 + dAlpha) / T) * (x_n - x_n1)) - dAlpha * x_d_n1; 95 } 96 97 /** hysteresis function dM/dt */ 98 template <typename Float> 99 static inline Float hysteresisFunc (Float M, Float H, Float H_d, HysteresisState& hp) noexcept 100 { 101 hp.Q = (H + M * HysteresisState::alpha) * (1.0 / hp.a); 102 hp.oneOverQ = (Float) 1.0 / hp.Q; 103 hp.oneOverQSq = hp.oneOverQ * hp.oneOverQ; 104 hp.oneOverQCubed = hp.oneOverQ * hp.oneOverQSq; 105 106 #if HYSTERESIS_USE_SIMD 107 hp.coth = (Float) 1.0 / xsimd::tanh (hp.Q); 108 hp.nearZero = (hp.Q < 0.001) && (hp.Q > -0.001); 109 #else 110 hp.coth = 1.0 / std::tanh (hp.Q); 111 hp.nearZero = hp.Q < 0.001 && hp.Q > -0.001; 112 #endif 113 114 hp.cothSq = hp.coth * hp.coth; 115 hp.M_diff = langevin<Float> (hp) * hp.M_s - M; 116 117 #if HYSTERESIS_USE_SIMD 118 const auto delta = xsimd::select (H_d >= 0.0, (Float) 1, (Float) -1); 119 const auto delta_M = chowdsp::Math::sign (delta) == chowdsp::Math::sign (hp.M_diff); 120 hp.kap1 = xsimd::select (delta_M, (Float) hp.nc, (Float) 0); 121 #else 122 const auto delta = (Float) ((H_d >= 0.0) - (H_d < 0.0)); 123 const auto delta_M = (Float) (sign (delta) == sign (hp.M_diff)); 124 hp.kap1 = (Float) hp.nc * delta_M; 125 #endif 126 127 hp.L_prime = langevinD<Float> (hp); 128 129 hp.f1Denom = ((Float) hp.nc * delta) * hp.k - (Float) HysteresisState::alpha * hp.M_diff; 130 hp.oneOverF1Denom = (Float) 1.0 / hp.f1Denom; 131 hp.f1 = hp.kap1 * hp.M_diff / hp.f1Denom; 132 hp.f2 = hp.L_prime * hp.M_s_oa_tc; 133 hp.f3 = (Float) 1.0 - (hp.L_prime * hp.M_s_oa_tc_talpha); 134 hp.oneOverF3 = (Float) 1.0 / hp.f3; 135 136 return H_d * (hp.f1 + hp.f2) * hp.oneOverF3; 137 } 138 139 // derivative of hysteresis func w.r.t M (depends on cached values from computing hysteresisFunc) 140 template <typename Float> 141 static inline Float hysteresisFuncPrime (Float H_d, Float dMdt, HysteresisState& hp) noexcept 142 { 143 const Float L_prime2 = langevinD2<Float> (hp); 144 const Float M_diff2 = hp.L_prime * hp.M_s_oa_talpha - 1.0; 145 146 Float f1_p = (M_diff2 * hp.oneOverF1Denom); 147 f1_p += hp.M_diff * HysteresisState::alpha * M_diff2 * (hp.oneOverF1Denom * hp.oneOverF1Denom); 148 f1_p *= hp.kap1; 149 const Float f2_p = L_prime2 * hp.M_s_oaSq_tc_talpha; 150 const Float f3_p = L_prime2 * (-hp.M_s_oaSq_tc_talphaSq); 151 152 return (H_d * (f1_p + f2_p) - dMdt * f3_p) * hp.oneOverF3; 153 } 154 155 } // namespace HysteresisOps