Cleanup of Harmonic modelling, only 1 frequency currently allowed.

akvo/gui/akvoGUI.py

@@ -168,6 +168,11 @@ class ApplicationWindow(QtWidgets.QMainWindow):
         self.ui.calcQGO.pressed.connect( self.calcQ )
         self.ui.FDSmartStackGO.pressed.connect( self.FDSmartStack )
         self.ui.harmonicGO.pressed.connect( self.harmonicModel )       
+
+        self.ui.f0K1Spin.valueChanged.connect( self.LCDHarmonics )
+        self.ui.f0KNSpin.valueChanged.connect( self.LCDHarmonics )
+        self.ui.f0KsSpin.valueChanged.connect( self.LCDHarmonics )
+        self.ui.f0Spin.valueChanged.connect( self.LCDHarmonics )
  
         self.ui.plotQD.setEnabled(False) 
         self.ui.plotQD.pressed.connect( self.plotQD )
@@ -228,7 +233,7 @@ class ApplicationWindow(QtWidgets.QMainWindow):
         self.ui.loopTableWidget.setDragEnabled(False)
         #self.ui.loopTableWidget.setDragDropMode(QtWidgets.QAbstractItemView.InternalMove)
 
-        self.loops = {}        
+        self.loops = {}          
 
         ##########################################################################
         # layer Table 
@@ -254,6 +259,14 @@ class ApplicationWindow(QtWidgets.QMainWindow):
                 self.ui.layerTableWidget.setItem(ir, ic, pCell)
         self.ui.layerTableWidget.cellChanged.connect(self.sigmaCellChanged) 
 
+    def LCDHarmonics(self):
+        self.ui.lcdH1F.setEnabled(True)
+        self.ui.lcdH1F.display( self.ui.f0Spin.value() * self.ui.f0K1Spin.value() )
+        self.ui.lcdHNF.setEnabled(True)
+        self.ui.lcdHNF.display( self.ui.f0Spin.value() * self.ui.f0KNSpin.value() )
+        self.ui.lcdf0NK.setEnabled(True)
+        self.ui.lcdf0NK.display( (self.ui.f0KNSpin.value()+1-self.ui.f0K1Spin.value()) * self.ui.f0KsSpin.value() )
+
     def closeTabs(self):
         #self.ui.ProcTabs.removeTab(idx)    
         self.ui.ProcTabs.clear( )    
@@ -1084,12 +1097,12 @@ class ApplicationWindow(QtWidgets.QMainWindow):
 
     def harmonicModel(self):
 
-
         if "Harmonic modelling" not in self.YamlNode.Processing.keys():
             self.YamlNode.Processing["Harmonic modelling"] = {}
             self.YamlNode.Processing["Harmonic modelling"]["NF"] = str( self.ui.NHarmonicsFreqsSpin.value() ) 
-            self.YamlNode.Processing["Harmonic modelling"]["Nk"] = str( self.ui.NKHarmonicsSpin.value() )
-            self.YamlNode.Processing["Harmonic modelling"]["f0"] = str( self.ui.f0Spin.value() )
+            self.YamlNode.Processing["Harmonic modelling"]["f0K1"] = str( self.ui.f0K1Spin.value() )
+            self.YamlNode.Processing["Harmonic modelling"]["f0KN"] = str( self.ui.f0KNSpin.value() )
+            self.YamlNode.Processing["Harmonic modelling"]["f0Ks"] = str( self.ui.f0KsSpin.value() )
             self.YamlNode.Processing["Harmonic modelling"]["f1"] = str( self.ui.f1Spin.value() )
             self.YamlNode.Processing["Harmonic modelling"]["Segments"] = str( 1 ) # Future 
             self.Log()
@@ -1104,11 +1117,14 @@ class ApplicationWindow(QtWidgets.QMainWindow):
         thread.start_new_thread(self.RAWDataProc.harmonicModel, \
                 ( \
                  self.ui.NHarmonicsFreqsSpin.value(), \
-                 self.ui.NKHarmonicsSpin.value(), \
                  self.ui.f0Spin.value(), \
+                 self.ui.f0K1Spin.value(), \
+                 self.ui.f0KNSpin.value(), \
+                 self.ui.f0KsSpin.value(), \
                  self.ui.f1Spin.value(), \
                  self.ui.plotHarmonic.isChecked(), \
-                 self.ui.mplwidget) \
+                 self.ui.mplwidget \
+                ) \
         )
 
     def FDSmartStack(self):

akvo/gui/main.ui

@@ -1005,7 +1005,7 @@ background: dark grey;
               <property name="maximumSize">
                <size>
                 <width>16777215</width>
-                <height>200</height>
+                <height>350</height>
                </size>
               </property>
               <property name="toolTip">
@@ -1021,26 +1021,121 @@ background: dark grey;
                <bool>true</bool>
               </property>
               <layout class="QGridLayout" name="gridLayout_9">
-               <item row="2" column="0">
-                <widget class="QLabel" name="label_25">
+               <item row="1" column="2">
+                <widget class="Line" name="line_2">
+                 <property name="orientation">
+                  <enum>Qt::Horizontal</enum>
+                 </property>
+                </widget>
+               </item>
+               <item row="8" column="1">
+                <widget class="Line" name="line_8">
+                 <property name="orientation">
+                  <enum>Qt::Horizontal</enum>
+                 </property>
+                </widget>
+               </item>
+               <item row="5" column="0">
+                <widget class="QLabel" name="label_33">
+                 <property name="text">
+                  <string>Sub-harmonics </string>
+                 </property>
+                </widget>
+               </item>
+               <item row="7" column="1">
+                <widget class="QDoubleSpinBox" name="f1Spin">
+                 <property name="maximum">
+                  <double>25000.000000000000000</double>
+                 </property>
+                 <property name="value">
+                  <double>60.000000000000000</double>
+                 </property>
+                </widget>
+               </item>
+               <item row="8" column="2">
+                <widget class="Line" name="line_9">
+                 <property name="orientation">
+                  <enum>Qt::Horizontal</enum>
+                 </property>
+                </widget>
+               </item>
+               <item row="3" column="0">
+                <widget class="QLabel" name="label_11">
                  <property name="text">
-                  <string>Base freq.</string>
+                  <string>First harmonic 1</string>
+                 </property>
+                </widget>
+               </item>
+               <item row="6" column="1">
+                <widget class="Line" name="line_5">
+                 <property name="orientation">
+                  <enum>Qt::Horizontal</enum>
                  </property>
                 </widget>
                </item>
                <item row="1" column="1">
-                <widget class="QSpinBox" name="NKHarmonicsSpin">
-                 <property name="minimum">
-                  <number>5</number>
+                <widget class="Line" name="line">
+                 <property name="orientation">
+                  <enum>Qt::Horizontal</enum>
                  </property>
-                 <property name="maximum">
-                  <number>50</number>
+                </widget>
+               </item>
+               <item row="8" column="0">
+                <widget class="Line" name="line_7">
+                 <property name="orientation">
+                  <enum>Qt::Horizontal</enum>
+                 </property>
+                </widget>
+               </item>
+               <item row="10" column="1">
+                <widget class="QCheckBox" name="plotHarmonic">
+                 <property name="text">
+                  <string>Plot</string>
+                 </property>
+                 <property name="checked">
+                  <bool>true</bool>
+                 </property>
+                </widget>
+               </item>
+               <item row="6" column="2">
+                <widget class="Line" name="line_6">
+                 <property name="orientation">
+                  <enum>Qt::Horizontal</enum>
+                 </property>
+                </widget>
+               </item>
+               <item row="4" column="2">
+                <widget class="QLCDNumber" name="lcdHNF">
+                 <property name="styleSheet">
+                  <string notr="true">#lcdHNF {
+  color: green;
+  background: black;
+}
+
+#lcdHNF:disabled {
+  color: grey;
+  background: dark grey;
+}</string>
+                 </property>
+                </widget>
+               </item>
+               <item row="4" column="1">
+                <widget class="QSpinBox" name="f0KNSpin">
+                 <property name="minimum">
+                  <number>2</number>
                  </property>
                  <property name="value">
                   <number>40</number>
                  </property>
                 </widget>
                </item>
+               <item row="2" column="0">
+                <widget class="QLabel" name="label_25">
+                 <property name="text">
+                  <string>Base freq. 1</string>
+                 </property>
+                </widget>
+               </item>
                <item row="2" column="1">
                 <widget class="QDoubleSpinBox" name="f0Spin">
                  <property name="maximum">
@@ -1054,20 +1149,25 @@ background: dark grey;
                  </property>
                 </widget>
                </item>
-               <item row="3" column="1">
-                <widget class="QDoubleSpinBox" name="f1Spin">
-                 <property name="maximum">
-                  <double>25000.000000000000000</double>
-                 </property>
-                 <property name="value">
-                  <double>60.000000000000000</double>
+               <item row="3" column="2">
+                <widget class="QLCDNumber" name="lcdH1F">
+                 <property name="styleSheet">
+                  <string notr="true">#lcdH1F {
+  color: green;
+  background: black;
+}
+
+#lcdH1F:disabled {
+  color: grey;
+  background: dark grey;
+}</string>
                  </property>
                 </widget>
                </item>
-               <item row="0" column="0">
-                <widget class="QLabel" name="label_26">
+               <item row="7" column="0">
+                <widget class="QLabel" name="label_27">
                  <property name="text">
-                  <string>N freqs</string>
+                  <string>Base freq. 2 </string>
                  </property>
                 </widget>
                </item>
@@ -1077,28 +1177,51 @@ background: dark grey;
                   <number>1</number>
                  </property>
                  <property name="maximum">
-                  <number>2</number>
+                  <number>1</number>
                  </property>
                  <property name="value">
                   <number>1</number>
                  </property>
                 </widget>
                </item>
-               <item row="1" column="0">
-                <widget class="QLabel" name="label_11">
-                 <property name="text">
-                  <string>N harmonics</string>
+               <item row="6" column="0">
+                <widget class="Line" name="line_4">
+                 <property name="orientation">
+                  <enum>Qt::Horizontal</enum>
                  </property>
                 </widget>
                </item>
-               <item row="3" column="0">
-                <widget class="QLabel" name="label_27">
+               <item row="4" column="0">
+                <widget class="QLabel" name="label_31">
                  <property name="text">
-                  <string>Base freq. 2 </string>
+                  <string>Last harmonic 1</string>
                  </property>
                 </widget>
                </item>
-               <item row="3" column="2">
+               <item row="5" column="1">
+                <widget class="QSpinBox" name="f0KsSpin">
+                 <property name="toolTip">
+                  <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;Would you like to calculate subharmonics? For instance, setting this to 1, will calculate only the exact harmonics, setting this to 2 will calculate 1/2 step subharmonics (i.e. if the baseline frequency is 60 Hz, this will result in calcualtion of 30 Hz subharmonics)&lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
+                 </property>
+                 <property name="minimum">
+                  <number>1</number>
+                 </property>
+                 <property name="maximum">
+                  <number>3</number>
+                 </property>
+                 <property name="value">
+                  <number>1</number>
+                 </property>
+                </widget>
+               </item>
+               <item row="1" column="0">
+                <widget class="Line" name="line_3">
+                 <property name="orientation">
+                  <enum>Qt::Horizontal</enum>
+                 </property>
+                </widget>
+               </item>
+               <item row="10" column="2">
                 <widget class="QPushButton" name="harmonicGO">
                  <property name="layoutDirection">
                   <enum>Qt::LeftToRight</enum>
@@ -1112,13 +1235,41 @@ background: dark grey;
                  </property>
                 </widget>
                </item>
-               <item row="2" column="2">
-                <widget class="QCheckBox" name="plotHarmonic">
+               <item row="3" column="1">
+                <widget class="QSpinBox" name="f0K1Spin">
+                 <property name="minimum">
+                  <number>1</number>
+                 </property>
+                 <property name="maximum">
+                  <number>50</number>
+                 </property>
+                 <property name="value">
+                  <number>1</number>
+                 </property>
+                </widget>
+               </item>
+               <item row="0" column="0">
+                <widget class="QLabel" name="label_26">
                  <property name="text">
-                  <string>Plot</string>
+                  <string>N freqs</string>
                  </property>
-                 <property name="checked">
-                  <bool>true</bool>
+                </widget>
+               </item>
+               <item row="5" column="2">
+                <widget class="QLCDNumber" name="lcdf0NK">
+                 <property name="styleSheet">
+                  <string notr="true">#lcdf0NK { 
+ color: green;
+ background: black;
+}
+
+#lcdf0NK:disabled {
+  color: grey;
+  background: dark grey;
+}</string>
+                 </property>
+                 <property name="intValue" stdset="0">
+                  <number>0</number>
                  </property>
                 </widget>
                </item>
@@ -2273,8 +2424,8 @@ background: dark grey;
               <rect>
                <x>0</x>
                <y>0</y>
-               <width>100</width>
-               <height>30</height>
+               <width>96</width>
+               <height>26</height>
               </rect>
              </property>
              <attribute name="label">
@@ -2286,8 +2437,8 @@ background: dark grey;
               <rect>
                <x>0</x>
                <y>0</y>
-               <width>411</width>
-               <height>77</height>
+               <width>96</width>
+               <height>26</height>
               </rect>
              </property>
              <attribute name="label">

akvo/tressel/harmonic.py

@@ -4,18 +4,23 @@ from scipy.optimize import minimize
 from scipy.linalg import lstsq as sclstsq
 import scipy.linalg as lin
 
-def harmonicEuler ( f0, sN, fs, nK, t ): 
+#def harmonicEuler ( f0, sN, fs, nK, t ): 
+def harmonicEuler ( sN, fs, t, f0, k1, kN, ks ): 
     """
     Performs inverse calculation of harmonics contaminating a signal. 
     Args:
-        f0 = base frequency of the sinusoidal noise 
         sN = signal containing noise 
         fs = sampling frequency
-        nK = number of harmonics to calculate 
         t = time samples 
+        f0 = base frequency of the sinusoidal noise 
+        nK = number of harmonics to calculate 
+
     """
     
-    A = np.exp(1j* np.tile( np.arange(1,nK+1),(len(t), 1)) * 2*np.pi* (f0/fs) * np.tile( np.arange(1, len(t)+1, 1),(nK,1)).T  )
+    #A = np.exp(1j* np.tile( np.arange(1,nK+1),(len(t), 1)) * 2*np.pi* (f0/fs) * np.tile( np.arange(1, len(t)+1, 1),(nK,1)).T  )
+    KK = np.arange(k1, kN+1, 1/ks )
+    nK = len(KK)
+    A = np.exp(1j* np.tile(KK,(len(t), 1)) * 2*np.pi* (f0/fs) * np.tile( np.arange(1, len(t)+1, 1),(nK,1)).T  )
 
     v = np.linalg.lstsq(A, sN, rcond=None) 
     alpha = np.real(v[0]) 
@@ -25,45 +30,23 @@ def harmonicEuler ( f0, sN, fs, nK, t ):
     phase = np.angle(v[0]) 
 
     h = np.zeros(len(t))
-    for ik in range(nK):
-        h +=  2*amp[ik] * np.cos( 2.*np.pi*(ik+1) * (f0/fs) * np.arange(1, len(t)+1, 1 )  + phase[ik] )
+    #for ik in range(nK):
+    #    h +=  2*amp[ik] * np.cos( 2.*np.pi*(ik+1) * (f0/fs) * np.arange(1, len(t)+1, 1 )  + phase[ik] )
+    for ik, k in enumerate(KK):
+        h +=  2*amp[ik] * np.cos( 2.*np.pi*(k) * (f0/fs) * np.arange(1, len(t)+1, 1 )  + phase[ik] )
     
     return sN-h
     
     res = sN-h # residual 
 
-    # calculate jacobian 
-    #J = jacEuler(f0, sN, fs, nK, t)
-    #plt.matshow(np.real(J), aspect='auto')
-    
-    #Jv = J * np.tile(v[0], (len(t) ,1)) 
-    #plt.matshow(np.real(Jv), aspect='auto')
-    
-    #print ("shape J",  np.shape(J))
-    #print ("shape v",  np.shape(v[0]))
-    #print ("shape Jv",  np.shape(Jv))
-
-    #plt.figure()
-    #plt.plot(v[0])
-    #plt.figure()
-    #plt.plot(Jv)
-    #plt.show()
-    #exit()
-
-    #jac =  np.linalg.norm( np.dot(Jv.T, res)) 
-    #print("norm ||", jac , "||", " at freq", f0)
+def harmonicNorm (f0, sN, fs, t, k1, kN, ks): 
+    return np.linalg.norm( harmonicEuler(sN, fs, t, f0, k1, kN, ks)) 
 
-    #return res, jac 
-
-
-def jacEuler( f0, sN, fs, nK, t):
-    print("building Jacobian matrix ")
-    J = np.zeros( (len(t),  nK), dtype=np.complex64 )
-    for it, tt in enumerate(t):
-        for ik, k in enumerate( np.arange(0, nK) ):
-            c = 1j*2.*np.pi*(ik+1.)
-            J[it, ik] = (c*np.exp(it*c*f0/fs)) / fs 
-    return J
+def minHarmonic(sN, fs, t, f0, k1, kN, ks):
+    # CG, BFGS, Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov, trust-exact and trust-constr
+    res = minimize(harmonicNorm, np.array((f0)), args=(sN, fs, t, k1, kN, ks), jac='2-point', method='BFGS') # hess=None, bounds=None )
+    print(res)
+    return harmonicEuler(sN, fs, t, res.x[0], k1, kN, ks)#[0]
 
 def harmonicEuler2 ( f0, f1, sN, fs, nK, t ): 
     """
@@ -93,19 +76,18 @@ def harmonicEuler2 ( f0, f1, sN, fs, nK, t ):
 
     return sN-h
 
-def harmonicNorm ( f0, sN, fs, nK, t ): 
-    return np.linalg.norm( harmonicEuler(f0, sN, fs, nK, t) )
-
 def harmonic2Norm ( f0, sN, fs, nK, t ): 
     return np.linalg.norm(harmonicEuler2(f0[0], f0[1], sN, fs, nK, t))
 
-def minHarmonic(f0, sN, fs, nK, t):
-    f02 = guessf0(sN, fs)
-    print("minHarmonic", f0, fs, nK, " guess=", f02)
-    # CG, BFGS, Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov, trust-exact and trust-constr
-    res = minimize(harmonicNorm, np.array((f0)), args=(sN, fs, nK, t), jac='2-point', method='BFGS') #, jac=jacEuler) #, hess=None, bounds=None )
-    print(res)
-    return harmonicEuler(res.x[0], sN, fs, nK, t)#[0]
+#def minHarmonic(f0, sN, fs, nK, t):
+#    f02 = guessf0(sN, fs)
+#    print("minHarmonic", f0, fs, nK, " guess=", f02)
+#    # CG, BFGS, Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov, trust-exact and trust-constr
+#    res = minimize(harmonicNorm, np.array((f0)), args=(sN, fs, nK, t), jac='2-point', method='BFGS') #, jac=jacEuler) #, hess=None, bounds=None )
+#    print(res)
+#    return harmonicEuler(res.x[0], sN, fs, nK, t)#[0]
+
+
 
 def minHarmonic2(f1, f2, sN, fs, nK, t):
     #f02 = guessf0(sN, fs)

akvo/tressel/mrsurvey.py

@@ -517,7 +517,18 @@ class GMRDataProcessor(SNMRDataProcessor):
         canvas.draw()
         self.doneTrigger.emit() 
 
-    def harmonicModel(self, nF, nK, f0, f1, plot, canvas):
+    #def harmonicModel(self, nF, nK, f0, f1, plot, canvas):
+    
+    def harmonicModel(self, nF, f0, f0K1, f0KN, f0Ks, f1, plot, canvas):
+        """ nF = number of base frequencies, must be 1 or 2 
+            f0 = first base frequency  
+            f0K1 = first harmonic to model for first base frequency 
+            f0KN = last harmonic to model for the first base frequency 
+            f0Ks = subharmonic spacing, set to 1 for no subharmonics.
+            plot = should Akvo plot the results 
+            f1 = second base frequency  
+            canvas = mpl plotting axis      
+        """
         #print("harmonic modelling...", f0)
         #plot = True
         if plot:
@@ -543,7 +554,9 @@ class GMRDataProcessor(SNMRDataProcessor):
                                 label = "orig " +  pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " rchan="  + str(ichan))
 
                         if nF == 1:
-                            self.DATADICT[pulse][ichan][ipm][istack] = harmonic.minHarmonic( f0, self.DATADICT[pulse][ichan][ipm][istack], self.samp, nK, self.DATADICT[pulse]["TIMES"] ) 
+                            #self.DATADICT[pulse][ichan][ipm][istack] = harmonic.minHarmonic( f0, self.DATADICT[pulse][ichan][ipm][istack], self.samp, nK, self.DATADICT[pulse]["TIMES"] ) 
+                            self.DATADICT[pulse][ichan][ipm][istack] = harmonic.minHarmonic( self.DATADICT[pulse][ichan][ipm][istack], self.samp,  self.DATADICT[pulse]["TIMES"], \
+                                f0, f0K1, f0KN, f0Ks ) 
                         elif nF == 2:
                             self.DATADICT[pulse][ichan][ipm][istack] = harmonic.minHarmonic2( f0-1e-2, f1+1e-2, self.DATADICT[pulse][ichan][ipm][istack], self.samp, nK, self.DATADICT[pulse]["TIMES"] ) 
 
@@ -559,7 +572,9 @@ class GMRDataProcessor(SNMRDataProcessor):
                                 label = "orig " +  pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " chan="  + str(ichan))
                         
                         if nF == 1:
-                            self.DATADICT[pulse][ichan][ipm][istack] = harmonic.minHarmonic( f0, self.DATADICT[pulse][ichan][ipm][istack], self.samp, nK, self.DATADICT[pulse]["TIMES"] ) 
+                            #self.DATADICT[pulse][ichan][ipm][istack] = harmonic.minHarmonic( f0, self.DATADICT[pulse][ichan][ipm][istack], self.samp, nK, self.DATADICT[pulse]["TIMES"] ) 
+                            self.DATADICT[pulse][ichan][ipm][istack] = harmonic.minHarmonic( self.DATADICT[pulse][ichan][ipm][istack], self.samp,  self.DATADICT[pulse]["TIMES"], \
+                                f0, f0K1, f0KN, f0Ks ) 
                         elif nF == 2:
                             self.DATADICT[pulse][ichan][ipm][istack] = harmonic.minHarmonic2( f0-1e-2, f1+1e-2, self.DATADICT[pulse][ichan][ipm][istack], self.samp, nK, self.DATADICT[pulse]["TIMES"] ) 
                         #self.DATADICT[pulse][ichan][ipm][istack] = harmonic.harmonicEuler( f0, self.DATADICT[pulse][ichan][ipm][istack], self.samp, 20, self.DATADICT[pulse]["TIMES"] ) 

setup.py

@@ -21,7 +21,7 @@ with open("README.md", "r") as fh:
     long_description = fh.read()
 
 setup(name='Akvo',
-      version='1.0.23',
+      version='1.2.1',
       description='Surface nuclear magnetic resonance workbench',
       long_description=long_description,
       long_description_content_type='text/markdown',