Trevor Irons 5 лет назад
Родитель
Сommit
d0b39ca721
4 измененных файлов: 231 добавлений и 23 удалений
  1. 8
    1
      akvo/gui/akvoGUI.py
  2. 31
    2
      akvo/gui/main.ui
  3. 124
    20
      akvo/tressel/harmonic.py
  4. 68
    0
      akvo/tressel/mrsurvey.py

+ 8
- 1
akvo/gui/akvoGUI.py Просмотреть файл

@@ -165,7 +165,8 @@ class ApplicationWindow(QtWidgets.QMainWindow):
165 165
         self.ui.gateIntegrateGO.pressed.connect( self.gateIntegrate )
166 166
         self.ui.calcQGO.pressed.connect( self.calcQ )
167 167
         self.ui.FDSmartStackGO.pressed.connect( self.FDSmartStack )
168
-        
168
+        self.ui.harmonicGO.pressed.connect( self.harmonicModel )       
169
+ 
169 170
         self.ui.plotQD.setEnabled(False) 
170 171
         self.ui.plotQD.pressed.connect( self.plotQD )
171 172
         
@@ -1061,6 +1062,12 @@ class ApplicationWindow(QtWidgets.QMainWindow):
1061 1062
                 (self.ui.CentralVSpinBox.value(), \
1062 1063
                 self.ui.mplwidget))
1063 1064
 
1065
+    def harmonicModel(self):
1066
+        self.lock("harmonic noise modelling")
1067
+        thread.start_new_thread(self.RAWDataProc.harmonicModel, \
1068
+                (self.ui.f0Spin.value(), \
1069
+                self.ui.mplwidget))
1070
+
1064 1071
     def FDSmartStack(self):
1065 1072
 
1066 1073
         if "TD stack" not in self.YamlNode.Processing.keys():

+ 31
- 2
akvo/gui/main.ui Просмотреть файл

@@ -978,6 +978,35 @@ background: dark grey;
978 978
               <property name="checked">
979 979
                <bool>false</bool>
980 980
               </property>
981
+              <widget class="QPushButton" name="harmonicGO">
982
+               <property name="geometry">
983
+                <rect>
984
+                 <x>380</x>
985
+                 <y>120</y>
986
+                 <width>88</width>
987
+                 <height>34</height>
988
+                </rect>
989
+               </property>
990
+               <property name="text">
991
+                <string>PushButton</string>
992
+               </property>
993
+              </widget>
994
+              <widget class="QDoubleSpinBox" name="f0Spin">
995
+               <property name="geometry">
996
+                <rect>
997
+                 <x>160</x>
998
+                 <y>120</y>
999
+                 <width>71</width>
1000
+                 <height>32</height>
1001
+                </rect>
1002
+               </property>
1003
+               <property name="singleStep">
1004
+                <double>0.100000000000000</double>
1005
+               </property>
1006
+               <property name="value">
1007
+                <double>60.000000000000000</double>
1008
+               </property>
1009
+              </widget>
981 1010
              </widget>
982 1011
             </item>
983 1012
             <item>
@@ -2196,8 +2225,8 @@ background: dark grey;
2196 2225
               <rect>
2197 2226
                <x>0</x>
2198 2227
                <y>0</y>
2199
-               <width>100</width>
2200
-               <height>30</height>
2228
+               <width>96</width>
2229
+               <height>26</height>
2201 2230
               </rect>
2202 2231
              </property>
2203 2232
              <attribute name="label">

+ 124
- 20
akvo/tressel/harmonic.py Просмотреть файл

@@ -1,12 +1,55 @@
1 1
 import numpy as np 
2 2
 from scipy.optimize import least_squares 
3
+from scipy.optimize import minimize
4
+from scipy.linalg import lstsq as sclstsq
3 5
 
4
-def harmonic ( sN, f0, fs, nK, t  ): 
6
+def harmonic2 ( f1, f2, sN, fs, nK, t ): 
5 7
     """
6
-    Performs inverse calculation of harmonics contaminating a signal. 
8
+    Performs inverse calculation of two harmonics contaminating a signal. 
7 9
     Args:
10
+        f01 = base frequency of the first sinusoidal noise 
11
+        f02 = base frequency of the second sinusoidal noise 
8 12
         sN = signal containing noise 
13
+        fs = sampling frequency
14
+        nK = number of harmonics to calculate 
15
+        t = time samples 
16
+    """
17
+    print("building matrix ")
18
+    A = np.zeros( (len(t),  4*nK) )
19
+    #f1 = f1MHz * 1e-3
20
+    #f2 = f2MHz * 1e-3
21
+    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 )
25
+        A[irow, 2*nK+1::2] = np.sin( np.arange(nK)*2*np.pi*(f2/fs)*irow )
26
+
27
+    v = np.linalg.lstsq(A, sN, rcond=1e-8)
28
+    #v = sclstsq(A, sN) #, rcond=1e-6)
29
+
30
+    alpha = v[0][0:2*nK:2]
31
+    beta  = v[0][1:2*nK:2]
32
+    amp = np.sqrt( alpha**2 + beta**2 )
33
+    phase = np.arctan(- beta/alpha)
34
+    
35
+    alpha2 = v[0][2*nK::2]
36
+    beta2  = v[0][2*nK+1::2]
37
+    amp2 = np.sqrt( alpha2**2 + beta2**2 )
38
+    phase2 = np.arctan(- beta2/alpha2)
39
+
40
+    h = np.zeros(len(t))
41
+    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
+           + 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
+    return sN-h
46
+
47
+def harmonic ( f0, sN, fs, nK, t ): 
48
+    """
49
+    Performs inverse calculation of harmonics contaminating a signal. 
50
+    Args:
9 51
         f0 = base frequency of the sinusoidal noise 
52
+        sN = signal containing noise 
10 53
         fs = sampling frequency
11 54
         nK = number of harmonics to calculate 
12 55
         t = time samples 
@@ -16,10 +59,6 @@ def harmonic ( sN, f0, fs, nK, t  ):
16 59
     for irow, tt in enumerate(t): 
17 60
         A[irow, 0::2] = np.cos( np.arange(nK)*2*np.pi*(f0/fs)*irow )
18 61
         A[irow, 1::2] = np.sin( np.arange(nK)*2*np.pi*(f0/fs)*irow )
19
-        # brutal 
20
-        #for k, ik in enumerate( np.arange(0, 2*nK, 2) ):
21
-        #    A[irow, ik  ] = np.cos( k*2*np.pi*(f0/fs)*irow )
22
-        #    A[irow, ik+1] = np.sin( k*2*np.pi*(f0/fs)*irow )
23 62
 
24 63
     v = np.linalg.lstsq(A, sN, rcond=None) #, rcond=1e-8)
25 64
 
@@ -29,9 +68,11 @@ def harmonic ( sN, f0, fs, nK, t  ):
29 68
     amp = np.sqrt( alpha**2 + beta**2 )
30 69
     phase = np.arctan(- beta/alpha)
31 70
 
71
+    #print("amp:", amp, " phase", phase)
72
+
32 73
     h = np.zeros(len(t))
33 74
     for ik in range(nK):
34
-        h +=  np.sqrt(alpha[ik]**2 + beta[ik]**2) * np.cos( 2.*np.pi*ik * (f0/fs) * np.arange(0, len(t), 1 )  + phase[ik] )
75
+        h += np.sqrt(alpha[ik]**2 + beta[ik]**2) * np.cos( 2.*np.pi*ik * (f0/fs) * np.arange(0, len(t), 1 )  + phase[ik] )
35 76
 
36 77
     #plt.matshow(A, aspect='auto')
37 78
     #plt.colorbar()
@@ -44,36 +85,99 @@ def harmonic ( sN, f0, fs, nK, t  ):
44 85
     #plt.figure()
45 86
     #plt.plot(h)
46 87
     #plt.title("modelled noise")
88
+    return sN-h
89
+
90
+def jacobian( f0, sN, fs, nK, t):
91
+    print("building Jacobian matrix ")
92
+    A = np.zeros( (len(t),  2*nK) )
93
+    for irow, tt in enumerate(t): 
94
+        #A[irow, 0::2] = np.cos( np.arange(nK)*2*np.pi*(f0/fs)*irow )
95
+        #A[irow, 1::2] = np.sin( np.arange(nK)*2*np.pi*(f0/fs)*irow )
96
+        # brutal 
97
+        for k, ik in enumerate( np.arange(0, 2*nK, 2) ):
98
+            #A[irow, ik  ] = np.cos( k*2*np.pi*(f0/fs)*irow )
99
+            #A[irow, ik+1] = np.sin( k*2*np.pi*(f0/fs)*irow )    
100
+            A[irow, ik  ] = - (2.*np.pi*k*irow * sin((2.*np.pi*irow*f0)/fs)) / fs
101
+            A[irow, ik+1] =   (2.*np.pi*k*irow * cos((2.*np.pi*irow*f0)/fs)) / fs
102
+
103
+
104
+def harmonicNorm ( f0, sN, fs, nK, t ): 
105
+    return np.linalg.norm( harmonic(f0, sN, fs, nK, t))
47 106
 
48
-    return h
107
+def harmonic2Norm ( f0, sN, fs, nK, t ): 
108
+    return np.linalg.norm(harmonic2(f0[0], f0[1], sN, fs, nK, t))
109
+
110
+def minHarmonic(f0, sN, fs, nK, t):
111
+    f02 = guessf0(sN, fs)
112
+    print("minHarmonic", f0, fs, nK, " guess=", f02)
113
+    res = minimize( harmonicNorm, np.array((f0)), args=(sN, fs, nK, t)) #, method='Nelder-Mead' )# jac=None, hess=None, bounds=None )
114
+    print(res)
115
+    return harmonic(res.x[0], sN, fs, nK, t)
116
+
117
+def minHarmonic2(f1, f2, sN, fs, nK, t):
118
+    #f02 = guessf0(sN, fs)
119
+    #print("minHarmonic2", f0, fs, nK, " guess=", f02)
120
+    #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' )
122
+    print(res)
123
+    return harmonic2(res.x[0], res.x[1], sN, fs, nK, t) 
124
+
125
+def guessf0( sN, fs ):
126
+    S = np.fft.fft(sN)
127
+    w = np.fft.fftfreq( len(sN), 1/fs )
128
+    imax = np.argmax( np.abs(S) )
129
+    #plt.plot( w, np.abs(S) )
130
+    #plt.show()
131
+    #print(w)
132
+    #print ( w[imax], w[imax+1] )
133
+    return abs(w[imax])
49 134
 
50 135
 if __name__ == "__main__":
136
+    
51 137
     import matplotlib.pyplot as plt 
52 138
 
53 139
     f0 = 60      # Hz
54
-    delta = np.random.rand()
55
-    fs = 50000  #1e4    
140
+    f1 = 60      # Hz
141
+    delta  = 0 #np.random.rand() 
142
+    delta2 = 0 #np.random.rand() 
143
+    print("delta", delta)
144
+    fs = 10000   # GMR 
56 145
     t = np.arange(0, 1, 1/fs)
57
-    phi = .234
58
-    A = 1.0
59
-    nK = 20
60
-    sN = A * np.sin( (delta+f0)*2*np.pi*t + phi ) + np.random.normal(0,.1,len(t)) 
61
-    sNc = A * np.sin( (delta+f0)*2*np.pi*t + phi ) 
62
-    h = harmonic(sN, f0, fs, nK, t)
146
+    phi = 0 #np.random.rand() 
147
+    phi2 = 0 # np.random.rand() 
148
+    A =  1.0
149
+    A2 = 0.0 
150
+    nK = 10
151
+    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 ) + \
154
+              np.random.normal(0,.1,len(t)) + \
155
+              + np.exp( -t/T2  ) 
156
+
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  ) 
160
+
161
+
162
+    guessf0(sN, fs)
163
+
164
+    #h = harmonic( f0, sN, fs, nK, t) 
165
+    #h = minHarmonic2( f0, f1, sN, fs, nK, t) 
166
+    h = harmonic2( f0, f1, sN, fs, nK, t) 
63 167
 
64 168
     plt.figure()
65 169
     plt.plot(t, sN, label="sN")
66
-    plt.plot(t, sN-h, label="sN-h")
170
+    #plt.plot(t, sN-h, label="sN-h")
67 171
     plt.plot(t, h, label='h')
68
-    plt.title("true noise")
172
+    plt.title("harmonic")
69 173
     plt.legend()
70 174
 
71 175
     plt.figure()
72 176
     plt.plot(t, sN-sNc, label='true noise')
73
-    plt.plot(t, sN-h, label='harmonic removal')
177
+    plt.plot(t, h, label='harmonic removal')
178
+    plt.plot(t, np.exp(-t/T2), label="nmr")
74 179
     plt.legend()
75 180
     plt.title("true noise")
76 181
     
77 182
     plt.show()
78 183
 
79
-    print("hello")

+ 68
- 0
akvo/tressel/mrsurvey.py Просмотреть файл

@@ -25,6 +25,7 @@ import akvo.tressel.decay as decay
25 25
 import akvo.tressel.pca as pca
26 26
 import akvo.tressel.rotate as rotate
27 27
 import akvo.tressel.cmaps as cmaps
28
+import akvo.tressel.harmonic as harmonic
28 29
 
29 30
 import cmocean # colormaps for geophysical data 
30 31
 plt.register_cmap(name='viridis', cmap=cmaps.viridis)
@@ -516,6 +517,73 @@ class GMRDataProcessor(SNMRDataProcessor):
516 517
         canvas.draw()
517 518
         self.doneTrigger.emit() 
518 519
 
520
+    def harmonicModel(self, f0, canvas):
521
+        print("harmonic modelling...", f0)
522
+        plot = True
523
+        if plot:
524
+            canvas.reAx2()
525
+            canvas.ax1.tick_params(axis='both', which='major', labelsize=8)
526
+            canvas.ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y')  
527
+            canvas.ax2.tick_params(axis='both', which='major', labelsize=8)
528
+            canvas.ax2.ticklabel_format(style='sci', scilimits=(0,0), axis='y')  
529
+
530
+        # Data
531
+        iFID = 0
532
+        for pulse in self.DATADICT["PULSES"]:
533
+            self.DATADICT[pulse]["TIMES"] =  self.DATADICT[pulse]["TIMES"]
534
+            for ipm in range(self.DATADICT["nPulseMoments"]):
535
+                for istack in self.DATADICT["stacks"]:
536
+                    canvas.ax1.clear()
537
+                    canvas.ax2.clear()
538
+                    #for ichan in np.append(self.DATADICT[pulse]["chan"], self.DATADICT[pulse]["rchan"]):
539
+                    for ichan in self.DATADICT[pulse]["rchan"]:
540
+                        
541
+                        if plot:
542
+                            canvas.ax1.plot( self.DATADICT[pulse]["TIMES"], 1e9*self.DATADICT[pulse][ichan][ipm][istack], \
543
+                                label = "orig " +  pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " rchan="  + str(ichan))
544
+
545
+                        #self.DATADICT[pulse][ichan][ipm][istack] = harmonic.minHarmonic( f0, self.DATADICT[pulse][ichan][ipm][istack], self.samp, 40, self.DATADICT[pulse]["TIMES"] ) 
546
+                        #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"] ) 
547
+
548
+                        # plot
549
+                        #if plot:
550
+                        #    canvas.ax1.plot( self.DATADICT[pulse]["TIMES"], 1e9*self.DATADICT[pulse][ichan][ipm][istack], \
551
+                        #        label = pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " rchan="  + str(ichan))
552
+
553
+                    for ichan in self.DATADICT[pulse]["chan"]:
554
+                        
555
+                        if plot:
556
+                            canvas.ax2.plot( self.DATADICT[pulse]["TIMES"], 1e9*self.DATADICT[pulse][ichan][ipm][istack], \
557
+                                label = "orig " +  pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " chan="  + str(ichan))
558
+                        
559
+                        self.DATADICT[pulse][ichan][ipm][istack] = harmonic.minHarmonic( f0, self.DATADICT[pulse][ichan][ipm][istack], self.samp, 40, self.DATADICT[pulse]["TIMES"] ) 
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"] ) 
562
+               
563
+                        # plot
564
+                        if plot:
565
+                            canvas.ax2.plot( self.DATADICT[pulse]["TIMES"], 1e9*self.DATADICT[pulse][ichan][ipm][istack], \
566
+                                label = "data " + pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " chan="  + str(ichan))
567
+
568
+                    if plot:
569
+                        canvas.ax1.set_xlabel(r"time [s]", fontsize=8)
570
+                        canvas.ax1.set_ylabel(r"signal [nV]", fontsize=8)
571
+                        canvas.ax2.set_xlabel(r"time [s]", fontsize=8)
572
+                        canvas.ax2.set_ylabel(r"signal [nV]", fontsize=8)
573
+                        canvas.ax1.legend(prop={'size':6})
574
+                        canvas.ax2.legend(prop={'size':6})
575
+                        canvas.draw() 
576
+                    
577
+                percent = (int)(1e2*((float)(iFID*self.DATADICT["nPulseMoments"]+(ipm))/(len(self.DATADICT["PULSES"])*self.nPulseMoments)))
578
+                self.progressTrigger.emit(percent)  
579
+            iFID += 1
580
+        self.doneTrigger.emit() 
581
+        self.updateProcTrigger.emit()  
582
+
583
+
584
+        self.doneTrigger.emit() 
585
+        
586
+    
519 587
     def FDSmartStack(self, outlierTest, MADcutoff, canvas):
520 588
         
521 589
         print("FFT stuff")

Загрузка…
Отмена
Сохранить