Browse Source

cleanup for FastTIMES

tags/1.6.1
Trevor Irons 4 years ago
parent
commit
67a79d8c0f
10 changed files with 602 additions and 668 deletions
  1. 96
    69
      akvo/gui/addCircularLoop.ui
  2. 249
    0
      akvo/gui/addFigure8Loop.ui
  3. 81
    86
      akvo/gui/akvoGUI.py
  4. 27
    0
      akvo/gui/logo.py
  5. 7
    4
      akvo/gui/main.ui
  6. 78
    487
      akvo/tressel/adapt.py
  7. 2
    2
      akvo/tressel/harmonic.py
  8. 53
    17
      akvo/tressel/mrsurvey.py
  9. 4
    0
      pyuic.json
  10. 5
    3
      setup.py

+ 96
- 69
akvo/gui/addCircularLoop.ui View File

7
     <x>0</x>
7
     <x>0</x>
8
     <y>0</y>
8
     <y>0</y>
9
     <width>400</width>
9
     <width>400</width>
10
-    <height>395</height>
10
+    <height>438</height>
11
    </rect>
11
    </rect>
12
   </property>
12
   </property>
13
   <property name="windowTitle">
13
   <property name="windowTitle">
14
    <string>Dialog</string>
14
    <string>Dialog</string>
15
   </property>
15
   </property>
16
   <layout class="QGridLayout" name="gridLayout">
16
   <layout class="QGridLayout" name="gridLayout">
17
+   <item row="2" column="2">
18
+    <widget class="QDoubleSpinBox" name="loopHeight">
19
+     <property name="toolTip">
20
+      <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;Akvo uses a positive down convention, a slight negative value to the loop height improves numerical stability using digital filtering hankel transforms. &lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
21
+     </property>
22
+     <property name="decimals">
23
+      <number>3</number>
24
+     </property>
25
+     <property name="minimum">
26
+      <double>-99.000000000000000</double>
27
+     </property>
28
+     <property name="value">
29
+      <double>-0.001000000000000</double>
30
+     </property>
31
+    </widget>
32
+   </item>
33
+   <item row="1" column="2">
34
+    <widget class="QDoubleSpinBox" name="centreEast">
35
+     <property name="minimum">
36
+      <double>-99999999.000000000000000</double>
37
+     </property>
38
+     <property name="maximum">
39
+      <double>99999999.989999994635582</double>
40
+     </property>
41
+    </widget>
42
+   </item>
17
    <item row="2" column="0">
43
    <item row="2" column="0">
18
     <widget class="QLabel" name="label_3">
44
     <widget class="QLabel" name="label_3">
19
      <property name="text">
45
      <property name="text">
21
      </property>
47
      </property>
22
     </widget>
48
     </widget>
23
    </item>
49
    </item>
24
-   <item row="6" column="2">
50
+   <item row="7" column="2">
25
     <widget class="QDoubleSpinBox" name="dip">
51
     <widget class="QDoubleSpinBox" name="dip">
26
      <property name="enabled">
52
      <property name="enabled">
27
       <bool>false</bool>
53
       <bool>false</bool>
34
      </property>
60
      </property>
35
     </widget>
61
     </widget>
36
    </item>
62
    </item>
37
-   <item row="7" column="0">
38
-    <widget class="QLabel" name="label_8">
39
-     <property name="text">
40
-      <string>azimuth (°)</string>
63
+   <item row="5" column="2">
64
+    <widget class="QSpinBox" name="loopTurns">
65
+     <property name="minimum">
66
+      <number>1</number>
67
+     </property>
68
+     <property name="maximum">
69
+      <number>100000</number>
41
      </property>
70
      </property>
42
     </widget>
71
     </widget>
43
    </item>
72
    </item>
44
-   <item row="0" column="0">
45
-    <widget class="QLabel" name="label">
73
+   <item row="11" column="2">
74
+    <widget class="QDialogButtonBox" name="buttonBox">
75
+     <property name="orientation">
76
+      <enum>Qt::Horizontal</enum>
77
+     </property>
78
+     <property name="standardButtons">
79
+      <set>QDialogButtonBox::Cancel|QDialogButtonBox::Ok</set>
80
+     </property>
81
+    </widget>
82
+   </item>
83
+   <item row="5" column="0">
84
+    <widget class="QLabel" name="label_6">
46
      <property name="text">
85
      <property name="text">
47
-      <string>centre northing (m)</string>
86
+      <string>turns</string>
48
      </property>
87
      </property>
49
     </widget>
88
     </widget>
50
    </item>
89
    </item>
58
      </property>
97
      </property>
59
     </widget>
98
     </widget>
60
    </item>
99
    </item>
61
-   <item row="6" column="0">
62
-    <widget class="QLabel" name="label_7">
63
-     <property name="text">
64
-      <string>dip (°)</string>
100
+   <item row="6" column="2">
101
+    <widget class="QSpinBox" name="segments">
102
+     <property name="toolTip">
103
+      <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;Currently Akvo/Merlin calculates circular loops using segments of wire, forming a polygon. Analytic circular loops may be added in the future. &lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
65
      </property>
104
      </property>
66
-    </widget>
67
-   </item>
68
-   <item row="10" column="2">
69
-    <widget class="QDialogButtonBox" name="buttonBox">
70
-     <property name="orientation">
71
-      <enum>Qt::Horizontal</enum>
105
+     <property name="minimum">
106
+      <number>5</number>
72
      </property>
107
      </property>
73
-     <property name="standardButtons">
74
-      <set>QDialogButtonBox::Cancel|QDialogButtonBox::Ok</set>
108
+     <property name="maximum">
109
+      <number>200</number>
110
+     </property>
111
+     <property name="value">
112
+      <number>15</number>
75
      </property>
113
      </property>
76
     </widget>
114
     </widget>
77
    </item>
115
    </item>
91
      </property>
129
      </property>
92
     </widget>
130
     </widget>
93
    </item>
131
    </item>
94
-   <item row="2" column="2">
95
-    <widget class="QDoubleSpinBox" name="loopHeight">
96
-     <property name="toolTip">
97
-      <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;Akvo uses a positive down convention, a slight negative value to the loop height improves numerical stability using digital filtering hankel transforms. &lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
98
-     </property>
99
-     <property name="decimals">
100
-      <number>3</number>
101
-     </property>
102
-     <property name="minimum">
103
-      <double>-99.000000000000000</double>
132
+   <item row="7" column="0">
133
+    <widget class="QLabel" name="label_7">
134
+     <property name="text">
135
+      <string>dip (°)</string>
104
      </property>
136
      </property>
105
-     <property name="value">
106
-      <double>-0.001000000000000</double>
137
+    </widget>
138
+   </item>
139
+   <item row="3" column="0">
140
+    <widget class="QLabel" name="label_4">
141
+     <property name="text">
142
+      <string>radius (m)</string>
107
      </property>
143
      </property>
108
     </widget>
144
     </widget>
109
    </item>
145
    </item>
110
-   <item row="1" column="0">
111
-    <widget class="QLabel" name="label_2">
146
+   <item row="8" column="0">
147
+    <widget class="QLabel" name="label_8">
112
      <property name="text">
148
      <property name="text">
113
-      <string>centre easting (m)</string>
149
+      <string>azimuth (°)</string>
114
      </property>
150
      </property>
115
     </widget>
151
     </widget>
116
    </item>
152
    </item>
117
-   <item row="7" column="2">
153
+   <item row="8" column="2">
118
     <widget class="QDoubleSpinBox" name="az">
154
     <widget class="QDoubleSpinBox" name="az">
119
      <property name="enabled">
155
      <property name="enabled">
120
       <bool>false</bool>
156
       <bool>false</bool>
121
      </property>
157
      </property>
122
     </widget>
158
     </widget>
123
    </item>
159
    </item>
124
-   <item row="3" column="0">
125
-    <widget class="QLabel" name="label_4">
160
+   <item row="1" column="0">
161
+    <widget class="QLabel" name="label_2">
126
      <property name="text">
162
      <property name="text">
127
-      <string>radius (m)</string>
163
+      <string>centre easting (m)</string>
128
      </property>
164
      </property>
129
     </widget>
165
     </widget>
130
    </item>
166
    </item>
131
-   <item row="4" column="0">
132
-    <widget class="QLabel" name="label_6">
167
+   <item row="0" column="0">
168
+    <widget class="QLabel" name="label">
133
      <property name="text">
169
      <property name="text">
134
-      <string>turns</string>
135
-     </property>
136
-    </widget>
137
-   </item>
138
-   <item row="1" column="2">
139
-    <widget class="QDoubleSpinBox" name="centreEast">
140
-     <property name="minimum">
141
-      <double>-99999999.000000000000000</double>
142
-     </property>
143
-     <property name="maximum">
144
-      <double>99999999.989999994635582</double>
170
+      <string>centre northing (m)</string>
145
      </property>
171
      </property>
146
     </widget>
172
     </widget>
147
    </item>
173
    </item>
148
-   <item row="5" column="0">
174
+   <item row="6" column="0">
149
     <widget class="QLabel" name="label_5">
175
     <widget class="QLabel" name="label_5">
150
      <property name="text">
176
      <property name="text">
151
       <string>segments</string>
177
       <string>segments</string>
152
      </property>
178
      </property>
153
     </widget>
179
     </widget>
154
    </item>
180
    </item>
155
-   <item row="4" column="2">
156
-    <widget class="QSpinBox" name="loopTurns">
157
-     <property name="minimum">
158
-      <number>1</number>
159
-     </property>
160
-     <property name="maximum">
161
-      <number>100000</number>
181
+   <item row="4" column="0">
182
+    <widget class="QLabel" name="label_9">
183
+     <property name="text">
184
+      <string>current flow</string>
162
      </property>
185
      </property>
163
     </widget>
186
     </widget>
164
    </item>
187
    </item>
165
-   <item row="5" column="2">
166
-    <widget class="QSpinBox" name="segments">
167
-     <property name="toolTip">
168
-      <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;Currently Akvo/Merlin calculates circular loops using segments of wire, forming a polygon. Analytic circular loops may be added in the future. &lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
169
-     </property>
170
-     <property name="minimum">
171
-      <number>5</number>
172
-     </property>
188
+   <item row="4" column="2">
189
+    <widget class="QComboBox" name="cwiseBox">
190
+     <item>
191
+      <property name="text">
192
+       <string>clockwise</string>
193
+      </property>
194
+     </item>
195
+     <item>
196
+      <property name="text">
197
+       <string>anti-clockwise</string>
198
+      </property>
199
+     </item>
173
     </widget>
200
     </widget>
174
    </item>
201
    </item>
175
   </layout>
202
   </layout>

+ 249
- 0
akvo/gui/addFigure8Loop.ui View File

1
+<?xml version="1.0" encoding="UTF-8"?>
2
+<ui version="4.0">
3
+ <class>figure8LoopAdd</class>
4
+ <widget class="QDialog" name="figure8LoopAdd">
5
+  <property name="geometry">
6
+   <rect>
7
+    <x>0</x>
8
+    <y>0</y>
9
+    <width>400</width>
10
+    <height>456</height>
11
+   </rect>
12
+  </property>
13
+  <property name="windowTitle">
14
+   <string>Dialog</string>
15
+  </property>
16
+  <layout class="QGridLayout" name="gridLayout">
17
+   <item row="12" column="2">
18
+    <widget class="QDoubleSpinBox" name="doubleSpinBox_3"/>
19
+   </item>
20
+   <item row="7" column="0">
21
+    <widget class="QLabel" name="label_8">
22
+     <property name="text">
23
+      <string>centre easting 2 (m)</string>
24
+     </property>
25
+    </widget>
26
+   </item>
27
+   <item row="14" column="2">
28
+    <widget class="QSpinBox" name="segments">
29
+     <property name="toolTip">
30
+      <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;Currently Akvo/Merlin calculates circular loops using segments of wire, forming a polygon. Analytic circular loops may be added in the future. &lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
31
+     </property>
32
+     <property name="minimum">
33
+      <number>5</number>
34
+     </property>
35
+     <property name="maximum">
36
+      <number>200</number>
37
+     </property>
38
+     <property name="value">
39
+      <number>15</number>
40
+     </property>
41
+    </widget>
42
+   </item>
43
+   <item row="10" column="2">
44
+    <widget class="QDoubleSpinBox" name="loopHeight">
45
+     <property name="toolTip">
46
+      <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;Akvo uses a positive down convention, a slight negative value to the loop height improves numerical stability using digital filtering hankel transforms. &lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
47
+     </property>
48
+     <property name="decimals">
49
+      <number>3</number>
50
+     </property>
51
+     <property name="minimum">
52
+      <double>-99.000000000000000</double>
53
+     </property>
54
+     <property name="value">
55
+      <double>-0.001000000000000</double>
56
+     </property>
57
+    </widget>
58
+   </item>
59
+   <item row="13" column="0">
60
+    <widget class="QLabel" name="label_6">
61
+     <property name="text">
62
+      <string>turns</string>
63
+     </property>
64
+    </widget>
65
+   </item>
66
+   <item row="6" column="2">
67
+    <widget class="QDoubleSpinBox" name="centreNorth2"/>
68
+   </item>
69
+   <item row="8" column="1">
70
+    <widget class="Line" name="line_4">
71
+     <property name="orientation">
72
+      <enum>Qt::Horizontal</enum>
73
+     </property>
74
+    </widget>
75
+   </item>
76
+   <item row="5" column="0">
77
+    <widget class="Line" name="line">
78
+     <property name="orientation">
79
+      <enum>Qt::Horizontal</enum>
80
+     </property>
81
+    </widget>
82
+   </item>
83
+   <item row="1" column="2">
84
+    <widget class="QDoubleSpinBox" name="centreEast1">
85
+     <property name="minimum">
86
+      <double>-99999999.000000000000000</double>
87
+     </property>
88
+     <property name="maximum">
89
+      <double>99999999.989999994635582</double>
90
+     </property>
91
+    </widget>
92
+   </item>
93
+   <item row="12" column="0">
94
+    <widget class="QLabel" name="label_9">
95
+     <property name="text">
96
+      <string>TextLabel</string>
97
+     </property>
98
+    </widget>
99
+   </item>
100
+   <item row="1" column="0">
101
+    <widget class="QLabel" name="label_2">
102
+     <property name="text">
103
+      <string>centre easting 1 (m)</string>
104
+     </property>
105
+    </widget>
106
+   </item>
107
+   <item row="8" column="2">
108
+    <widget class="Line" name="line_3">
109
+     <property name="orientation">
110
+      <enum>Qt::Horizontal</enum>
111
+     </property>
112
+    </widget>
113
+   </item>
114
+   <item row="0" column="0">
115
+    <widget class="QLabel" name="label">
116
+     <property name="text">
117
+      <string>centre northing 1 (m)</string>
118
+     </property>
119
+    </widget>
120
+   </item>
121
+   <item row="6" column="0">
122
+    <widget class="QLabel" name="label_7">
123
+     <property name="text">
124
+      <string>centre northing 2 (m)</string>
125
+     </property>
126
+    </widget>
127
+   </item>
128
+   <item row="14" column="0">
129
+    <widget class="QLabel" name="label_5">
130
+     <property name="text">
131
+      <string>segments</string>
132
+     </property>
133
+    </widget>
134
+   </item>
135
+   <item row="17" column="2">
136
+    <widget class="QDialogButtonBox" name="buttonBox">
137
+     <property name="orientation">
138
+      <enum>Qt::Horizontal</enum>
139
+     </property>
140
+     <property name="standardButtons">
141
+      <set>QDialogButtonBox::Cancel|QDialogButtonBox::Ok</set>
142
+     </property>
143
+    </widget>
144
+   </item>
145
+   <item row="5" column="2">
146
+    <widget class="Line" name="line_2">
147
+     <property name="orientation">
148
+      <enum>Qt::Horizontal</enum>
149
+     </property>
150
+    </widget>
151
+   </item>
152
+   <item row="7" column="2">
153
+    <widget class="QDoubleSpinBox" name="centreEast2"/>
154
+   </item>
155
+   <item row="10" column="0">
156
+    <widget class="QLabel" name="label_3">
157
+     <property name="text">
158
+      <string>height (m)</string>
159
+     </property>
160
+    </widget>
161
+   </item>
162
+   <item row="0" column="2">
163
+    <widget class="QDoubleSpinBox" name="centreNorth1">
164
+     <property name="minimum">
165
+      <double>-999999999.000000000000000</double>
166
+     </property>
167
+     <property name="maximum">
168
+      <double>9999999999.000000000000000</double>
169
+     </property>
170
+    </widget>
171
+   </item>
172
+   <item row="13" column="2">
173
+    <widget class="QSpinBox" name="loopTurns">
174
+     <property name="minimum">
175
+      <number>1</number>
176
+     </property>
177
+     <property name="maximum">
178
+      <number>100000</number>
179
+     </property>
180
+    </widget>
181
+   </item>
182
+   <item row="8" column="0">
183
+    <widget class="Line" name="line_5">
184
+     <property name="orientation">
185
+      <enum>Qt::Horizontal</enum>
186
+     </property>
187
+    </widget>
188
+   </item>
189
+   <item row="9" column="0">
190
+    <widget class="QLabel" name="label_4">
191
+     <property name="text">
192
+      <string>radius  (m)</string>
193
+     </property>
194
+    </widget>
195
+   </item>
196
+   <item row="9" column="2">
197
+    <widget class="QDoubleSpinBox" name="loopRadius">
198
+     <property name="toolTip">
199
+      <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;Radius of the loop&lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
200
+     </property>
201
+     <property name="minimum">
202
+      <double>0.100000000000000</double>
203
+     </property>
204
+     <property name="maximum">
205
+      <double>600.000000000000000</double>
206
+     </property>
207
+     <property name="value">
208
+      <double>50.000000000000000</double>
209
+     </property>
210
+    </widget>
211
+   </item>
212
+  </layout>
213
+ </widget>
214
+ <resources/>
215
+ <connections>
216
+  <connection>
217
+   <sender>buttonBox</sender>
218
+   <signal>accepted()</signal>
219
+   <receiver>figure8LoopAdd</receiver>
220
+   <slot>accept()</slot>
221
+   <hints>
222
+    <hint type="sourcelabel">
223
+     <x>248</x>
224
+     <y>254</y>
225
+    </hint>
226
+    <hint type="destinationlabel">
227
+     <x>157</x>
228
+     <y>274</y>
229
+    </hint>
230
+   </hints>
231
+  </connection>
232
+  <connection>
233
+   <sender>buttonBox</sender>
234
+   <signal>rejected()</signal>
235
+   <receiver>figure8LoopAdd</receiver>
236
+   <slot>reject()</slot>
237
+   <hints>
238
+    <hint type="sourcelabel">
239
+     <x>316</x>
240
+     <y>260</y>
241
+    </hint>
242
+    <hint type="destinationlabel">
243
+     <x>286</x>
244
+     <y>274</y>
245
+    </hint>
246
+   </hints>
247
+  </connection>
248
+ </connections>
249
+</ui>

+ 81
- 86
akvo/gui/akvoGUI.py View File

1
 #/usr/bin/env python 
1
 #/usr/bin/env python 
2
 import sys
2
 import sys
3
-#import readline
4
-try: 
5
-    from akvo.gui.main_ui import Ui_MainWindow
6
-    uicerr = False
7
-except: # Fallback 
8
-    from akvo.gui.mainui import Ui_MainWindow
9
-    uicerr = """
10
-USING THE DEFAULT GUI FILES, AKVO MAY NOT WORK CORRECTLY!
11
-
12
-See INSTALL.txt for details regarding GUI configuration 
13
-if you are encountering problems.     
14
-
15
-Clicking ignore will prevent this warning from showing 
16
-each time you launch Akvo.                  
17
-"""
18
-
19
-from akvo.gui.addCircularLoop_ui import Ui_circularLoopAdd
20
 
3
 
21
 import matplotlib
4
 import matplotlib
22
 matplotlib.use("QT5Agg")
5
 matplotlib.use("QT5Agg")
23
-from PyQt5 import QtCore, QtGui, QtWidgets
24
-
6
+from PyQt5 import QtCore, QtGui, QtWidgets #, uic
25
 import numpy as np
7
 import numpy as np
26
 import time
8
 import time
27
 import os
9
 import os
28
 from copy import deepcopy
10
 from copy import deepcopy
29
 from matplotlib.backends.backend_qt4 import NavigationToolbar2QT #as NavigationToolbar
11
 from matplotlib.backends.backend_qt4 import NavigationToolbar2QT #as NavigationToolbar
30
 import datetime, time
12
 import datetime, time
13
+import pkg_resources  # part of setuptools
14
+from collections import OrderedDict
15
+from ruamel import yaml
31
 
16
 
17
+from akvo.gui.main_ui import Ui_MainWindow
18
+from akvo.gui.addCircularLoop_ui import Ui_circularLoopAdd
19
+from akvo.gui.addFigure8Loop_ui import Ui_figure8LoopAdd
32
 from akvo.tressel import mrsurvey 
20
 from akvo.tressel import mrsurvey 
33
 
21
 
34
 from pyLemma import LemmaCore  
22
 from pyLemma import LemmaCore  
35
 from pyLemma import FDEM1D 
23
 from pyLemma import FDEM1D 
36
 from pyLemma import Merlin
24
 from pyLemma import Merlin
37
 
25
 
38
-import pkg_resources  # part of setuptools
39
-version = pkg_resources.require("Akvo")[0].version
26
+VERSION = pkg_resources.require("Akvo")[0].version
40
 
27
 
41
-from collections import OrderedDict
42
-
43
-from ruamel import yaml
44
-#import ruamel.yaml 
45
-#yaml = ruamel.yaml.YAML()
46
-#yaml.indent(mapping=4)
47
-
48
-#import yaml
49
 # Writes out numpy arrays into Eigen vectors as serialized by Lemma
28
 # Writes out numpy arrays into Eigen vectors as serialized by Lemma
50
 class MatrixXr(yaml.YAMLObject):
29
 class MatrixXr(yaml.YAMLObject):
51
     yaml_tag = u'MatrixXr'
30
     yaml_tag = u'MatrixXr'
75
 class AkvoYamlNode(yaml.YAMLObject):
54
 class AkvoYamlNode(yaml.YAMLObject):
76
     yaml_tag = u'AkvoData'
55
     yaml_tag = u'AkvoData'
77
     def __init__(self):
56
     def __init__(self):
78
-        self.Akvo_VERSION = version
57
+        self.Akvo_VERSION = VERSION
79
         self.Import = OrderedDict() # {}
58
         self.Import = OrderedDict() # {}
80
         self.Processing = [] # OrderedDict() 
59
         self.Processing = [] # OrderedDict() 
81
         self.Stacking = OrderedDict()
60
         self.Stacking = OrderedDict()
82
         self.META = OrderedDict() 
61
         self.META = OrderedDict() 
83
-    #def __init__(self, node):
84
-    #    self.Akvo_VERSION = node["version"]
85
-    #    self.Import = OrderedDict( node["Import"] ) # {}
86
-    #    self.Processing = OrderedDict( node["Processing"] ) 
87
     def __repr__(self):
62
     def __repr__(self):
88
         return "%s(name=%r, Akvo_VERSION=%r, Import=%r, Processing=%r, self.Stacking=%r, self.META=%r)" % (
63
         return "%s(name=%r, Akvo_VERSION=%r, Import=%r, Processing=%r, self.Stacking=%r, self.META=%r)" % (
89
             self.__class__.__name__, self.Akvo_VERSION, self.Import, self.Processing, self.Stacking, self.META )
64
             self.__class__.__name__, self.Akvo_VERSION, self.Import, self.Processing, self.Stacking, self.META )
90
-            #self.__class__.__name__, self.Akvo_VERSION, self.Import, OrderedDict(self.Processing) ) 
91
     
65
     
92
 try:    
66
 try:    
93
     import thread 
67
     import thread 
107
     
81
     
108
     def __init__(self):
82
     def __init__(self):
109
         
83
         
110
-        QtWidgets.QMainWindow.__init__(self)
111
-        self.setAttribute(QtCore.Qt.WA_DeleteOnClose)       
112
- 
113
-        akvohome = os.path.expanduser("~") + "/.akvo"
114
-        if not os.path.exists(akvohome):
115
-            os.makedirs(akvohome)
84
+        super().__init__() 
85
+        #QtWidgets.QMainWindow.__init__(self)
86
+        self.setAttribute(QtCore.Qt.WA_DeleteOnClose)      
87
+
88
+        # alternative to calling pyuic 
89
+        #self.ui = uic.loadUi('main.ui', self) 
116
 
90
 
117
         self.ui = Ui_MainWindow()
91
         self.ui = Ui_MainWindow()
118
         self.ui.setupUi(self)
92
         self.ui.setupUi(self)
119
-
120
-        if uicerr != False and not os.path.exists(akvohome+"/pyuic-warned"):
121
-            reply = QtGui.QMessageBox.warning(self, 'Warning', uicerr, QtGui.QMessageBox.Ok, QtGui.QMessageBox.Ignore) 
122
-            if reply == 1024:      # "0x400" in hex
123
-                pass
124
-            elif reply == 1048576: # "0x100000" in hex
125
-                warn = open( akvohome+"/pyuic-warned" ,"w" )
126
-                warn.write("Gui files were not compiled locally using pyuic! Further warnings have been supressed")
127
-                warn.close()
93
+    
128
         self.RAWDataProc = None 
94
         self.RAWDataProc = None 
129
-
130
         self.YamlNode = AkvoYamlNode()           
95
         self.YamlNode = AkvoYamlNode()           
131
  
96
  
132
         # initialise some stuff
97
         # initialise some stuff
159
         ###########
124
         ###########
160
         # Buttons #
125
         # Buttons #
161
         ###########
126
         ###########
162
-#         #QtCore.QObject.connect(self.ui.fullWorkflowPushButton, QtCore.SIGNAL("clicked()"), self.preprocess )
163
         self.ui.loadDataPushButton.pressed.connect(self.loadRAW) 
127
         self.ui.loadDataPushButton.pressed.connect(self.loadRAW) 
164
         self.ui.sumDataGO.pressed.connect( self.sumDataChans )
128
         self.ui.sumDataGO.pressed.connect( self.sumDataChans )
165
         self.ui.bandPassGO.pressed.connect( self.bandPassFilter )
129
         self.ui.bandPassGO.pressed.connect( self.bandPassFilter )
167
         self.ui.fdDesignPushButton.pressed.connect( self.designFDFilter )
131
         self.ui.fdDesignPushButton.pressed.connect( self.designFDFilter )
168
         self.ui.downSampleGO.pressed.connect( self.downsample )
132
         self.ui.downSampleGO.pressed.connect( self.downsample )
169
         self.ui.windowFilterGO.pressed.connect( self.windowFilter )
133
         self.ui.windowFilterGO.pressed.connect( self.windowFilter )
170
-#        self.ui.despikeGO.pressed.connect( self.despikeFilter ) # use smart stack instead 
171
         self.ui.adaptGO.pressed.connect( self.adaptFilter )
134
         self.ui.adaptGO.pressed.connect( self.adaptFilter )
172
         self.ui.adaptFDGO.pressed.connect( self.adaptFilterFD )
135
         self.ui.adaptFDGO.pressed.connect( self.adaptFilterFD )
173
         self.ui.qdGO.pressed.connect( self.quadDet )
136
         self.ui.qdGO.pressed.connect( self.quadDet )
207
         self.ui.dateEdit.dateChanged.connect( self.logSite )
170
         self.ui.dateEdit.dateChanged.connect( self.logSite )
208
         # this may call the yaml stuff too often... 
171
         # this may call the yaml stuff too often... 
209
         self.ui.txtComments.textChanged.connect( self.logSite )
172
         self.ui.txtComments.textChanged.connect( self.logSite )
210
-        
211
         self.ui.plotLoops.pressed.connect( self.plotLoops2 )       
173
         self.ui.plotLoops.pressed.connect( self.plotLoops2 )       
174
+        self.ui.removeLoopButton.pressed.connect( self.removeLoop )       
212
  
175
  
213
         # Loops
176
         # Loops
214
         self.ui.addLoopButton.pressed.connect( self.loopAdd )
177
         self.ui.addLoopButton.pressed.connect( self.loopAdd )
373
             Error.setDetailedText("Each loop label must be unique and comprise at least one character. Leading and trailing whitespace will be trimmed.")
336
             Error.setDetailedText("Each loop label must be unique and comprise at least one character. Leading and trailing whitespace will be trimmed.")
374
             Error.exec_()
337
             Error.exec_()
375
         else:
338
         else:
339
+            
340
+            ### Circular loop
376
             if self.ui.loopGeom.currentText() == "Circular":
341
             if self.ui.loopGeom.currentText() == "Circular":
377
                 dialog = QtWidgets.QDialog()
342
                 dialog = QtWidgets.QDialog()
378
                 dialog.ui = Ui_circularLoopAdd()
343
                 dialog.ui = Ui_circularLoopAdd()
379
                 dialog.ui.setupUi(dialog)
344
                 dialog.ui.setupUi(dialog)
380
                 dialog.exec_()
345
                 dialog.exec_()
381
                 dialog.show()
346
                 dialog.show()
347
+
382
                 if dialog.result():
348
                 if dialog.result():
383
                     cn = dialog.ui.centreNorth.value()
349
                     cn = dialog.ui.centreNorth.value()
384
                     ce = dialog.ui.centreEast.value()
350
                     ce = dialog.ui.centreEast.value()
386
                     rad = dialog.ui.loopRadius.value()
352
                     rad = dialog.ui.loopRadius.value()
387
                     turns = dialog.ui.loopTurns.value()
353
                     turns = dialog.ui.loopTurns.value()
388
                     ns = dialog.ui.segments.value()
354
                     ns = dialog.ui.segments.value()
389
-                    #print("dip", dialog.ui.dip.value())
390
-                    #print("azimuth", dialog.ui.az.value())
355
+                    cwise = dialog.ui.cwiseBox.currentIndex()
356
+                    print("cwise", cwise)
357
+                    #dip = dialog.ui.dip.value()
358
+                    #azimuth = dialog.ui.az.value()
391
                     
359
                     
392
                     self.loops[self.ui.loopLabel.text()] = FDEM1D.PolygonalWireAntenna()
360
                     self.loops[self.ui.loopLabel.text()] = FDEM1D.PolygonalWireAntenna()
393
                     self.loops[self.ui.loopLabel.text()].SetNumberOfPoints( dialog.ui.segments.value() + 1 )
361
                     self.loops[self.ui.loopLabel.text()].SetNumberOfPoints( dialog.ui.segments.value() + 1 )
395
 
363
 
396
                     points = np.linspace(0, 2*np.pi, dialog.ui.segments.value()+1)
364
                     points = np.linspace(0, 2*np.pi, dialog.ui.segments.value()+1)
397
                     for iseg, ipt in enumerate(points):
365
                     for iseg, ipt in enumerate(points):
398
-                        #self.loops[self.ui.loopLabel.text()].SetPoint(iseg,( dialog.ui.centreNorth.value(),  dialog.ui.centreEast.value(), dialog.ui.loopHeight.value()))
399
-                        self.loops[self.ui.loopLabel.text()].SetPoint(iseg, ( cn+rad*np.sin(ipt), ce+rad*np.cos(ipt), ht ))
400
-                        #self.loops[self.ui.loopLabel.text()].SetPoint(1,(150,  50, -1e-3))
401
-                    
366
+                        if cwise == 0:
367
+                            self.loops[self.ui.loopLabel.text()].SetPoint(iseg, (  cn+rad*np.sin(ipt), ce+rad*np.cos(ipt), ht) )
368
+                        else:
369
+                            self.loops[self.ui.loopLabel.text()].SetPoint(iseg, ( -cn+rad*np.sin(ipt), ce+rad*np.cos(ipt), ht) )
402
                     self.loops[self.ui.loopLabel.text()].SetNumberOfFrequencies(1)
370
                     self.loops[self.ui.loopLabel.text()].SetNumberOfFrequencies(1)
403
-                    #self.loops[self.ui.loopLabel.text].SetFrequency(0,2246) 
404
                     self.loops[self.ui.loopLabel.text()].SetCurrent(1.)
371
                     self.loops[self.ui.loopLabel.text()].SetCurrent(1.)
405
-                   
406
-                    # update the table 
407
-                    self.ui.txRxTable.setRowCount( len(self.loops.keys()) )
372
+            
373
+            if self.ui.loopGeom.currentText() == "figure-8":
374
+                dialog = QtWidgets.QDialog()
375
+                dialog.ui = Ui_figure8LoopAdd()
376
+                dialog.ui.setupUi(dialog)
377
+                dialog.exec_()
378
+                dialog.show()
379
+                
380
+                if dialog.result():
381
+                    cn1 = dialog.ui.centreNorth1.value()
382
+                    ce1 = dialog.ui.centreEast1.value()
383
+                    cn2 = dialog.ui.centreNorth2.value()
384
+                    ce2 = dialog.ui.centreEast2.value()
385
+                    ht = dialog.ui.loopHeight.value()
386
+                    rad = dialog.ui.loopRadius.value()
387
+                    turns = dialog.ui.loopTurns.value()
388
+                    ns = dialog.ui.segments.value()
389
+                    #cwise = dialog.ui.cwiseBox.currentIndex()
390
+                    print(cn1, ce1, cn2, ce2, ht, rad, turns, ns)   
391
+ 
392
+            # general across all types   
393
+            if dialog.result():
394
+                # update the table 
395
+                self.ui.txRxTable.setRowCount( len(self.loops.keys()) )
408
                     
396
                     
409
-                    pCell = QtWidgets.QTableWidgetItem()
410
-                    pCell.setText( self.ui.loopLabel.text() ) 
411
-                    pCell.setFlags( QtCore.Qt.ItemIsSelectable | QtCore.Qt.ItemIsEnabled )
412
-                    self.ui.txRxTable.setItem( len(self.loops.keys())-1, 0, pCell)
397
+                pCell = QtWidgets.QTableWidgetItem()
398
+                pCell.setText( self.ui.loopLabel.text() ) 
399
+                pCell.setFlags( QtCore.Qt.ItemIsSelectable | QtCore.Qt.ItemIsEnabled )
400
+                self.ui.txRxTable.setItem( len(self.loops.keys())-1, 0, pCell)
413
                     
401
                     
414
-                    gCell = QtWidgets.QTableWidgetItem()
415
-                    gCell.setText( self.ui.loopGeom.currentText() ) 
416
-                    gCell.setFlags( QtCore.Qt.ItemIsEnabled )
417
-                    self.ui.txRxTable.setItem( len(self.loops.keys())-1, 1, gCell)
402
+                gCell = QtWidgets.QTableWidgetItem()
403
+                gCell.setText( self.ui.loopGeom.currentText() ) 
404
+                gCell.setFlags( QtCore.Qt.ItemIsEnabled )
405
+                self.ui.txRxTable.setItem( len(self.loops.keys())-1, 1, gCell)
418
                     
406
                     
419
-                    tCell = QtWidgets.QTableWidgetItem()
420
-                    tCell.setText( str(dialog.ui.loopTurns.value()) ) 
421
-                    tCell.setFlags( QtCore.Qt.ItemIsEnabled )
422
-                    self.ui.txRxTable.setItem( len(self.loops.keys())-1, 2, tCell)
407
+                tCell = QtWidgets.QTableWidgetItem()
408
+                tCell.setText( str(dialog.ui.loopTurns.value()) ) 
409
+                tCell.setFlags( QtCore.Qt.ItemIsEnabled )
410
+                self.ui.txRxTable.setItem( len(self.loops.keys())-1, 2, tCell)
423
                     
411
                     
424
-                    txCell = QtWidgets.QTableWidgetItem()
425
-                    txCell.setText( str(self.ui.loopType.currentText()) ) 
426
-                    txCell.setFlags( QtCore.Qt.ItemIsEnabled )
427
-                    self.ui.txRxTable.setItem( len(self.loops.keys())-1, 3, txCell)
428
-
412
+                txCell = QtWidgets.QTableWidgetItem()
413
+                txCell.setText( str(self.ui.loopType.currentText()) ) 
414
+                txCell.setFlags( QtCore.Qt.ItemIsEnabled )
415
+                self.ui.txRxTable.setItem( len(self.loops.keys())-1, 3, txCell)
429
 
416
 
430
     def headerBoxShrink(self):
417
     def headerBoxShrink(self):
431
         #self.ui.headerFileBox.setVisible(False)
418
         #self.ui.headerFileBox.setVisible(False)
536
         #self.ui.loopTableWidget.itemClicked.connect(self.loopCellClicked) 
523
         #self.ui.loopTableWidget.itemClicked.connect(self.loopCellClicked) 
537
 
524
 
538
     def loopCellChanged(self):
525
     def loopCellChanged(self):
539
-        self.ui.loopTableWidget.cellChanged.disconnect(self.loopCellChanged) 
540
         
526
         
527
+        self.ui.loopTableWidget.cellChanged.disconnect(self.loopCellChanged) 
541
         jj = self.ui.loopTableWidget.currentColumn()
528
         jj = self.ui.loopTableWidget.currentColumn()
542
         ii = self.ui.loopTableWidget.currentRow()
529
         ii = self.ui.loopTableWidget.currentRow()
543
 
530
 
581
 
568
 
582
         for loop in self.loops:
569
         for loop in self.loops:
583
             POINTS = self.loops[loop].GetPoints().T
570
             POINTS = self.loops[loop].GetPoints().T
584
-            self.ui.mplwidget.ax1.plot( POINTS[:,0], POINTS[:,1] )
585
-        
571
+            self.ui.mplwidget.ax1.plot( POINTS[:,1], POINTS[:,0], label=loop )
572
+       
573
+        self.ui.mplwidget.ax1.spines['right'].set_visible(False)
574
+        self.ui.mplwidget.ax1.spines['top'].set_visible(False) 
575
+        self.ui.mplwidget.ax1.set_xlabel("easting (m)")
576
+        self.ui.mplwidget.ax1.set_ylabel("northing (m)")
577
+        self.ui.mplwidget.ax1.legend()
586
         self.ui.mplwidget.ax1.set_aspect('equal') #, adjustable='box')
578
         self.ui.mplwidget.ax1.set_aspect('equal') #, adjustable='box')
587
         self.ui.mplwidget.draw()
579
         self.ui.mplwidget.draw()
580
+    
581
+    def removeLoop(self):
582
+        del self.loops[ self.ui.txRxTable.item( self.ui.txRxTable.currentRow(), 0).text() ]
583
+        self.ui.txRxTable.removeRow(self.ui.txRxTable.currentRow()) 
588
 
584
 
589
     def plotLoops(self):
585
     def plotLoops(self):
590
        
586
        
1639
         splash.show()
1635
         splash.show()
1640
     
1636
     
1641
     aw = ApplicationWindow()
1637
     aw = ApplicationWindow()
1642
-
1643
     #img=mpimg.imread(logo)
1638
     #img=mpimg.imread(logo)
1644
     for ax in [ aw.ui.mplwidget ]: 
1639
     for ax in [ aw.ui.mplwidget ]: 
1645
         ax.fig.clear()
1640
         ax.fig.clear()
1672
         splash.finish(aw)
1667
         splash.finish(aw)
1673
         #time.sleep(1) 
1668
         #time.sleep(1) 
1674
 
1669
 
1675
-    aw.setWindowTitle("Akvo v"+str(version)) 
1670
+    aw.setWindowTitle("Akvo v"+str(VERSION)) 
1676
     aw.show()
1671
     aw.show()
1677
     qApp.setWindowIcon(QtGui.QIcon(logo2))
1672
     qApp.setWindowIcon(QtGui.QIcon(logo2))
1678
     sys.exit(qApp.exec_())
1673
     sys.exit(qApp.exec_())

+ 27
- 0
akvo/gui/logo.py View File

38
     ax.plot(x[o],  y4[o], linewidth=3, color='black')
38
     ax.plot(x[o],  y4[o], linewidth=3, color='black')
39
     ax.plot(x[o], -y4[o], linewidth=3, color='black')
39
     ax.plot(x[o], -y4[o], linewidth=3, color='black')
40
     ax.fill_between(x[o],  y4[o], y2=-y4[o],  where=y4[o]<=1, interpolate=True, linewidth=0, alpha=.95, color='maroon')
40
     ax.fill_between(x[o],  y4[o], y2=-y4[o],  where=y4[o]<=1, interpolate=True, linewidth=0, alpha=.95, color='maroon')
41
+
42
+if __name__ == "__main__":
43
+    import matplotlib.pyplot as plt 
44
+
45
+    fig = plt.figure( figsize=(6,3) )
46
+    #ax = fig.add_axes([.1,.1,.8,.8])
47
+    ax = fig.add_subplot(211)
48
+
49
+    fig.patch.set_facecolor( None )
50
+    fig.patch.set_alpha( .0 )
51
+    ax.axis('off')
52
+    plotLogo(ax)
53
+    ax.xaxis.set_major_locator(plt.NullLocator()) 
54
+    ax.yaxis.set_major_locator(plt.NullLocator()) 
55
+
56
+    subplot2 = fig.add_subplot(212)
57
+    subplot2.text(0.5, 1.,'surface NMR workbench',
58
+            horizontalalignment='center',
59
+            verticalalignment='center',
60
+            size=22,
61
+            transform = subplot2.transAxes)
62
+    subplot2.xaxis.set_major_locator(plt.NullLocator()) 
63
+    subplot2.yaxis.set_major_locator(plt.NullLocator()) 
64
+    subplot2.axis('off')
65
+    plt.savefig("logo.pdf")
66
+    plt.show()
67
+
41
 #ax.fill_between(x[o], -y4[o], y2=0, where=-y4[o]<=1, interpolate=True, linewidth=0, alpha=.5, color='black')
68
 #ax.fill_between(x[o], -y4[o], y2=0, where=-y4[o]<=1, interpolate=True, linewidth=0, alpha=.5, color='black')
42
 #ax.plot(x[o], y2[o], linewidth=3, color='black')
69
 #ax.plot(x[o], y2[o], linewidth=3, color='black')
43
 #ax.plot(x[o],-y2[o], linewidth=3, color='black')
70
 #ax.plot(x[o],-y2[o], linewidth=3, color='black')

+ 7
- 4
akvo/gui/main.ui View File

760
            <enum>Qt::LeftToRight</enum>
760
            <enum>Qt::LeftToRight</enum>
761
           </property>
761
           </property>
762
           <property name="currentIndex">
762
           <property name="currentIndex">
763
-           <number>3</number>
763
+           <number>0</number>
764
           </property>
764
           </property>
765
           <property name="elideMode">
765
           <property name="elideMode">
766
            <enum>Qt::ElideLeft</enum>
766
            <enum>Qt::ElideLeft</enum>
2253
                  <property name="toolTip">
2253
                  <property name="toolTip">
2254
                   <string>Forgetting factor, how quickly does the filter adapt.</string>
2254
                   <string>Forgetting factor, how quickly does the filter adapt.</string>
2255
                  </property>
2255
                  </property>
2256
+                 <property name="decimals">
2257
+                  <number>3</number>
2258
+                 </property>
2256
                  <property name="minimum">
2259
                  <property name="minimum">
2257
-                  <double>0.200000000000000</double>
2260
+                  <double>0.950000000000000</double>
2258
                  </property>
2261
                  </property>
2259
                  <property name="maximum">
2262
                  <property name="maximum">
2260
                   <double>1.000000000000000</double>
2263
                   <double>1.000000000000000</double>
3874
               <rect>
3877
               <rect>
3875
                <x>0</x>
3878
                <x>0</x>
3876
                <y>0</y>
3879
                <y>0</y>
3877
-               <width>411</width>
3878
-               <height>67</height>
3880
+               <width>96</width>
3881
+               <height>26</height>
3879
               </rect>
3882
               </rect>
3880
              </property>
3883
              </property>
3881
              <attribute name="label">
3884
              <attribute name="label">

+ 78
- 487
akvo/tressel/adapt.py View File

23
         with an adaptive filter that utilizes 2 reference signals instead of  
23
         with an adaptive filter that utilizes 2 reference signals instead of  
24
         the standard method which allows for only 1 reference signal.         
24
         the standard method which allows for only 1 reference signal.         
25
         Author: Rob Clemens              Date: 3/16/06                       
25
         Author: Rob Clemens              Date: 3/16/06                       
26
-        Modified and ported to Python, now takes arbitray number of reference points  
26
+        Modified and ported to Python, now takes arbitray number of reference points 
27
+        Original public domain source 
28
+        https://www.mathworks.com/matlabcentral/fileexchange/10447-noise-canceling-adaptive-filter 
29
+            x = data array 
30
+            R = reference array 
31
+            M = number of taps 
32
+            mu = forgetting factor 
33
+            PCA = Perform PCA 
27
         """
34
         """
28
         #from akvo.tressel import pca 
35
         #from akvo.tressel import pca 
29
         import akvo.tressel.pca as pca 
36
         import akvo.tressel.pca as pca 
46
             #H0 = H0[0:2*np.shape(x)[0]]
53
             #H0 = H0[0:2*np.shape(x)[0]]
47
 
54
 
48
         if all(H0) == 0:
55
         if all(H0) == 0:
56
+            # corrects for dimensionality issues if a simple 0 is passed
49
             H = np.zeros( (len(R)*M))
57
             H = np.zeros( (len(R)*M))
50
-            #print ("resetting filter")
51
         else:
58
         else:
52
             H = H0
59
             H = H0
53
-            #H = np.zeros( (len(R)*M))
54
 
60
 
55
         Rn = np.ones(len(R)*M) / mu 
61
         Rn = np.ones(len(R)*M) / mu 
56
         
62
         
57
         r_ = np.zeros( (len(R), M) ) 
63
         r_ = np.zeros( (len(R), M) ) 
58
-        e = np.zeros(len(x)) # error, desired output
64
+        e = np.zeros(len(x))           # error, in our case the desired output
59
         ilambda = lambda2**-1
65
         ilambda = lambda2**-1
60
 
66
 
61
-        for z in range(0, len(x)):
67
+        for ix in range(0, len(x)):
62
             # Only look forwards, to avoid distorting the lates times 
68
             # Only look forwards, to avoid distorting the lates times 
63
             # (run backwards, if opposite and you don't care about distorting very late time.)
69
             # (run backwards, if opposite and you don't care about distorting very late time.)
64
-            for ir in range(len(R)):
65
-                if z < M:
66
-                    r_[ir,0:z] = R[ir][0:z]
67
-                    r_[ir,z:M] = 0 
70
+            for ir in range(len(R)):  # number of reference channels 
71
+                if ix < M:
72
+                    r_[ir,0:ix] = R[ir][0:ix]
73
+                    r_[ir,ix:M] = 0 
68
                 else:
74
                 else:
69
-                    # TODO, use np.delete and np.append to speed this up
70
-                    r_[ir,:] = R[ir][z-M:z]
75
+                    r_[ir,:] = R[ir][ix-M:ix]
76
+
71
             # reshape            
77
             # reshape            
72
-            r_n = np.reshape(r_, -1) #concatenate((r_v, r_h ))
78
+            r_n = np.reshape(r_, -1) # concatenate the ref channels in to a 1D array 
73
 
79
 
74
-            #K      = np.dot( np.diag(Rn,0), r_n) / (lambda2 + np.dot(r_n*Rn, r_n))  # Create/update K
75
-            K      = (Rn* r_n) / (lambda2 + np.dot(r_n*Rn, r_n))  # Create/update K
76
-            e[z]   = x[z] - np.dot(r_n.T, H)             # e is the filtered signal, input - r(n) * Filter Coefs
77
-            H     += K*e[z];                             # Update Filter Coefficients
78
-            Rn     = ilambda*Rn - ilambda*np.dot(np.dot(K, r_n.T), Rn)     # Update R(n)
80
+            K      = (Rn* r_n) / (lambda2 + np.dot(r_n*Rn, r_n))       # Create/update K
81
+            e[ix]  = x[ix] - np.dot(r_n.T, H)                          # e is the filtered signal, input - r(n) * Filter Coefs
82
+            H     += K*e[ix];                                          # Update Filter Coefficients
83
+            Rn     = ilambda*Rn - ilambda*np.dot(np.dot(K, r_n.T), Rn) # Update R(n)
79
         return e, H
84
         return e, H
80
 
85
 
81
     def transferFunctionFFT(self, D, R, reg=1e-2):
86
     def transferFunctionFFT(self, D, R, reg=1e-2):
82
-        from akvo.tressel import pca
83
-        """
84
-            Computes the transfer function (H) between a Data channel and 
85
-            a number of Reference channels. The Matrices D and R are 
86
-            expected to be in the frequency domain on input.
87
-            | R1'R1 R1'R2 R1'R3|   |h1|   |R1'D|
88
-            | R2'R1 R2'R2 R2'R3| * |h2| = |R2'D|
89
-            | R3'R1 R3'R2 R3'R3|   |h3|   |R3'D|
90
-
91
-            Returns the corrected array 
92
-        """
93
-
94
-        # PCA decomposition on ref channels so signals are less related
95
-        #transMatrix, K, means = pca.pca( np.array([rx0, rx1]))   
96
-        #RR = np.zeros(( np.shape(R[0])[0]*np.shape(R[0])[1], len(R)))
97
-#         RR = np.zeros(( len(R), np.shape(R[0])[0]*np.shape(R[0])[1] ))
98
-#         for ir in range(len(R)):
99
-#             RR[ir,:] = np.reshape(R[ir], -1)
100
-#         transMatrix, K, means = pca.pca(RR)    
101
-#         #R rx0 = transMatrix[0,:]
102
-#         # rx1 = transMatrix[1,:]
103
-#         for ir in range(len(R)):
104
-#             R[ir] = transMatrix[ir,0]
105
-
106
-        import scipy.linalg 
107
-        import akvo.tressel.pca as pca 
108
-        # Compute as many transfer functions as len(R)
109
-        # A*H = B
110
-        nref = len(R)
111
-        H = np.zeros( (np.shape(D)[1], len(R)), dtype=complex )
112
-        for iw in range(np.shape(D)[1]):
113
-            A = np.zeros( (nref, nref), dtype=complex )
114
-            B = np.zeros( (nref) , dtype=complex)
115
-            for ii in range(nref):
116
-                for jj in range(nref):
117
-                    # build A
118
-                    A[ii,jj] = np.dot(R[ii][:,iw], R[jj][:,iw])                 
119
-
120
-                # build B
121
-                B[ii] = np.dot( R[ii][:,iw], D[:,iw] )
122
-
123
-            # compute H(iw)
124
-            #linalg.solve(a,b) if a is square
125
-            #print "A", A
126
-            #print "B", B
127
-            # TODO, regularise this solve step? So as to not fit the spurious noise
128
-            #print np.shape(B), np.shape(A) 
129
-            #H[iw, :] = scipy.linalg.solve(A,B)
130
-            H[iw, :] = scipy.linalg.lstsq(A,B,cond=reg)[0]
131
-            #print "lstqt", np.shape(scipy.linalg.lstsq(A,B))
132
-            #print "solve", scipy.linalg.solve(A,B)
133
-            #H[iw,:]  = scipy.linalg.lstsq(A,B) # otherwise 
134
-                #H = np.zeros( (np.shape(D)[1], )   )
135
-        #print H #A, B
136
-        Error = np.zeros(np.shape(D), dtype=complex)
137
-        for ir in range(nref):
138
-            for q in range( np.shape(D)[0] ):
139
-                #print "dimcheck", np.shape(H[:,ir]), np.shape(R[ir][q,:] )
140
-                Error[q,:] += H[:,ir]*R[ir][q,:]
141
-        return D - Error
142
-
143
-    def adapt_filt_tworefFreq(self, x, rx0, rx1, M, lambda2=0.95):
144
-        """ Frequency domain version of above
145
-        """
146
-        from akvo.tressel import pca 
147
-
148
-        pylab.figure()
149
-        pylab.plot(rx0)
150
-        pylab.plot(rx1)
151
-
152
-        # PCA decomposition on ref channels so signals are less related
153
-        transMatrix, K, means = pca.pca( np.array([rx0, rx1]))    
154
-        rx0 = transMatrix[:,0]
155
-        rx1 = transMatrix[:,1]
156
-        
157
-        pylab.plot(rx0)
158
-        pylab.plot(rx1)
159
-        pylab.show()
160
-        exit()
161
-
162
-        if np.shape(x) != np.shape(rx0) or np.shape(x) != np.shape(rx1):
163
-            print ("Error, non aligned")
164
-            exit(1)
165
-
166
-        wx = fft.rfft(x)
167
-        wr0 = fft.rfft(rx0)
168
-        wr1 = fft.rfft(rx1)
87
+         from akvo.tressel import pca
88
+         """
89
+             Computes the transfer function (H) between a Data channel and 
90
+             a number of Reference channels. The Matrices D and R are 
91
+             expected to be in the frequency domain on input.
92
+             | R1'R1 R1'R2 R1'R3|   |h1|   |R1'D|
93
+             | R2'R1 R2'R2 R2'R3| * |h2| = |R2'D|
94
+             | R3'R1 R3'R2 R3'R3|   |h3|   |R3'D|
169
  
95
  
170
-        H = np.zeros( (2*M), dtype=complex ) 
171
-        ident_mat = np.eye((2*M))
172
-        Rn = ident_mat / 0.1
173
-        r_v = np.zeros( (M), dtype=complex ) 
174
-        r_h = np.zeros( (M), dtype=complex ) 
175
-        e = np.zeros(len(x), dtype=complex )
176
-        ilambda = lambda2**-1
177
-
178
-        for z in range(0, len(wx)):
179
-            # TODO Padd with Zeros or truncate if M >,< arrays 
180
-            r_v = wr0[::-1][:M] 
181
-            r_h = wr1[::-1][:M] 
182
-            r_n = np.concatenate((r_v, r_h ))
183
-            K      = np.dot(Rn, r_n) / (lambda2 + np.dot(np.dot(r_n.T, Rn), r_n))  # Create/update K
184
-            e[z]   = wx[z] - np.dot(r_n,H)        # e is the filtered signal, input - r(n) * Filter Coefs
185
-            H     += K * e[z];                    # Update Filter Coefficients
186
-            Rn     = ilambda*Rn - ilambda*K*r_n.T*Rn  # Update R(n)
187
-        
188
-        return fft.irfft(e)
189
-
190
-    def iter_filt_refFreq(self, x, rx0, Ahat=.05, Bhat=.5, k=0.05):
191
-
192
-        X = np.fft.rfft(x)
193
-        X0 = np.copy(X)
194
-        RX0 = np.fft.rfft(rx0)
195
-
196
-        # step 0
197
-        Abs2HW = []
198
-        alphai =  k * (np.abs(Ahat)**2 / np.abs(Bhat)**2) 
199
-        betai  =  k * (1. / (np.abs(Bhat)**2) ) 
200
-        Hw     =  ((1.+alphai) * np.abs(X)**2 ) / (np.abs(X)**2 + betai*(np.abs(RX0)**2))
201
-        H      =  np.abs(Hw)**2
202
-        pylab.ion()
203
-        pylab.figure()
204
-        for i in range(10):
205
-            #print "alphai", alphai
206
-            #print "betai", betai
207
-            #print "Hw", Hw
208
-            alphai = k * (np.abs(Ahat)**2 / np.abs(Bhat)**2) * np.product(H, axis=0)
209
-            betai  = k * (1. / np.abs(Bhat)**2) * np.product(H, axis=0)
210
-            # update signal
211
-            Hw   =  ((1.+alphai) * np.abs(X)**2) / (np.abs(X)**2 + betai*np.abs(RX0)**2)
212
-            Hw = np.nan_to_num(Hw)
213
-            X *= Hw
214
-            H = np.vstack( (H, np.abs(Hw)**2) )
215
-            #print "Hw", Hw
216
-            pylab.cla()
217
-            pylab.plot(Hw)
218
-            #pylab.plot(np.abs(X))
219
-            #pylab.plot(np.abs(RX0))
220
-            pylab.draw()
221
-            raw_input("wait")
222
-
223
-        pylab.cla()
224
-        pylab.ioff()
225
-        #return np.fft.irfft(X0-X)
226
-        return np.fft.irfft(X)
227
-
228
-    def iter_filt_refFreq(self, x, rx0, rx1, Ahat=.1, Bhat=1., k=0.001):
229
-
230
-        X = np.fft.rfft(x)
231
-        X0 = np.copy(X)
232
-        RX0 = np.fft.rfft(rx0)
233
-        RX1 = np.fft.rfft(rx1)
234
-
235
-        # step 0
236
-        alphai =  k * (np.abs(Ahat)**2 / np.abs(Bhat)**2) 
237
-        betai  =  k * (1. / (np.abs(Bhat)**2) ) 
238
-        #Hw     =  ((1.+alphai) * np.abs(X)**2 ) / (np.abs(X)**2 + betai*(np.abs(RX0)**2))
239
-        H      =  np.ones(len(X)) # abs(Hw)**2
240
-        #pylab.ion()
241
-        #pylab.figure(334)
242
-        for i in range(1000):
243
-            #print "alphai", alphai
244
-            #print "betai", betai
245
-            #print "Hw", Hw
246
-            alphai = k * (np.abs(Ahat)**2 / np.abs(Bhat)**2) * np.product(H, axis=0)
247
-            betai  = k * (1. / np.abs(Bhat)**2) * np.product(H, axis=0)
248
-            # update signal
249
-            Hw   =  ((1.+alphai) * np.abs(X)**2) / (np.abs(X)**2 + betai*np.abs(RX0)**2)
250
-            Hw = np.nan_to_num(Hw)
251
-            X *= Hw #.conjugate
252
-            #H = np.vstack((H, np.abs(Hw)**2) )
253
-            H = np.vstack((H, np.abs(Hw)) )
254
-            #print "Hw", Hw
255
-            #pylab.cla()
256
-            #pylab.plot(Hw)
257
-            #pylab.plot(np.abs(X))
258
-            #pylab.plot(np.abs(RX0))
259
-            #pylab.draw()
260
-            #raw_input("wait")
261
-
262
-        #pylab.cla()
263
-        #pylab.ioff()
264
-        return np.fft.irfft(X0-X)
265
-        #return np.fft.irfft(X)
266
-
267
-    def Tdomain_DFT(self, desired, input, S):
268
-        """ Lifted from Adaptive filtering toolbox. Modefied to accept more than one input 
269
-            vector
270
-        """
271
-        
272
-        # Initialisation Procedure
273
-        nCoefficients =   S["filterOrderNo"]/2+1
274
-        nIterations   =   len(desired)
275
-
276
-        # Pre Allocations
277
-        errorVector  = np.zeros(nIterations, dtype='complex')
278
-        outputVector = np.zeros(nIterations, dtype='complex')
279
-        
280
-        # Initial State
281
-        coefficientVectorDFT =   np.fft.rfft(S["initialCoefficients"])/np.sqrt(float(nCoefficients))
282
-        desiredDFT           =   np.fft.rfft(desired)
283
-        powerVector          =   S["initialPower"]*np.ones(nCoefficients)
284
-
285
-        # Improve source code regularity, pad with zeros
286
-        # TODO, confirm zeros(nCoeffics) not nCoeffics-1
287
-        prefixedInput  =   np.concatenate([np.zeros(nCoefficients-1), np.array(input)])
288
-
289
-        # Body
290
-        pylab.ion()
291
-        pylab.figure(11)
292
-        for it in range(nIterations): # = 1:nIterations,
293
-            
294
-            regressorDFT = np.fft.rfft(prefixedInput[it:it+nCoefficients][::-1]) /\
295
-                           np.sqrt(float(nCoefficients))
296
-
297
-            # Summing two column vectors
298
-            powerVector = S["alpha"] * (regressorDFT*np.conjugate(regressorDFT)) + \
299
-                                  (1.-S["alpha"])*(powerVector)
300
-
301
-            pylab.cla()
302
-            #pylab.plot(prefixedInput[::-1], 'b')
303
-            #pylab.plot(prefixedInput[it:it+nCoefficients][::-1], 'g', linewidth=3)
304
-            #pylab.plot(regressorDFT.real)
305
-            #pylab.plot(regressorDFT.imag)
306
-            pylab.plot(powerVector.real)
307
-            pylab.plot(powerVector.imag)
308
-            #pylab.plot(outputVector)
309
-            #pylab.plot(errorVector.real)
310
-            #pylab.plot(errorVector.imag)
311
-            pylab.draw()
312
-            #raw_input("wait")
313
-
314
-            outputVector[it] = np.dot(coefficientVectorDFT.T, regressorDFT)
315
-
316
-            #errorVector[it] = desired[it] - outputVector[it]
317
-            errorVector[it] = desiredDFT[it] - outputVector[it]
318
-
319
-            #print errorVector[it], desired[it], outputVector[it]
320
-
321
-            # Vectorized
322
-            coefficientVectorDFT += (S["step"]*np.conjugate(errorVector[it])*regressorDFT) /\
323
-                                    (S['gamma']+powerVector)
324
-
325
-        return np.real(np.fft.irfft(errorVector))
326
-        #coefficientVector = ifft(coefficientVectorDFT)*sqrt(nCoefficients);
327
-
328
-    def Tdomain_DCT(self, desired, input, S):
329
-        """ Lifted from Adaptive filtering toolbox. Modefied to accept more than one input 
330
-            vector. Uses cosine transform
331
-        """
332
-        from scipy.fftpack import dct
96
+             Returns the corrected array 
97
+         """
333
  
98
  
334
-        # Initialisation Procedure
335
-        nCoefficients =   S["filterOrderNo"]+1
336
-        nIterations   =   len(desired)
337
-
338
-        # Pre Allocations
339
-        errorVector  = np.zeros(nIterations)
340
-        outputVector = np.zeros(nIterations)
341
-        
342
-        # Initial State
343
-        coefficientVectorDCT =   dct(S["initialCoefficients"]) #/np.sqrt(float(nCoefficients))
344
-        desiredDCT           =   dct(desired)
345
-        powerVector          =   S["initialPower"]*np.ones(nCoefficients)
346
-
347
-        # Improve source code regularity, pad with zeros
348
-        prefixedInput  =   np.concatenate([np.zeros(nCoefficients-1), np.array(input)])
349
-        
350
-        # Body
351
-        #pylab.figure(11)
352
-        #pylab.ion()
353
-        for it in range(0, nIterations): # = 1:nIterations,
354
-            
355
-            regressorDCT = dct(prefixedInput[it:it+nCoefficients][::-1], type=2) 
356
-            #regressorDCT = dct(prefixedInput[it+nCoefficients:it+nCoefficients*2+1])#[::-1]) 
357
-
358
-            # Summing two column vectors
359
-            powerVector = S["alpha"]*(regressorDCT) + (1.-S["alpha"])*(powerVector)
360
-            #pylab.cla()
361
-            #pylab.plot(powerVector)
362
-            #pylab.draw()
363
-
364
-            outputVector[it] = np.dot(coefficientVectorDCT.T, regressorDCT)
365
-            #errorVector[it] = desired[it] - outputVector[it]
366
-            errorVector[it] = desiredDCT[it] - outputVector[it]
367
-
368
-            # Vectorized
369
-            coefficientVectorDCT += (S["step"]*errorVector[it]*regressorDCT) #/\
370
-                                    #(S['gamma']+powerVector)
371
-
372
-        #pylab.plot(errorVector)
373
-        #pylab.show()
374
-        return dct(errorVector, type=3)
375
-        #coefficientVector = ifft(coefficientVectorDCT)*sqrt(nCoefficients);
376
-
377
-
378
-
379
-    def Tdomain_CORR(self, desired, input, S):
380
-
381
-        from scipy.linalg import toeplitz
382
-        from scipy.signal import correlate
383
-
384
-        # Autocorrelation
385
-        ac = np.correlate(input, input, mode='full')
386
-        ac = ac[ac.size/2:]
387
-        R = toeplitz(ac)
388
-        
389
-        # cross correllation
390
-        r = np.correlate(desired, input, mode='full')
391
-        r = r[r.size/2:]
392
-        
393
-        #r = np.correlate(desired, input, mode='valid')
394
-        print ("R", np.shape(R))
395
-        print ("r", np.shape(r))
396
-        print ("solving")
397
-        #H = np.linalg.solve(R,r)
398
-        H = np.linalg.lstsq(R,r,rcond=.01)[0]
399
-        #return desired - np.dot(H,input)
400
-        print ("done solving")
401
-        pylab.figure()
402
-        pylab.plot(H)
403
-        pylab.title("H")
404
-        #return desired - np.convolve(H, input, mode='valid')
405
-        #return desired - np.convolve(H, input, mode='same')
406
-        #return np.convolve(H, input, mode='same')
407
-        return desired - np.dot(toeplitz(H), input)
408
-        #return np.dot(R, H)
409
-
410
-#         T = toeplitz(input)
411
-#         print "shapes", np.shape(T), np.shape(desired)
412
-#         h = np.linalg.lstsq(T, desired)[0]
413
-#         print "shapes", np.shape(h), np.shape(input)
414
-#         #return np.convolve(h, input, mode='same')
415
-#         return desired - np.dot(T,h)
99
+         # PCA decomposition on ref channels so signals are less related
100
+         #transMatrix, K, means = pca.pca( np.array([rx0, rx1]))   
101
+         #RR = np.zeros(( np.shape(R[0])[0]*np.shape(R[0])[1], len(R)))
102
+ #         RR = np.zeros(( len(R), np.shape(R[0])[0]*np.shape(R[0])[1] ))
103
+ #         for ir in range(len(R)):
104
+ #             RR[ir,:] = np.reshape(R[ir], -1)
105
+ #         transMatrix, K, means = pca.pca(RR)    
106
+ #         #R rx0 = transMatrix[0,:]
107
+ #         # rx1 = transMatrix[1,:]
108
+ #         for ir in range(len(R)):
109
+ #             R[ir] = transMatrix[ir,0]
416
  
110
  
417
-    def Fdomain_CORR(self, desired, input, dt, freq):
418
-        
419
-        from scipy.linalg import toeplitz
420
-        
421
-        # Fourier domain
422
-        Input = np.fft.rfft(input)
423
-        Desired = np.fft.rfft(desired)
424
-
425
-        T = toeplitz(Input)
426
-        #H = np.linalg.solve(T, Desired)
427
-        H = np.linalg.lstsq(T, Desired)[0]
428
-#         ac = np.correlate(Input, Input, mode='full')
429
-#         ac = ac[ac.size/2:]
430
-#         R = toeplitz(ac)
431
-#         
432
-#         r = np.correlate(Desired, Input, mode='full')
433
-#         r = r[r.size/2:]
434
-#         
435
-#         #r = np.correlate(desired, input, mode='valid')
436
-#         print "R", np.shape(R)
437
-#         print "r", np.shape(r)
438
-#         print "solving"
439
-#         H = np.linalg.solve(R,r)
440
-#         #H = np.linalg.lstsq(R,r)
441
-#         #return desired - np.dot(H,input)
442
-#         print "done solving"
443
-        pylab.figure()
444
-        pylab.plot(H.real)
445
-        pylab.plot(H.imag)
446
-        pylab.plot(Input.real)
447
-        pylab.plot(Input.imag)
448
-        pylab.plot(Desired.real)
449
-        pylab.plot(Desired.imag)
450
-        pylab.legend(["hr","hi","ir","ii","dr","di"])
451
-        pylab.title("H")
452
-        #return desired - np.fft.irfft(Input*H)
453
-        return np.fft.irfft(H*Input)
454
-
455
-    def Tdomain_RLS(self, desired, input, S):
456
-        """
457
-            A DFT is first performed on the data. Than a RLS algorithm is carried out 
458
-            for noise cancellation. Related to the RLS_Alt Algoritm 5.3 in  Diniz book.
459
-            The desired and input signals are assummed to be real time series data.
460
-        """
461
-
462
-        # Transform data into frequency domain
463
-        Input = np.fft.rfft(input)
464
-        Desired = np.fft.rfft(desired)
465
-
466
-        # Initialisation Procedure
467
-        nCoefficients = S["filterOrderNo"]+1
468
-        nIterations   = len(Desired)
469
-
470
-        # Pre Allocations
471
-        errorVector  = np.zeros(nIterations, dtype="complex")
472
-        outputVector = np.zeros(nIterations, dtype="complex")
473
-        errorVectorPost  = np.zeros(nIterations, dtype="complex")
474
-        outputVectorPost = np.zeros(nIterations, dtype="complex")
475
-        coefficientVector = np.zeros( (nCoefficients, nIterations+1), dtype="complex" )        
476
-
477
-        # Initial State
478
-        coefficientVector[:,1] = S["initialCoefficients"]  
479
-        S_d                    = S["delta"]*np.eye(nCoefficients)
480
-
481
-        # Improve source code regularity, pad with zeros
482
-        prefixedInput = np.concatenate([np.zeros(nCoefficients-1, dtype="complex"), 
483
-                                np.array(Input)])
484
-        invLambda = 1./S["lambda"]
485
-        
486
-        # Body
487
-        pylab.ion()
488
-        pylab.figure(11)
489
-
490
-        for it in range(nIterations):
491
-            
492
-            regressor = prefixedInput[it:it+nCoefficients][::-1]
493
-
494
-            # a priori estimated output
495
-            outputVector[it] = np.dot(coefficientVector[:,it].T, regressor)
496
-       
497
-            # a priori error
498
-            errorVector[it] = Desired[it] - outputVector[it]
499
-
500
-            psi             = np.dot(S_d, regressor)
501
-            if np.isnan(psi).any():
502
-                print ("psi", psi)
503
-                exit(1)
504
-            
505
-            pylab.cla()
506
-            #pylab.plot(psi)
507
-            pylab.plot(regressor.real)
508
-            pylab.plot(regressor.imag)
509
-            pylab.plot(coefficientVector[:,it].real)
510
-            pylab.plot(coefficientVector[:,it].imag)
511
-            pylab.legend(["rr","ri", "cr", "ci"])
512
-            pylab.draw()
513
-            raw_input("paws")
514
-
515
-            S_d             = invLambda * (S_d - np.dot(psi, psi.T)  /\
516
-                                S["lambda"] + np.dot(psi.T, regressor))
517
-
518
-            coefficientVector[:,it+1] = coefficientVector[:,it] + \
519
-                                        np.conjugate(errorVector[it])*np.dot(S_d, regressor)
520
-            # A posteriori estimated output
521
-            outputVectorPost[it]  =  np.dot(coefficientVector[:,it+1].T, regressor)
522
-
523
-            # A posteriori error
524
-            errorVectorPost[it] = Desired[it] - outputVectorPost[it]
111
+         import scipy.linalg 
112
+         import akvo.tressel.pca as pca 
113
+         # Compute as many transfer functions as len(R)
114
+         # A*H = B
115
+         nref = len(R)
116
+         H = np.zeros( (np.shape(D)[1], len(R)), dtype=complex )
117
+         for iw in range(np.shape(D)[1]):
118
+             A = np.zeros( (nref, nref), dtype=complex )
119
+             B = np.zeros( (nref) , dtype=complex)
120
+             for ii in range(nref):
121
+                 for jj in range(nref):
122
+                     # build A
123
+                     A[ii,jj] = np.dot(R[ii][:,iw], R[jj][:,iw])                 
525
  
124
  
526
-        errorVectorPost = np.nan_to_num(errorVectorPost)
527
-
528
-        pylab.figure(11)
529
-        print (np.shape(errorVectorPost))
530
-        pylab.plot(errorVectorPost.real)
531
-        pylab.plot(errorVectorPost.imag)
532
-        pylab.show()
533
-        print(errorVectorPost)
534
-        #return np.fft.irfft(Desired)
535
-        return np.fft.irfft(errorVectorPost)
536
-
537
-if __name__ == "__main__":
538
-
539
-    def noise(nu, t, phi):
540
-        return np.sin(nu*2.*np.pi*t + phi)
125
+                 # build B
126
+                 B[ii] = np.dot( R[ii][:,iw], D[:,iw] )
127
+ 
128
+             # compute H(iw)
129
+             #linalg.solve(a,b) if a is square
130
+             #print "A", A
131
+             #print "B", B
132
+             # TODO, regularise this solve step? So as to not fit the spurious noise
133
+             #print np.shape(B), np.shape(A) 
134
+             #H[iw, :] = scipy.linalg.solve(A,B)
135
+             H[iw, :] = scipy.linalg.lstsq(A,B,cond=reg)[0]
136
+             #print "lstqt", np.shape(scipy.linalg.lstsq(A,B))
137
+             #print "solve", scipy.linalg.solve(A,B)
138
+             #H[iw,:]  = scipy.linalg.lstsq(A,B) # otherwise 
139
+                 #H = np.zeros( (np.shape(D)[1], )   )
140
+         #print H #A, B
141
+         Error = np.zeros(np.shape(D), dtype=complex)
142
+         for ir in range(nref):
143
+             for q in range( np.shape(D)[0] ):
144
+                 #print "dimcheck", np.shape(H[:,ir]), np.shape(R[ir][q,:] )
145
+                 Error[q,:] += H[:,ir]*R[ir][q,:]
146
+         return D - Error
541
 
147
 
542
-    import matplotlib.pyplot as plt
543
-    print("Test driver for adaptive filtering")
544
-    Filt = AdaptiveFilter(.1)
545
-    t = np.arange(0, .5, 1e-4)
546
-    omega = 2000 * 2.*np.pi
547
-    T2 = .100
548
-    n1 = noise(60, t, .2   )
549
-    n2 = noise(61, t, .514 )
550
-    x = np.sin(omega*t)* np.exp(-t/T2) + 2.3*noise(60, t, .34) + 1.783*noise(31, t, 2.1)
551
-    e = Filt.adapt_filt_tworef(x, n1, n2, 200, .98)
552
-    plt.plot(t,  x)
553
-    plt.plot(t, n1)
554
-    plt.plot(t, n2)
555
-    plt.plot(t,  e)
556
-    plt.show()

+ 2
- 2
akvo/tressel/harmonic.py View File

43
     if Nsearch != False:
43
     if Nsearch != False:
44
         kNs = k1+Nsearch    
44
         kNs = k1+Nsearch    
45
     if Bounds == 0:
45
     if Bounds == 0:
46
-        print("UNbounded search from ", k1, " to ", kNs) # for f0 with fN=10 in search", f0)
47
         # CG, BFGS, Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov, trust-exact and trust-constr
46
         # CG, BFGS, Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov, trust-exact and trust-constr
48
         res = minimize(harmonicNorm, np.array((f0)), args=(sN, fs, t, k1, kNs, ks), jac='2-point', method='BFGS') # hess=None, bounds=None )
47
         res = minimize(harmonicNorm, np.array((f0)), args=(sN, fs, t, k1, kNs, ks), jac='2-point', method='BFGS') # hess=None, bounds=None )
48
+        print("UNbounded search from ", k1, " to ", kNs, res.x[0]) # for f0 with fN=10 in search", f0)
49
     
49
     
50
     else:
50
     else:
51
         bnds = ( (f0-Bounds, f0+Bounds), )
51
         bnds = ( (f0-Bounds, f0+Bounds), )
52
-        print("bounded ( +-", Bounds, ") search from ", k1, "to", kNs) # for f0 with fN=10 in search", f0)
53
         res = minimize(harmonicNorm, (f0,), args=(sN, fs, t, k1, kNs, ks), jac='2-point', method='L-BFGS-B', bounds=bnds ) # hess=None, bounds=None )
52
         res = minimize(harmonicNorm, (f0,), args=(sN, fs, t, k1, kNs, ks), jac='2-point', method='L-BFGS-B', bounds=bnds ) # hess=None, bounds=None )
53
+        print("bounded ( +-", Bounds, ") search from ", k1, "to", kNs, res.x[0]) # for f0 with fN=10 in search", f0)
54
 
54
 
55
     return harmonicEuler(sN, fs, t, res.x[0], k1, kN, ks), res.x[0]#[0]
55
     return harmonicEuler(sN, fs, t, res.x[0], k1, kN, ks), res.x[0]#[0]
56
 
56
 

+ 53
- 17
akvo/tressel/mrsurvey.py View File

19
 import multiprocessing 
19
 import multiprocessing 
20
 import itertools 
20
 import itertools 
21
 
21
 
22
+import padasip as pa
23
+
22
 import akvo.tressel.adapt as adapt
24
 import akvo.tressel.adapt as adapt
23
 #import akvo.tressel.cadapt as adapt # cython for more faster
25
 #import akvo.tressel.cadapt as adapt # cython for more faster
24
 import akvo.tressel.decay as decay
26
 import akvo.tressel.decay as decay
608
                                     X = np.fft.rfft(self.DATADICT[pulse][ichan][ipm][istack])
610
                                     X = np.fft.rfft(self.DATADICT[pulse][ichan][ipm][istack])
609
                                     canvas.ax1.plot(np.abs(ww[0:len(X)]), np.abs(X), alpha=.5)
611
                                     canvas.ax1.plot(np.abs(ww[0:len(X)]), np.abs(X), alpha=.5)
610
                             canvas.ax1.set_prop_cycle(None)
612
                             canvas.ax1.set_prop_cycle(None)
611
-                            #canvas.ax1.set_ylim(-mmaxr, mmaxr) 
613
+                            canvas.ax1.set_ylim(-mmaxr, mmaxr) 
612
                         for ichan in self.DATADICT[pulse]["chan"]:
614
                         for ichan in self.DATADICT[pulse]["chan"]:
613
                             if TDPlot:
615
                             if TDPlot:
614
                                 canvas.ax2.plot( self.DATADICT[pulse]["TIMES"], 1e9*self.DATADICT[pulse][ichan][ipm][istack], alpha=.5) 
616
                                 canvas.ax2.plot( self.DATADICT[pulse]["TIMES"], 1e9*self.DATADICT[pulse][ichan][ipm][istack], alpha=.5) 
618
                                 X = np.fft.rfft(self.DATADICT[pulse][ichan][ipm][istack])
620
                                 X = np.fft.rfft(self.DATADICT[pulse][ichan][ipm][istack])
619
                                 canvas.ax2.plot(np.abs(ww[0:len(X)]), np.abs(X), alpha=.5)
621
                                 canvas.ax2.plot(np.abs(ww[0:len(X)]), np.abs(X), alpha=.5)
620
                         canvas.ax2.set_prop_cycle(None)
622
                         canvas.ax2.set_prop_cycle(None)
621
-                        #canvas.ax2.set_ylim(-mmaxd, mmaxd)
623
+                        canvas.ax2.set_ylim(-mmaxd, mmaxd)
622
                     if procRefs: 
624
                     if procRefs: 
623
                         for ichan in self.DATADICT[pulse]["rchan"]:
625
                         for ichan in self.DATADICT[pulse]["rchan"]:
624
                             if nF == 1:
626
                             if nF == 1:
1621
 
1623
 
1622
                     # calculate effective current/moment
1624
                     # calculate effective current/moment
1623
                     I0 = np.abs(X)/len(X) 
1625
                     I0 = np.abs(X)/len(X) 
1624
-                    qeff = I0 #* (self.DATADICT[pulse]["PULSE_TIMES"][-1]-self.DATADICT[pulse]["PULSE_TIMES"][0])
1626
+                    qeff = I0 * (self.DATADICT[pulse]["PULSE_TIMES"][-1]-self.DATADICT[pulse]["PULSE_TIMES"][0])
1625
 
1627
 
1626
                     # frequency plot
1628
                     # frequency plot
1627
                     #canvas.ax1.set_title(r"pulse moment index " +str(ipm), fontsize=10)
1629
                     #canvas.ax1.set_title(r"pulse moment index " +str(ipm), fontsize=10)
1708
         for pulse in self.DATADICT["PULSES"]:
1710
         for pulse in self.DATADICT["PULSES"]:
1709
             H[pulse] = {}
1711
             H[pulse] = {}
1710
             for ichan in self.DATADICT[pulse]["chan"]:
1712
             for ichan in self.DATADICT[pulse]["chan"]:
1711
-                H[pulse][ichan] = np.zeros(M)
1713
+                H[pulse][ichan] = np.zeros(M*len( self.DATADICT[pulse]["rchan"] ))
1712
         
1714
         
1713
         iFID = 0 
1715
         iFID = 0 
1714
         # original ordering... 
1716
         # original ordering... 
1726
                         mmax = max(mmax, np.max(1e9*self.DATADICT[pulse][ichan][ipm][istack])) 
1728
                         mmax = max(mmax, np.max(1e9*self.DATADICT[pulse][ichan][ipm][istack])) 
1727
                     canvas.ax2.set_ylim(-mmax, mmax) 
1729
                     canvas.ax2.set_ylim(-mmax, mmax) 
1728
                     canvas.ax2.set_prop_cycle(None)
1730
                     canvas.ax2.set_prop_cycle(None)
1729
-
1730
                     for ichan in self.DATADICT[pulse]["chan"]:
1731
                     for ichan in self.DATADICT[pulse]["chan"]:
1731
                         #H = np.zeros(M)
1732
                         #H = np.zeros(M)
1732
                         RX = []
1733
                         RX = []
1733
                         for irchan in self.DATADICT[pulse]["rchan"]:
1734
                         for irchan in self.DATADICT[pulse]["rchan"]:
1734
                             RX.append(self.DATADICT[pulse][irchan][ipm][istack][::-1])
1735
                             RX.append(self.DATADICT[pulse][irchan][ipm][istack][::-1])
1735
-                        if all(H[pulse][ichan]) == 0:
1736
-                            # call twice to reconcile filter wind-up
1737
-                            [e,H[pulse][ichan]] = Filt.adapt_filt_Ref( self.DATADICT[pulse][ichan][ipm][istack][::-1],\
1738
-                                                        RX,\
1739
-                                                        M, mu, PCA, flambda, H[pulse][ichan])
1740
-                            [e,H[pulse][ichan]] = Filt.adapt_filt_Ref( self.DATADICT[pulse][ichan][ipm][istack][::-1],\
1736
+                        # Reset each time? 
1737
+                        #H[pulse][ichan] *= 0
1738
+                        #if all(H[pulse][ichan]) == 0:
1739
+                        if False: 
1740
+                            ####################################################################################
1741
+                            # Padasip adaptive filter implimentations, do not allow for variable filter length
1742
+                            ####################################################################################
1743
+                            # identification                                                                   #
1744
+                            #f = pa.filters.FilterRLS(n=len(self.DATADICT[pulse]["rchan"]), mu=0.99, w="zeros") #
1745
+                            #f = pa.filters.FilterGNGD(n=len(self.DATADICT[pulse]["rchan"]), mu=0.1)           #                          # Nope
1746
+                            #f = pa.filters.FilterLMS(n=len(self.DATADICT[pulse]["rchan"]), mu=0.1)            #                          # NOPE
1747
+                            #f = pa.filters.AdaptiveFilter(model="NLMS", n=len(self.DATADICT[pulse]["rchan"]), mu=0.1, w="random")        # NOPE  
1748
+                            #f = pa.filters.AdaptiveFilter(model="GNGD", n=len(self.DATADICT[pulse]["rchan"]), mu=0.1)                    # horrendous
1749
+                            #f = pa.filters.FilterNLMF(n=len(self.DATADICT[pulse]["rchan"]), mu=0.005, w="random")                        # BAD
1750
+                            
1751
+                            #f = pa.filters.FilterSSLMS(n=len(self.DATADICT[pulse]["rchan"]), mu=0.01, w="zeros")                         # pretty good
1752
+                            f = pa.filters.FilterNSSLMS(n=len(self.DATADICT[pulse]["rchan"]), mu=0.1, w="zeros")                          # pretty good 
1753
+
1754
+                            y, e, H[pulse][ichan] = f.run(self.DATADICT[pulse][ichan][ipm][istack][::-1], np.array(RX).T)    #
1755
+                            ####################################################################################
1756
+                            e = self.DATADICT[pulse][ichan][ipm][istack][::-1] - y
1757
+                        elif True:
1758
+                            # check for change in filter coefficients and rerun if things are changing too rapidly, 
1759
+                            #       this is especially true for the first run 
1760
+                            hm1 = np.copy(H[pulse][ichan]) 
1761
+                            [e, H[pulse][ichan]] = Filt.adapt_filt_Ref( self.DATADICT[pulse][ichan][ipm][istack][::-1],\
1762
+                                                     RX,\
1763
+                                                     M, mu, PCA, flambda, H[pulse][ichan])
1764
+                            iloop = 0
1765
+                            #while False: 
1766
+                            while (np.linalg.norm( H[pulse][ichan] - hm1) > .05):
1767
+                                hm1 = np.copy(H[pulse][ichan]) 
1768
+                                [e, H[pulse][ichan]] = Filt.adapt_filt_Ref( self.DATADICT[pulse][ichan][ipm][istack][::-1],\
1741
                                                         RX,\
1769
                                                         RX,\
1742
                                                         M, mu, PCA, flambda, H[pulse][ichan])
1770
                                                         M, mu, PCA, flambda, H[pulse][ichan])
1771
+                                iloop += 1
1772
+                            print("Recalled ", iloop, "times with norm=", np.linalg.norm(hm1-H[pulse][ichan]))
1743
                         else:
1773
                         else:
1744
                             [e,H[pulse][ichan]] = Filt.adapt_filt_Ref( self.DATADICT[pulse][ichan][ipm][istack][::-1],\
1774
                             [e,H[pulse][ichan]] = Filt.adapt_filt_Ref( self.DATADICT[pulse][ichan][ipm][istack][::-1],\
1745
                                                         RX,\
1775
                                                         RX,\
1753
                             canvas.ax2.plot( self.DATADICT[pulse]["TIMES"], 1e9* e[::-1],\
1783
                             canvas.ax2.plot( self.DATADICT[pulse]["TIMES"], 1e9* e[::-1],\
1754
                                 label = pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " ichan="  + str(ichan))
1784
                                 label = pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " ichan="  + str(ichan))
1755
                             self.DATADICT[pulse][ichan][ipm][istack] = e[::-1]
1785
                             self.DATADICT[pulse][ichan][ipm][istack] = e[::-1]
1756
-                            
1757
-                        canvas.ax1.plot( H[pulse][ichan] ) # , label="taps") 
1786
+                           
1787
+     
1788
+                        canvas.ax1.plot( H[pulse][ichan].reshape(-1, len(RX)) ) # , label="taps") 
1789
+
1758
                         canvas.ax2.legend(prop={'size':6}, loc='upper right')
1790
                         canvas.ax2.legend(prop={'size':6}, loc='upper right')
1759
                         #canvas.ax2.legend(prop={'size':6}, loc='upper right')
1791
                         #canvas.ax2.legend(prop={'size':6}, loc='upper right')
1760
 
1792
 
1793
+                        mh = np.max(np.abs( H[pulse][ichan] ))
1794
+                        canvas.ax1.set_ylim( -mh, mh )
1795
+
1761
                         canvas.ax2.set_xlabel(r"time [s]", fontsize=8)
1796
                         canvas.ax2.set_xlabel(r"time [s]", fontsize=8)
1762
                         canvas.ax2.set_ylabel(r"signal [nV]", fontsize=8)
1797
                         canvas.ax2.set_ylabel(r"signal [nV]", fontsize=8)
1763
 
1798
 
1930
         Filt = adapt.AdaptiveFilter(0.)
1965
         Filt = adapt.AdaptiveFilter(0.)
1931
         for pulse in self.DATADICT["PULSES"]:
1966
         for pulse in self.DATADICT["PULSES"]:
1932
             # Compute window function and dimensions      
1967
             # Compute window function and dimensions      
1933
-            [WINDOW, nd, wstart, wend, dead] = self.computeWindow(pulse, band, centre, ftype) 
1968
+            [WINDOW, nd, wstart, wend, dead, idead] = self.computeWindow(pulse, band, centre, ftype) 
1934
             for istack in self.DATADICT["stacks"]:
1969
             for istack in self.DATADICT["stacks"]:
1935
                 for ichan in self.DATADICT[pulse]["chan"]:
1970
                 for ichan in self.DATADICT[pulse]["chan"]:
1936
                     # FFT of stack
1971
                     # FFT of stack
2521
 
2556
 
2522
                     canvas.ax3.legend(prop={'size':6}, loc='upper right')
2557
                     canvas.ax3.legend(prop={'size':6}, loc='upper right')
2523
                     canvas.ax2.legend(prop={'size':6}, loc='upper right')
2558
                     canvas.ax2.legend(prop={'size':6}, loc='upper right')
2524
-                    
2525
-                    canvas.ax1.set_title("stack "+str(istack)+" pulse index " + str(ipm), fontsize=8)
2559
+
2526
                     canvas.ax1.set_ylabel("Current [A]", fontsize=8) 
2560
                     canvas.ax1.set_ylabel("Current [A]", fontsize=8) 
2527
                     canvas.ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y') 
2561
                     canvas.ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y') 
2562
+                    canvas.ax1.set_xlabel("time [s]", fontsize=8)
2528
                     
2563
                     
2564
+                    canvas.ax2.set_title("stack "+str(istack)+" pulse index " + str(ipm), fontsize=8)
2529
                     canvas.ax2.set_ylabel("RAW signal [V]", fontsize=8)
2565
                     canvas.ax2.set_ylabel("RAW signal [V]", fontsize=8)
2530
                     canvas.ax3.set_ylabel("RAW signal [V]", fontsize=8)
2566
                     canvas.ax3.set_ylabel("RAW signal [V]", fontsize=8)
2531
                     canvas.ax2.tick_params(axis='both', which='major', labelsize=8)
2567
                     canvas.ax2.tick_params(axis='both', which='major', labelsize=8)
2532
                     canvas.ax2.tick_params(axis='both', which='minor', labelsize=6)
2568
                     canvas.ax2.tick_params(axis='both', which='minor', labelsize=6)
2533
-                    canvas.ax2.set_xlabel("time [s]", fontsize=8)
2569
+
2534
                     canvas.draw()
2570
                     canvas.draw()
2535
 
2571
 
2536
             percent = (int) (1e2*((float)((iistack*self.nPulseMoments+ipm+1))  / (len(procStacks)*self.nPulseMoments)))
2572
             percent = (int) (1e2*((float)((iistack*self.nPulseMoments+ipm+1))  / (len(procStacks)*self.nPulseMoments)))

+ 4
- 0
pyuic.json View File

7
         [
7
         [
8
             "akvo/gui/addCircularLoop.ui",
8
             "akvo/gui/addCircularLoop.ui",
9
             "akvo/gui"
9
             "akvo/gui"
10
+        ],
11
+        [
12
+            "akvo/gui/addFigure8Loop.ui",
13
+            "akvo/gui"
10
         ]
14
         ]
11
     ],
15
     ],
12
     "pyrcc": "pyrcc5",
16
     "pyrcc": "pyrcc5",

+ 5
- 3
setup.py View File

21
     long_description = fh.read()
21
     long_description = fh.read()
22
 
22
 
23
 setup(name='Akvo',
23
 setup(name='Akvo',
24
-      version='1.3.0',
24
+      version='1.3.1',
25
+      python_requires='>3.7.0', # due to pyLemma
25
       description='Surface nuclear magnetic resonance workbench',
26
       description='Surface nuclear magnetic resonance workbench',
26
       long_description=long_description,
27
       long_description=long_description,
27
       long_description_content_type='text/markdown',
28
       long_description_content_type='text/markdown',
43
 #          'rpy2',
44
 #          'rpy2',
44
           'matplotlib',
45
           'matplotlib',
45
           'scipy',
46
           'scipy',
47
+          'padasip',
46
           'numpy',
48
           'numpy',
47
           'pyqt5',
49
           'pyqt5',
48
           'pyyaml',
50
           'pyyaml',
68
       # Mechanism to include auxiliary files
70
       # Mechanism to include auxiliary files
69
       include_package_data=True,
71
       include_package_data=True,
70
       package_data={
72
       package_data={
71
-        'akvo.gui': ['*.png'],  #All .r files 
72
-        'akvo.lemma': ['pyLemmaCore.so']
73
+        'akvo.gui': ['*.png'], 
74
+        'akvo.gui': ['*.ui']
73
       },
75
       },
74
       classifiers=[
76
       classifiers=[
75
         "Programming Language :: Python :: 3",
77
         "Programming Language :: Python :: 3",

Loading…
Cancel
Save