瀏覽代碼

Euler approach to harmonic noise cancel

tags/1.6.1
Trevor Irons 5 年之前
父節點
當前提交
627dc90d14
共有 3 個檔案被更改,包括 91 行新增27 行删除
  1. 1
    1
      akvo/gui/akvoGUI.py
  2. 89
    25
      akvo/tressel/harmonic.py
  3. 1
    1
      akvo/tressel/mrsurvey.py

+ 1
- 1
akvo/gui/akvoGUI.py 查看文件

@@ -1322,7 +1322,7 @@ def main():
1322 1322
     logo2 = pkg_resources.resource_filename(__name__, 'akvo2.png')
1323 1323
     qApp = QtWidgets.QApplication(sys.argv)
1324 1324
 
1325
-    ssplash = False
1325
+    ssplash = True
1326 1326
     if ssplash:
1327 1327
         pixmap = QtGui.QPixmap(logo)
1328 1328
         splash = QtWidgets.QSplashScreen(pixmap, QtCore.Qt.WindowStaysOnTopHint)

+ 89
- 25
akvo/tressel/harmonic.py 查看文件

@@ -2,6 +2,7 @@ import numpy as np
2 2
 from scipy.optimize import least_squares 
3 3
 from scipy.optimize import minimize
4 4
 from scipy.linalg import lstsq as sclstsq
5
+import scipy.linalg as lin
5 6
 
6 7
 def harmonic2 ( f1, f2, sN, fs, nK, t ): 
7 8
     """
@@ -14,32 +15,32 @@ def harmonic2 ( f1, f2, sN, fs, nK, t ):
14 15
         nK = number of harmonics to calculate 
15 16
         t = time samples 
16 17
     """
17
-    print("building matrix ")
18
+    print("building matrix 2")
18 19
     A = np.zeros( (len(t),  4*nK) )
19 20
     #f1 = f1MHz * 1e-3
20 21
     #f2 = f2MHz * 1e-3
21 22
     for irow, tt in enumerate(t): 
22
-        A[irow, 0:2*nK:2] = np.cos( np.arange(nK)*2*np.pi*(f1/fs)*irow )
23
-        A[irow, 1:2*nK:2] = np.sin( np.arange(nK)*2*np.pi*(f1/fs)*irow )
24
-        A[irow, 2*nK::2] = np.cos( np.arange(nK)*2*np.pi*(f2/fs)*irow )
23
+        A[irow, 0:2*nK:2]  = np.cos( np.arange(nK)*2*np.pi*(f1/fs)*irow )
24
+        A[irow, 1:2*nK:2]  = np.sin( np.arange(nK)*2*np.pi*(f1/fs)*irow )
25
+        A[irow, 2*nK::2]   = np.cos( np.arange(nK)*2*np.pi*(f2/fs)*irow )
25 26
         A[irow, 2*nK+1::2] = np.sin( np.arange(nK)*2*np.pi*(f2/fs)*irow )
26 27
 
27 28
     v = np.linalg.lstsq(A, sN, rcond=1e-8)
28 29
     #v = sclstsq(A, sN) #, rcond=1e-6)
29 30
 
30
-    alpha = v[0][0:2*nK:2]
31
-    beta  = v[0][1:2*nK:2]
31
+    alpha = v[0][0:2*nK:2] + 1e-16
32
+    beta  = v[0][1:2*nK:2] + 1e-16
32 33
     amp = np.sqrt( alpha**2 + beta**2 )
33 34
     phase = np.arctan(- beta/alpha)
34 35
     
35
-    alpha2 = v[0][2*nK::2]
36
-    beta2  = v[0][2*nK+1::2]
36
+    alpha2 = v[0][2*nK::2]   + 1e-16
37
+    beta2  = v[0][2*nK+1::2] + 1e-16
37 38
     amp2 = np.sqrt( alpha2**2 + beta2**2 )
38 39
     phase2 = np.arctan(- beta2/alpha2)
39 40
 
40 41
     h = np.zeros(len(t))
41 42
     for ik in range(nK):
42
-        h += np.sqrt(alpha[ik]**2 + beta[ik]**2) * np.cos( 2.*np.pi*ik * (f1/fs) * np.arange(0, len(t), 1 )  + phase[ik] ) \
43
+        h += np.sqrt( alpha[ik]**2 + beta[ik]**2)  * np.cos( 2.*np.pi*ik * (f1/fs) * np.arange(0, len(t), 1 )  + phase[ik] ) \
43 44
            + np.sqrt(alpha2[ik]**2 + beta2[ik]**2) * np.cos( 2.*np.pi*ik * (f2/fs) * np.arange(0, len(t), 1 )  + phase2[ik] )
44 45
 
45 46
     return sN-h
@@ -60,10 +61,17 @@ def harmonic ( f0, sN, fs, nK, t ):
60 61
         A[irow, 0::2] = np.cos( np.arange(nK)*2*np.pi*(f0/fs)*irow )
61 62
         A[irow, 1::2] = np.sin( np.arange(nK)*2*np.pi*(f0/fs)*irow )
62 63
 
63
-    v = np.linalg.lstsq(A, sN, rcond=None) #, rcond=1e-8)
64
-
64
+    v = np.linalg.lstsq(A, sN, rcond=1e-16) # rcond=None) #, rcond=1e-8)
65
+    #v = sclstsq(A, sN) #, rcond=1e-6)
65 66
     alpha = v[0][0::2]
66 67
     beta  = v[0][1::2]
68
+    
69
+    #print("Solving A A.T")
70
+    #v = lin.solve(np.dot(A,A.T).T, sN) #, rcond=1e-6)
71
+    #v = np.dot(A.T, v)
72
+    #v = np.dot(np.linalg.inv(np.dot(A.T, A)), np.dot(A.T, sN))
73
+    #alpha = v[0::2]
74
+    #beta  = v[1::2]
67 75
 
68 76
     amp = np.sqrt( alpha**2 + beta**2 )
69 77
     phase = np.arctan(- beta/alpha)
@@ -87,6 +95,59 @@ def harmonic ( f0, sN, fs, nK, t ):
87 95
     #plt.title("modelled noise")
88 96
     return sN-h
89 97
 
98
+
99
+def harmonicEuler ( f0, sN, fs, nK, t ): 
100
+    """
101
+    Performs inverse calculation of harmonics contaminating a signal. 
102
+    Args:
103
+        f0 = base frequency of the sinusoidal noise 
104
+        sN = signal containing noise 
105
+        fs = sampling frequency
106
+        nK = number of harmonics to calculate 
107
+        t = time samples 
108
+    """
109
+    print("building Euler matrix ")
110
+    A = np.zeros( (len(t),  nK), dtype=np.complex64)
111
+    for irow, tt in enumerate(t): 
112
+        #A[irow, 0::2] = np.cos( np.arange(nK)*2*np.pi*(f0/fs)*irow )
113
+        #A[irow, 1::2] = np.sin( np.arange(nK)*2*np.pi*(f0/fs)*irow )
114
+        A[irow,:] = np.exp( 1j* np.arange(nK)*2*np.pi*(f0/fs)*irow )
115
+
116
+
117
+    v = np.linalg.lstsq(A, sN, rcond=None) # rcond=None) #, rcond=1e-8)
118
+    #v = sclstsq(A, sN) #, rcond=1e-6)
119
+    alpha = np.real(v[0]) #[0::2]
120
+    beta  = np.imag(v[0]) #[1::2]
121
+    
122
+    #print("Solving A A.T")
123
+    #v = lin.solve(np.dot(A,A.T).T, sN) #, rcond=1e-6)
124
+    #v = np.dot(A.T, v)
125
+    #v = np.dot(np.linalg.inv(np.dot(A.T, A)), np.dot(A.T, sN))
126
+    #alpha = v[0::2]
127
+    #beta  = v[1::2]
128
+
129
+    amp = np.abs(v[0])     #np.sqrt( alpha**2 + beta**2 )
130
+    phase = np.angle(v[0]) # np.arctan(- beta/alpha)
131
+
132
+    #print("amp:", amp, " phase", phase)
133
+
134
+    h = np.zeros(len(t))
135
+    for ik in range(nK):
136
+        h +=  2*amp[ik] * np.cos( 2.*np.pi*ik * (f0/fs) * np.arange(0, len(t), 1 )  + phase[ik] )
137
+
138
+    #plt.matshow(np.imag(A), aspect='auto')
139
+    #plt.colorbar()
140
+
141
+    #plt.figure()
142
+    #plt.plot(alpha)
143
+    #plt.plot(beta)
144
+    #plt.plot(amp)
145
+
146
+    #plt.figure()
147
+    #plt.plot(h)
148
+    #plt.title("modelled noise")
149
+    return sN-h
150
+
90 151
 def jacobian( f0, sN, fs, nK, t):
91 152
     print("building Jacobian matrix ")
92 153
     A = np.zeros( (len(t),  2*nK) )
@@ -102,7 +163,7 @@ def jacobian( f0, sN, fs, nK, t):
102 163
 
103 164
 
104 165
 def harmonicNorm ( f0, sN, fs, nK, t ): 
105
-    return np.linalg.norm( harmonic(f0, sN, fs, nK, t))
166
+    return np.linalg.norm( harmonicEuler(f0, sN, fs, nK, t))
106 167
 
107 168
 def harmonic2Norm ( f0, sN, fs, nK, t ): 
108 169
     return np.linalg.norm(harmonic2(f0[0], f0[1], sN, fs, nK, t))
@@ -112,13 +173,13 @@ def minHarmonic(f0, sN, fs, nK, t):
112 173
     print("minHarmonic", f0, fs, nK, " guess=", f02)
113 174
     res = minimize( harmonicNorm, np.array((f0)), args=(sN, fs, nK, t)) #, method='Nelder-Mead' )# jac=None, hess=None, bounds=None )
114 175
     print(res)
115
-    return harmonic(res.x[0], sN, fs, nK, t)
176
+    return harmonicEuler(res.x[0], sN, fs, nK, t)
116 177
 
117 178
 def minHarmonic2(f1, f2, sN, fs, nK, t):
118 179
     #f02 = guessf0(sN, fs)
119 180
     #print("minHarmonic2", f0, fs, nK, " guess=", f02)
120 181
     #methods with bounds, L-BFGS-B, TNC, SLSQP
121
-    res = minimize( harmonic2Norm, np.array((f1,f2)), args=(sN, fs, nK, t)) #, bounds=((f1-1.,f1+1.0),(f2-1.0,f2+1.0)), method='SLSQP' )
182
+    res = minimize( harmonic2Norm, np.array((f1,f2)), args=(sN, fs, nK, t)) #, bounds=((f1-1.,f1+1.0),(f2-1.0,f2+1.0)), method='TNC' )
122 183
     print(res)
123 184
     return harmonic2(res.x[0], res.x[1], sN, fs, nK, t) 
124 185
 
@@ -138,32 +199,35 @@ if __name__ == "__main__":
138 199
 
139 200
     f0 = 60      # Hz
140 201
     f1 = 60      # Hz
141
-    delta  = 0 #np.random.rand() 
202
+    delta  = np.random.rand() 
142 203
     delta2 = 0 #np.random.rand() 
143 204
     print("delta", delta)
144 205
     fs = 10000   # GMR 
145 206
     t = np.arange(0, 1, 1/fs)
146
-    phi = 0 #np.random.rand() 
207
+    phi = np.random.rand() 
147 208
     phi2 = 0 # np.random.rand() 
148 209
     A =  1.0
149 210
     A2 = 0.0 
150
-    nK = 10
211
+    A3 = 0.0 
212
+    nK = 35
151 213
     T2 = .200
152
-    sN  = A * np.sin( ( 1*(delta  +f0))*2*np.pi*t + phi ) + \
153
-          A2* np.sin( ( 1*(delta2 +f1))*2*np.pi*t + phi2 ) + \
214
+    sN  = A *np.sin( ( 1*(delta  +f0))*2*np.pi*t + phi ) + \
215
+          A2*np.sin( ( 1*(delta2 +f1))*2*np.pi*t + phi2 ) + \
154 216
               np.random.normal(0,.1,len(t)) + \
155
-              + np.exp( -t/T2  ) 
217
+              + A3*np.exp( -t/T2  ) 
156 218
 
157
-    sNc = A * np.sin(  (1*(delta +f0))*2*np.pi*t + phi ) + \
158
-          A2* np.sin(  (1*(delta2+f1))*2*np.pi*t + phi2 ) + \
159
-              + np.exp( -t/T2  ) 
219
+    sNc = A *np.sin(  (1*(delta +f0))*2*np.pi*t + phi ) + \
220
+          A2*np.sin(  (1*(delta2+f1))*2*np.pi*t + phi2 ) + \
221
+              + A3*np.exp( -t/T2  ) 
160 222
 
161 223
 
162 224
     guessf0(sN, fs)
163 225
 
164
-    #h = harmonic( f0, sN, fs, nK, t) 
226
+    #h = harmonicEuler( f0, sN, fs, nK, t) 
165 227
     #h = minHarmonic2( f0, f1, sN, fs, nK, t) 
166
-    h = harmonic2( f0, f1, sN, fs, nK, t) 
228
+    #h = harmonic2( f0, f1, sN, fs, nK, t) 
229
+    
230
+    h = minHarmonic( f0, sN, fs, nK, t) 
167 231
 
168 232
     plt.figure()
169 233
     plt.plot(t, sN, label="sN")

+ 1
- 1
akvo/tressel/mrsurvey.py 查看文件

@@ -558,7 +558,7 @@ class GMRDataProcessor(SNMRDataProcessor):
558 558
                         
559 559
                         self.DATADICT[pulse][ichan][ipm][istack] = harmonic.minHarmonic( f0, self.DATADICT[pulse][ichan][ipm][istack], self.samp, 40, self.DATADICT[pulse]["TIMES"] ) 
560 560
                         #self.DATADICT[pulse][ichan][ipm][istack] = harmonic.minHarmonic2( f0-.25, f0+.25, self.DATADICT[pulse][ichan][ipm][istack], self.samp, 20, self.DATADICT[pulse]["TIMES"] ) 
561
-                        #self.DATADICT[pulse][ichan][ipm][istack] = harmonic.harmonic( f0, self.DATADICT[pulse][ichan][ipm][istack], self.samp, 50, self.DATADICT[pulse]["TIMES"] ) 
561
+                        #self.DATADICT[pulse][ichan][ipm][istack] = harmonic.harmonicEuler( f0, self.DATADICT[pulse][ichan][ipm][istack], self.samp, 20, self.DATADICT[pulse]["TIMES"] ) 
562 562
                
563 563
                         # plot
564 564
                         if plot:

Loading…
取消
儲存