jatinchowdhury18 / audio_dspy

Compare b92cf98 ... +0 ... 5d39eba


@@ -22,7 +22,6 @@
Loading
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 @@
Loading
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 @@
Loading
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 @@
Loading
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,21 +174,97 @@
Loading
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))
157 253
        M_n1 = 0
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 @@
Loading
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!')

Everything is accounted for!

No changes detected that need to be reviewed.
What changes does Codecov check for?
Lines, not adjusted in diff, that have changed coverage data.
Files that introduced coverage data that had none before.
Files that have missing coverage data that once were tracked.
Files Coverage
audio_dspy 0.57% 86.37%
tests 100.00%
setup.py 100.00%
Project Totals (27 files) 91.49%
Loading