HysteresisProcessing.h (6146B)
1 #ifndef HYSTERESISPROCESSING_H_INCLUDED 2 #define HYSTERESISPROCESSING_H_INCLUDED 3 4 #include "HysteresisOps.h" 5 #include "HysteresisSTN.h" 6 7 enum SolverType 8 { 9 RK2 = 0, 10 RK4, 11 NR4, 12 NR8, 13 STN, 14 NUM_SOLVERS 15 }; 16 17 /* 18 Hysteresis processing for a model of an analog tape machine. 19 For more information on the DSP happening here, see: 20 https://ccrma.stanford.edu/~jatin/420/tape/TapeModel_DAFx.pdf 21 */ 22 class HysteresisProcessing 23 { 24 public: 25 HysteresisProcessing(); 26 HysteresisProcessing (HysteresisProcessing&&) noexcept = default; 27 28 void reset(); 29 void setSampleRate (double newSR); 30 31 void cook (double drive, double width, double sat, bool v1); 32 33 /* Process a single sample */ 34 template <SolverType solver, typename Float> 35 inline Float process (Float H) noexcept 36 { 37 auto H_d = HysteresisOps::deriv (H, H_n1, H_d_n1, (Float) T); 38 39 Float M; 40 switch (solver) 41 { 42 case RK2: 43 M = RK2Solver (H, H_d); 44 break; 45 case RK4: 46 M = RK4Solver (H, H_d); 47 break; 48 case NR4: 49 M = NRSolver<4> (H, H_d); 50 break; 51 case NR8: 52 M = NRSolver<8> (H, H_d); 53 break; 54 case STN: 55 M = STNSolver (H, H_d); 56 break; 57 58 default: 59 M = 0.0; 60 }; 61 62 // check for instability 63 #if HYSTERESIS_USE_SIMD 64 auto notIllCondition = ! (xsimd::isnan (M) || (M > upperLim)); 65 M = xsimd::select (notIllCondition, M, (Float) 0.0); 66 H_d = xsimd::select (notIllCondition, H_d, (Float) 0.0); 67 #else 68 bool illCondition = std::isnan (M) || M > upperLim; 69 M = illCondition ? 0.0 : M; 70 H_d = illCondition ? 0.0 : H_d; 71 #endif 72 73 M_n1 = M; 74 H_n1 = H; 75 H_d_n1 = H_d; 76 77 return M; 78 } 79 80 private: 81 // runge-kutta solvers 82 template <typename Float> 83 inline Float RK2Solver (Float H, Float H_d) noexcept 84 { 85 const Float k1 = HysteresisOps::hysteresisFunc (M_n1, H_n1, H_d_n1, hpState) * T; 86 const Float k2 = HysteresisOps::hysteresisFunc (M_n1 + (k1 * 0.5), (H + H_n1) * 0.5, (H_d + H_d_n1) * 0.5, hpState) * T; 87 88 return M_n1 + k2; 89 } 90 91 template <typename Float> 92 inline Float RK4Solver (Float H, Float H_d) noexcept 93 { 94 const Float H_1_2 = (H + H_n1) * 0.5; 95 const Float H_d_1_2 = (H_d + H_d_n1) * 0.5; 96 97 const Float k1 = HysteresisOps::hysteresisFunc (M_n1, H_n1, H_d_n1, hpState) * T; 98 const Float k2 = HysteresisOps::hysteresisFunc (M_n1 + (k1 * 0.5), H_1_2, H_d_1_2, hpState) * T; 99 const Float k3 = HysteresisOps::hysteresisFunc (M_n1 + (k2 * 0.5), H_1_2, H_d_1_2, hpState) * T; 100 const Float k4 = HysteresisOps::hysteresisFunc (M_n1 + k3, H, H_d, hpState) * T; 101 102 constexpr double oneSixth = 1.0 / 6.0; 103 constexpr double oneThird = 1.0 / 3.0; 104 return M_n1 + k1 * oneSixth + k2 * oneThird + k3 * oneThird + k4 * oneSixth; 105 } 106 107 // newton-raphson solvers 108 template <int nIterations, typename Float> 109 inline Float NRSolver (Float H, Float H_d) noexcept 110 { 111 using namespace chowdsp::SIMDUtils; 112 113 Float M = M_n1; 114 const Float last_dMdt = HysteresisOps::hysteresisFunc (M_n1, H_n1, H_d_n1, hpState); 115 116 Float dMdt; 117 Float dMdtPrime; 118 Float deltaNR; 119 for (int n = 0; n < nIterations; ++n) 120 { 121 using namespace HysteresisOps; 122 dMdt = hysteresisFunc (M, H, H_d, hpState); 123 dMdtPrime = hysteresisFuncPrime (H_d, dMdt, hpState); 124 deltaNR = (M - M_n1 - (Float) Talpha * (dMdt + last_dMdt)) / (Float (1.0) - (Float) Talpha * dMdtPrime); 125 M -= deltaNR; 126 } 127 128 return M; 129 } 130 131 // state transition network solver 132 template <typename Float> 133 inline Float STNSolver (Float H, Float H_d) noexcept 134 { 135 #if HYSTERESIS_USE_SIMD 136 double H_arr alignas (xsimd::default_arch::alignment())[2]; 137 double H_d_arr alignas (xsimd::default_arch::alignment())[2]; 138 double H_n1_arr alignas (xsimd::default_arch::alignment())[2]; 139 double H_d_n1_arr alignas (xsimd::default_arch::alignment())[2]; 140 double M_n1_arr alignas (xsimd::default_arch::alignment())[2]; 141 double M_out alignas (xsimd::default_arch::alignment())[2]; 142 143 H.store_aligned ((double*) H_arr); 144 H_d.store_aligned ((double*) H_d_arr); 145 H_n1.store_aligned ((double*) H_n1_arr); 146 H_d_n1.store_aligned ((double*) H_d_n1_arr); 147 M_n1.store_aligned ((double*) M_n1_arr); 148 149 for (int ch = 0; ch < 2; ++ch) 150 { 151 double input alignas (xsimd::default_arch::alignment())[5] = { H_arr[ch], H_d_arr[ch], H_n1_arr[ch], H_d_n1_arr[ch], M_n1_arr[ch] }; 152 153 // scale derivatives 154 input[1] *= HysteresisSTN::diffMakeup; 155 input[3] *= HysteresisSTN::diffMakeup; 156 FloatVectorOperations::multiply (input, 0.7071 / hpState.a, 4); // scale by drive param 157 158 M_out[ch] = hysteresisSTN.process (input) + M_n1_arr[ch]; 159 } 160 161 return Float::load_aligned (M_out); 162 163 #else 164 double input alignas (xsimd::default_arch::alignment())[5] = { H, H_d, H_n1, H_d_n1, M_n1 }; 165 166 // scale derivatives 167 input[1] *= HysteresisSTN::diffMakeup; 168 input[3] *= HysteresisSTN::diffMakeup; 169 FloatVectorOperations::multiply (input, 0.7071 / hpState.a, 4); // scale by drive param 170 171 return hysteresisSTN.process (input) + M_n1; 172 #endif 173 } 174 175 // parameter values 176 double fs = 48000.0; 177 double T = 1.0 / fs; 178 double Talpha = T / 1.9; 179 double upperLim = 20.0; 180 181 // state variables 182 #if HYSTERESIS_USE_SIMD 183 xsimd::batch<double> M_n1 = 0.0; 184 xsimd::batch<double> H_n1 = 0.0; 185 xsimd::batch<double> H_d_n1 = 0.0; 186 #else 187 double M_n1 = 0.0; 188 double H_n1 = 0.0; 189 double H_d_n1 = 0.0; 190 #endif 191 192 HysteresisSTN hysteresisSTN; 193 HysteresisOps::HysteresisState hpState; 194 195 JUCE_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR (HysteresisProcessing) 196 }; 197 198 #endif