jatinchowdhury18 / audio_dspy
Showing 2 of 2 files from the diff.

@@ -22,7 +22,6 @@
 22 22 ` """Class to implement hysteresis processing"""` 23 23 24 24 ` def __init__(self, drive, sat, width, fs, dAlpha=1.0, mode='RK2'):` 25 - ` # def __init__(self, M_s, a, alpha, k, c, fs, dAlpha=1.0, mode='RK2'):` 26 25 ` """` 27 26 ` Parameters` 28 27 ` ----------`
@@ -36,6 +35,20 @@
 36 35 ` Sample rate` 37 36 ` dAlpha : float` 38 37 ` Alpha value used for the alpha transform` 38 + ` mode : {'RK2', 'RK4', 'NR*'}, optional` 39 + ` 'RK2'` 40 + ` By default, this class will use the 2nd-order` 41 + ` Runge-Kutta method for solving the Jiles-Atherton` 42 + ` equation` 43 + 44 + ` 'RK4'` 45 + ` Mode uses the 4th-order Runge-Kutta method` 46 + 47 + ` 'NR{X}'` 48 + ` Mode uses Newton-Raphson iterations, with 'X'` 49 + ` maximum iterations, i.e. 'NR10' corresponds to` 50 + ` a Newton-Raphson iteration that times out after` 51 + ` 10 iterations.` 39 52 ` """` 40 53 ` self.deriv = Differentiator(fs, dAlpha)` 41 54 ` self.T = 1.0 / fs`
@@ -46,23 +59,33 @@
 46 59 ` self.k = 30 * (1-0.5)**6 + 0.01 # Coercivity` 47 60 ` self.c = (1-width)**0.5 - 0.01 # Changes slope` 48 61 49 - ` assert mode == 'RK2' or mode == 'RK4', "Invalid mode!"` 62 + ` assert mode == 'RK2' or mode == 'RK4' or mode[:2] == 'NR', "Invalid mode!"` 50 63 ` self.mode = mode` 51 64 52 - ` def langevin(self, x):` 65 + ` @staticmethod` 66 + ` def langevin(x):` 53 67 ` """Langevin function: coth(x) - (1/x)"""` 54 68 ` if (abs(x) > 10 ** -4):` 55 69 ` return (1 / np.tanh(x)) - (1/x)` 56 70 ` else:` 57 71 ` return (x / 3)` 58 72 59 - ` def langevin_deriv(self, x):` 73 + ` @staticmethod` 74 + ` def langevin_deriv(x):` 60 75 ` """Derivative of the Langevin function: (1/x^2) - coth(x)^2 + 1"""` 61 76 ` if (abs(x) > 10 ** -4):` 62 77 ` return (1 / x ** 2) - (1 / np.tanh(x)) ** 2 + 1` 63 78 ` else:` 64 79 ` return (1 / 3)` 65 80 81 + ` @staticmethod` 82 + ` def langevin_deriv2(x):` 83 + ` """2nd derivative of the Langevin function: 2 coth(x) (coth(x)^2 - 1) - 2/x^3"""` 84 + ` if (abs(x) > 10 ** -3):` 85 + ` return 2 * (1 / np.tanh(x)) * ((1 / np.tanh(x)) ** 2 - 1) - (2 / x ** 3)` 86 + ` else:` 87 + ` return -2 * x / 15` 88 + 66 89 ` def dMdt(self, M, H, H_d):` 67 90 ` """Jiles-Atherton differential equation` 68 91
@@ -123,7 +146,7 @@
 123 146 ` return M_n1 + k2` 124 147 125 148 ` def RK4(self, M_n1, H, H_n1, H_d, H_d_n1):` 126 - ` """Compute hysteresis function with Runge-Kutta 2nd order` 149 + ` """Compute hysteresis function with Runge-Kutta 4th order` 127 150 128 151 ` Parameters` 129 152 ` ----------`
@@ -151,6 +174,79 @@
 151 174 ` k4 = self.T * self.dMdt(M_n1 + k3, H, H_d)` 152 175 ` return M_n1 + (k1 / 6) + (k2 / 3) + (k3 / 3) + (k4 / 6)` 153 176 177 + ` def dMdt_prime(self, M, H, H_d):` 178 + ` """Jiles-Atherton differential equation Jacobian` 179 + 180 + ` Parameters` 181 + ` ----------` 182 + ` M : float` 183 + ` Magnetisation` 184 + ` H : float` 185 + ` Magnetic field` 186 + ` H_d : float` 187 + ` Time derivative of magnetic field` 188 + 189 + ` Returns` 190 + ` -------` 191 + ` dMdt : float` 192 + ` Derivative of dMdt w.r.t. M` 193 + ` """` 194 + ` Q = (H + self.alpha * M) / self.a` 195 + ` M_diff = self.M_s * self.langevin(Q) - M` 196 + ` delta = 1 if H_d > 0 else -1` 197 + ` delta_M = 1 if np.sign(delta) == np.sign(M_diff) else 0` 198 + ` L_prime = self.langevin_deriv(Q)` 199 + ` L_prime2 = self.langevin_deriv2(Q)` 200 + ` M_diff2 = self.alpha * self.M_s * L_prime / self.a - 1` 201 + 202 + ` k1 = (1 - self.c) * delta_M` 203 + ` k2 = (1 - self.c) * delta * self.k` 204 + 205 + ` f1_denom = k2 - self.alpha * M_diff` 206 + ` f1 = H_d * k1 * M_diff / f1_denom` 207 + ` f2 = self.c * self.M_s * H_d * L_prime / self.a` 208 + ` f3 = 1 - self.c * self.alpha * self.M_s * L_prime / self.a` 209 + 210 + ` f1_p = k1 * H_d * ((M_diff2 / f1_denom) +` 211 + ` (M_diff * self.alpha * M_diff2 / f1_denom ** 2))` 212 + ` f2_p = self.c * self.alpha * self.M_s * H_d * L_prime2 / self.a ** 2` 213 + ` f3_p = -self.c * self.alpha ** 2 * self.M_s * L_prime2 / self.a ** 2` 214 + 215 + ` return ((f1_p + f2_p) * f3 - (f1 + f2) * f3_p) / f3 ** 2` 216 + 217 + ` def NR(self, M_n1, H, H_n1, H_d, H_d_n1):` 218 + ` """Compute hysteresis function with Newton-Raphson iteration` 219 + 220 + ` Parameters` 221 + ` ----------` 222 + ` M_n1 : float` 223 + ` Previous magnetisation` 224 + ` H : float` 225 + ` Magnetic field` 226 + ` H_n1 : float` 227 + ` Previous magnetic field` 228 + ` H_d : float` 229 + ` Magnetic field derivative` 230 + ` H_d_n1 : float` 231 + ` Previous magnetic field derivative` 232 + 233 + ` Returns` 234 + ` -------` 235 + ` M : float` 236 + ` Current magnetisation` 237 + ` """` 238 + ` N_iter = int(self.mode[2:])` 239 + 240 + ` M = M_n1` 241 + ` T_2 = self.T / 2.0` 242 + ` last_dMdt = self.dMdt(M_n1, H_n1, H_d_n1)` 243 + ` for _ in range(N_iter):` 244 + ` delta = (M - M_n1 - T_2 * (self.dMdt(M, H, H_d) + last_dMdt)` 245 + ` ) / (1 - T_2 * self.dMdt_prime(M, H, H_d))` 246 + ` M -= delta` 247 + 248 + ` return M` 249 + 154 250 ` def process_block(self, x):` 155 251 ` """Process block of samples"""` 156 252 ` M_out = np.zeros(len(x))`
@@ -158,14 +254,17 @@
 158 254 ` H_n1 = 0` 159 255 ` H_d_n1 = 0` 160 256 257 + ` if self.mode == 'RK2':` 258 + ` solver = self.RK2` 259 + ` elif self.mode == 'RK4':` 260 + ` solver = self.RK4` 261 + ` elif self.mode[:2] == 'NR':` 262 + ` solver = self.NR` 263 + 161 264 ` n = 0` 162 265 ` for H in x:` 163 266 ` H_d = self.deriv.differentiate(H)` 164 - 165 - ` if self.mode == 'RK2':` 166 - ` M = self.RK2(M_n1, H, H_n1, H_d, H_d_n1)` 167 - ` elif self.mode == 'RK4':` 168 - ` M = self.RK4(M_n1, H, H_n1, H_d, H_d_n1)` 267 + ` M = solver(M_n1, H, H_n1, H_d, H_d_n1)` 169 268 170 269 ` M_n1 = M` 171 270 ` H_n1 = H`

@@ -19,5 +19,11 @@
 19 19 ` hysteresis2 = adsp.Hysteresis(1.0, 1.0, 1.0, _FS_, mode='RK2')` 20 20 ` y2 = hysteresis2.process_block(x)` 21 21 22 + ` hysteresis3 = adsp.Hysteresis(1.0, 1.0, 1.0, _FS_, mode='NR10')` 23 + ` y3 = hysteresis3.process_block(x)` 24 + 22 25 ` self.assertTrue(np.sum(np.abs(y - y2)) / _N_ < 5.0e-6,` 23 - ` 'Hysteresis response is incorrect!')` 26 + ` 'RK2 Hysteresis response is incorrect!')` 27 + 28 + ` self.assertTrue(np.sum(np.abs(y - y3)) / _N_ < 5.0e-6,` 29 + ` 'NR10 Hysteresis response is incorrect!')`
Files Coverage
audio_dspy 86.37%
tests 100.00%
setup.py 100.00%
Project Totals (27 files) 91.49%