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

@@ -11,6 +11,7 @@
Loading
11 11
from .adaptive_filt import *
12 12
from .level_detector import level_detect
13 13
from .farina import Farina
14 +
from .hysteresis import Hysteresis
14 15
name = "audio_dspy"
15 16
16 17

@@ -0,0 +1,177 @@
Loading
1 +
import numpy as np
2 +
3 +
4 +
class Differentiator:
5 +
    """Time domain differentiation using the alpha transform"""
6 +
7 +
    def __init__(self, fs, alpha=1.0):
8 +
        self.T = 1.0 / fs
9 +
        self.alpha = alpha
10 +
        self.x_1 = 0.0
11 +
        self.xD_1 = 0.0
12 +
13 +
    def differentiate(self, x):
14 +
        xD = (((1 + self.alpha) / self.T) * (x - self.x_1)) - \
15 +
            self.alpha * self.xD_1
16 +
        self.x_1 = x
17 +
        self.xD_1 = xD
18 +
        return xD
19 +
20 +
21 +
class Hysteresis:
22 +
    """Class to implement hysteresis processing"""
23 +
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 +
        """
27 +
        Parameters
28 +
        ----------
29 +
        drive : float
30 +
            Hysteresis drive parameter
31 +
        sat : float
32 +
            Saturation parameter
33 +
        width : float
34 +
            Hysteresis width parameter
35 +
        fs : float
36 +
            Sample rate
37 +
        dAlpha : float
38 +
            Alpha value used for the alpha transform
39 +
        """
40 +
        self.deriv = Differentiator(fs, dAlpha)
41 +
        self.T = 1.0 / fs
42 +
43 +
        self.M_s = 0.5 + 1.5*(1-sat)  # saturation
44 +
        self.a = self.M_s / (0.01 + 6*drive)  # adjustable parameter
45 +
        self.alpha = 1.6e-3
46 +
        self.k = 30 * (1-0.5)**6 + 0.01  # Coercivity
47 +
        self.c = (1-width)**0.5 - 0.01  # Changes slope
48 +
49 +
        assert mode == 'RK2' or mode == 'RK4', "Invalid mode!"
50 +
        self.mode = mode
51 +
52 +
    def langevin(self, x):
53 +
        """Langevin function: coth(x) - (1/x)"""
54 +
        if (abs(x) > 10 ** -4):
55 +
            return (1 / np.tanh(x)) - (1/x)
56 +
        else:
57 +
            return (x / 3)
58 +
59 +
    def langevin_deriv(self, x):
60 +
        """Derivative of the Langevin function: (1/x^2) - coth(x)^2 + 1"""
61 +
        if (abs(x) > 10 ** -4):
62 +
            return (1 / x ** 2) - (1 / np.tanh(x)) ** 2 + 1
63 +
        else:
64 +
            return (1 / 3)
65 +
66 +
    def dMdt(self, M, H, H_d):
67 +
        """Jiles-Atherton differential equation
68 +
69 +
        Parameters
70 +
        ----------
71 +
        M : float
72 +
            Magnetisation
73 +
        H : float
74 +
            Magnetic field
75 +
        H_d : float
76 +
            Time derivative of magnetic field
77 +
78 +
        Returns
79 +
        -------
80 +
        dMdt : float
81 +
            Derivative of magnetisation w.r.t time
82 +
        """
83 +
        Q = (H + self.alpha * M) / self.a
84 +
        M_diff = self.M_s * self.langevin(Q) - M
85 +
        delta = 1 if H_d > 0 else -1
86 +
        delta_M = 1 if np.sign(delta) == np.sign(M_diff) else 0
87 +
        L_prime = self.langevin_deriv(Q)
88 +
89 +
        denominator = 1 - self.c * self.alpha * (self.M_s / self.a) * L_prime
90 +
91 +
        t1_num = (1 - self.c) * delta_M * M_diff
92 +
        t1_den = (1 - self.c) * delta * self.k - self.alpha * M_diff
93 +
        t1 = (t1_num / t1_den) * H_d
94 +
95 +
        t2 = self.c * (self.M_s / self.a) * H_d * L_prime
96 +
97 +
        return (t1 + t2) / denominator
98 +
99 +
    def RK2(self, M_n1, H, H_n1, H_d, H_d_n1):
100 +
        """Compute hysteresis function with Runge-Kutta 2nd order
101 +
102 +
        Parameters
103 +
        ----------
104 +
        M_n1 : float
105 +
            Previous magnetisation
106 +
        H : float
107 +
            Magnetic field
108 +
        H_n1 : float
109 +
            Previous magnetic field
110 +
        H_d : float
111 +
            Magnetic field derivative
112 +
        H_d_n1 : float
113 +
            Previous magnetic field derivative
114 +
115 +
        Returns
116 +
        -------
117 +
        M : float
118 +
            Current magnetisation
119 +
        """
120 +
        k1 = self.T * self.dMdt(M_n1, H_n1, H_d_n1)
121 +
        k2 = self.T * self.dMdt(M_n1 + k1/2, (H + H_n1) /
122 +
                                2, (H_d + H_d_n1) / 2)
123 +
        return M_n1 + k2
124 +
125 +
    def RK4(self, M_n1, H, H_n1, H_d, H_d_n1):
126 +
        """Compute hysteresis function with Runge-Kutta 2nd order
127 +
128 +
        Parameters
129 +
        ----------
130 +
        M_n1 : float
131 +
            Previous magnetisation
132 +
        H : float
133 +
            Magnetic field
134 +
        H_n1 : float
135 +
            Previous magnetic field
136 +
        H_d : float
137 +
            Magnetic field derivative
138 +
        H_d_n1 : float
139 +
            Previous magnetic field derivative
140 +
141 +
        Returns
142 +
        -------
143 +
        M : float
144 +
            Current magnetisation
145 +
        """
146 +
        k1 = self.T * self.dMdt(M_n1, H_n1, H_d_n1)
147 +
        k2 = self.T * self.dMdt(M_n1 + k1/2, (H + H_n1) /
148 +
                                2, (H_d + H_d_n1) / 2)
149 +
        k3 = self.T * self.dMdt(M_n1 + k2/2, (H + H_n1) /
150 +
                                2, (H_d + H_d_n1) / 2)
151 +
        k4 = self.T * self.dMdt(M_n1 + k3, H, H_d)
152 +
        return M_n1 + (k1 / 6) + (k2 / 3) + (k3 / 3) + (k4 / 6)
153 +
154 +
    def process_block(self, x):
155 +
        """Process block of samples"""
156 +
        M_out = np.zeros(len(x))
157 +
        M_n1 = 0
158 +
        H_n1 = 0
159 +
        H_d_n1 = 0
160 +
161 +
        n = 0
162 +
        for H in x:
163 +
            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)
169 +
170 +
            M_n1 = M
171 +
            H_n1 = H
172 +
            H_d_n1 = H_d
173 +
174 +
            M_out[n] = M
175 +
            n += 1
176 +
177 +
        return M_out

@@ -0,0 +1,23 @@
Loading
1 +
from unittest import TestCase
2 +
import numpy as np
3 +
4 +
import audio_dspy as adsp
5 +
6 +
_FS_ = 44100.0
7 +
_N_ = 1024
8 +
9 +
10 +
class TestHysteresis(TestCase):
11 +
    def test_hysteresis(self):
12 +
        freq = 100
13 +
        n = np.arange(_N_)
14 +
        x = np.sin(2 * np.pi * n * freq / _FS_) * (n/_N_)
15 +
16 +
        hysteresis = adsp.Hysteresis(1.0, 1.0, 1.0, _FS_, mode='RK4')
17 +
        y = hysteresis.process_block(x)
18 +
19 +
        hysteresis2 = adsp.Hysteresis(1.0, 1.0, 1.0, _FS_, mode='RK2')
20 +
        y2 = hysteresis2.process_block(x)
21 +
22 +
        self.assertTrue(np.sum(np.abs(y - y2)) / _N_ < 5.0e-6,
23 +
                        'Hysteresis response is incorrect!')
Files Coverage
audio_dspy 85.80%
tests 100.00%
setup.py 100.00%
Project Totals (27 files) 91.25%
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading