Surface NMR processing and inversion GUI
Du kannst nicht mehr als 25 Themen auswählen Themen müssen mit entweder einem Buchstaben oder einer Ziffer beginnen. Sie können Bindestriche („-“) enthalten und bis zu 35 Zeichen lang sein.

mrsurvey.py 110KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116
  1. from PyQt5.QtCore import *
  2. import numpy as np
  3. import scipy.signal as signal
  4. import pylab
  5. import sys
  6. import scipy
  7. import copy
  8. from scipy.io.matlab import mio
  9. from numpy import pi
  10. from math import floor
  11. import matplotlib as mpl
  12. from matplotlib.ticker import FuncFormatter
  13. import matplotlib.font_manager as fm
  14. import matplotlib.pyplot as plt
  15. import matplotlib.ticker
  16. from matplotlib.ticker import MaxNLocator
  17. import akvo.tressel.adapt as adapt
  18. #import akvo.tressel.cadapt as adapt # cython for more faster
  19. import akvo.tressel.decay as decay
  20. import akvo.tressel.pca as pca
  21. import akvo.tressel.rotate as rotate
  22. import akvo.tressel.cmaps as cmaps
  23. import cmocean # colormaps for geophysical data
  24. plt.register_cmap(name='viridis', cmap=cmaps.viridis)
  25. plt.register_cmap(name='inferno', cmap=cmaps.inferno)
  26. plt.register_cmap(name='inferno_r', cmap=cmaps.inferno_r)
  27. plt.register_cmap(name='magma', cmap=cmaps.magma)
  28. plt.register_cmap(name='magma_r', cmap=cmaps.magma_r)
  29. class SNMRDataProcessor(QObject):
  30. """ Revised class for preprocessing sNMR Data.
  31. Derived types can read GMR files
  32. """
  33. def __init__(self):
  34. QObject.__init__(self)
  35. self.numberOfMoments = 0
  36. self.numberOfPulsesPerMoment = 0
  37. self.pulseType = "NONE"
  38. self.transFreq = 0
  39. self.pulseLength = np.zeros(1)
  40. self.nPulseMoments = 0
  41. self.dt = 0
  42. def mfreqz(self, b,a=1):
  43. """ Plots the frequency response of a filter specified with a and b weights
  44. """
  45. import scipy.signal as signal
  46. pylab.figure(124)
  47. w,h = signal.freqz(b,a)
  48. w /= max(w)
  49. w *= .5/self.dt
  50. h_dB = 20 * pylab.log10 (abs(h))
  51. pylab.subplot(211)
  52. #pylab.plot(w/max(w),h_dB)
  53. pylab.plot(w,h_dB)
  54. pylab.ylim(-150, 5)
  55. pylab.ylabel('Magnitude (dB)')
  56. #pylab.xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
  57. pylab.xlabel(r'Hz')
  58. pylab.title(r'Frequency response')
  59. pylab.subplot(212)
  60. h_Phase = pylab.unwrap(pylab.arctan2(pylab.imag(h), pylab.real(h)))
  61. #pylab.plot(w/max(w),h_Phase)
  62. pylab.plot(w,h_Phase)
  63. pylab.ylabel('Phase (radians)')
  64. pylab.xlabel(r'Hz')
  65. #pylab.xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
  66. pylab.title(r'Phase response')
  67. pylab.subplots_adjust(hspace=0.5)
  68. def mfreqz2(self, b, a, canvas):
  69. "for analysing filt-filt"
  70. import scipy.signal as signal
  71. canvas.reAx2(False,False)
  72. canvas.ax1.tick_params(axis='both', which='major', labelsize=8)
  73. canvas.ax2.tick_params(axis='both', which='major', labelsize=8)
  74. #canvas.ax2.tick_params(axis='both', which='minor', labelsize=6)
  75. #pylab.figure(124)
  76. w,h = signal.freqz(b,a)
  77. w /= max(w)
  78. w *= .5/self.dt
  79. h_dB = 20 * pylab.log10(abs(h*h) + 1e-16)
  80. #ab.subplot(211)
  81. #pylab.plot(w/max(w),h_dB)
  82. canvas.ax1.plot(w,h_dB)
  83. canvas.ax1.set_ylim(-150, 5)
  84. canvas.ax1.set_ylabel('Magnitude [db]', fontsize=8)
  85. #pylab.xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
  86. canvas.ax1.set_xlabel(r'[Hz]', fontsize=8)
  87. canvas.ax1.set_title(r'Frequency response', fontsize=8)
  88. canvas.ax1.grid(True)
  89. tt = np.arange(0, .02, self.dt)
  90. impulse = signal.dimpulse((self.filt_z, self.filt_p, self.filt_k, self.dt), t=tt)
  91. #impulse = signal.dstep((self.filt_z, self.filt_p, self.filt_k, self.dt), t=tt)
  92. #print impulse
  93. #for ii in range(len(impulse[1])):
  94. impulse_dB = 20.*np.log10(np.abs(np.array(impulse[1][0])))
  95. #canvas.ax2.plot(np.array(impulse[0]), impulse_dB)
  96. canvas.ax2.plot(np.array(impulse[0]), impulse[1][0])
  97. #h_Phase = pylab.unwrap(pylab.arctan2(pylab.imag(h), pylab.real(h)))
  98. #canvas.ax2.plot(w,h_Phase)
  99. canvas.ax2.set_ylabel('response [%]', fontsize=8)
  100. canvas.ax2.set_xlabel(r'time [s]', fontsize=8)
  101. canvas.ax2.set_title(r'impulse response', fontsize=8)
  102. #canvas.ax2.grid(True)
  103. canvas.draw()
  104. # search for last
  105. return impulse #[np.where(impulse[1][0] > .01)[-1]]
  106. class GMRDataProcessor(SNMRDataProcessor):
  107. # slots
  108. progressTrigger = pyqtSignal("int")
  109. doneTrigger = pyqtSignal()
  110. enableDSPTrigger = pyqtSignal()
  111. updateProcTrigger = pyqtSignal()
  112. def __init__(self):
  113. SNMRDataProcessor.__init__(self)
  114. self.maxBusV = 0.
  115. self.samp = 50000. # sampling frequency
  116. self.dt = 2e-5 # sampling rate
  117. self.deadTime = .0055 # instrument dead time before measurement
  118. self.prePulseDelay = 0.05 # delay before pulse
  119. self.windead = 0. # FD window filter dead time
  120. self.pulseType = -1
  121. self.transFreq = -1
  122. self.maxBusV = -1
  123. self.pulseLength = -1
  124. self.interpulseDelay = -1 # for T2, Spin Echo
  125. self.repetitionDelay = -1 # delay between first pulse
  126. self.nPulseMoments = -1 # Number of pulse moments per stack
  127. self.TuneCapacitance = -1 # tuning capac in uF
  128. self.nTransVersion = -1 # Transmitter version
  129. self.nDAQVersion = -1 # DAQ software version
  130. self.nInterleaves = -1 # num interleaves
  131. # self.nReceiveChannels = 4 # Num receive channels
  132. self.RotatedAmplitude = False
  133. # self.DATA = np.zeros(1) # Numpy array to hold all data, dimensions resized based on experiment
  134. # self.PULSES = np.zeros(1) # Numpy array to hold all data, dimensions resized based on experiment
  135. def Print(self):
  136. print ("pulse type", self.pulseType)
  137. print ("maxBusV", self.maxBusV)
  138. print ("inner pulse delay", self.interpulseDelay)
  139. print ("tuning capacitance", self.TuneCapacitance)
  140. print ("sampling rate", self.samp)
  141. print ("dt", self.dt)
  142. print ("dead time", self.deadTime)
  143. print ("pre pulse delay", self.prePulseDelay)
  144. print ("number of pulse moments", self.nPulseMoments)
  145. print ("pulse Length", self.pulseLength)
  146. print ("trans freq", self.transFreq)
  147. def readHeaderFile(self, FileName):
  148. HEADER = np.loadtxt(FileName)
  149. pulseTypeDict = {
  150. 1 : lambda: "FID",
  151. 2 : lambda: "T1",
  152. 3 : lambda: "SPINECHO",
  153. 4 : lambda: "4PhaseT1"
  154. }
  155. pulseLengthDict = {
  156. 1 : lambda x: np.ones(1) * x,
  157. 2 : lambda x: np.ones(2) * x,
  158. 3 : lambda x: np.array([x, 2.*x]),
  159. 4 : lambda x: np.ones(2) * x
  160. }
  161. self.pulseType = pulseTypeDict.get((int)(HEADER[0]))()
  162. self.transFreq = HEADER[1]
  163. self.maxBusV = HEADER[2]
  164. self.pulseLength = pulseLengthDict.get((int)(HEADER[0]))(1e-3*HEADER[3])
  165. self.interpulseDelay = 1e-3*HEADER[4] # for T2, Spin Echo
  166. self.repetitionDelay = HEADER[5] # delay between first pulse
  167. self.nPulseMoments = (int)(HEADER[6]) # Number of pulse moments per stack
  168. self.TuneCapacitance = HEADER[7] # tuning capacitance in uF
  169. self.nTransVersion = HEADER[8] # Transmitter version
  170. self.nDAQVersion = HEADER[9] # DAQ software version
  171. self.nInterleaves = HEADER[10] # num interleaves
  172. self.gain()
  173. # default
  174. self.samp = 50000. # sampling frequency
  175. self.dt = 2e-5 # sampling rate
  176. # newer header files contain 64 entries
  177. if self.nDAQVersion >= 2:
  178. #print "new file format"
  179. self.samp = HEADER[14] # sampling frequency
  180. self.dt = 1./self.samp # sampling rate
  181. self.deadTime = .0055 # instrument dead time before measurement
  182. self.prePulseDelay = 0.05 # delay before pulse
  183. #exit()
  184. def gain(self):
  185. #######################################################
  186. # Circuit gain
  187. # From MRSMatlab
  188. w = 2*np.pi*self.transFreq
  189. # 1e6 due to uF of reported capacitance
  190. L_coil = 1e6/(self.TuneCapacitance*(w**2))
  191. R_coil = 1.
  192. Z1_in = .5 + 1j*.5*w
  193. Z2_in = 1./(1j*w*.000001616)
  194. Z_eq_inv = (1./Z1_in) + (1./Z2_in)
  195. Zeq = 1./Z_eq_inv
  196. Zsource = R_coil + 1j*w*L_coil
  197. voltage_in = Zeq / (Zsource + Zeq)
  198. self.circuitGain = np.abs(voltage_in)
  199. self.circuitPhase_deg = (180/np.pi)+np.angle(voltage_in)
  200. circuitImpedance_ohms = np.abs(Zsource + Zeq)
  201. ######################################################
  202. # PreAmp gain
  203. if self.nTransVersion == 4:
  204. self.PreAmpGain = 1000.
  205. elif self.nTransVersion == 1 or self.nTransVersion == 2 or self.nTransVersion == 3 or self.nTransVersion == 6:
  206. self.PreAmpGain = 500.
  207. else:
  208. print ("unsupported transmitter version")
  209. exit(1)
  210. # Total Receiver Gain
  211. self.RxGain = self.circuitGain * self.PreAmpGain
  212. #####################################################
  213. # Current gain
  214. if floor(self.nDAQVersion) == 1:
  215. self.CurrentGain = 150.
  216. elif floor(self.nDAQVersion) ==2:
  217. self.CurrentGain = 180.
  218. def updateProgress(self):
  219. pass
  220. def TDSmartStack(self, outlierTest, MADcutoff, canvas):
  221. Stack = {}
  222. # align for stacking and modulate
  223. for pulse in self.DATADICT["PULSES"]:
  224. stack = np.zeros(( len(self.DATADICT[pulse]["chan"]), self.DATADICT["nPulseMoments"],\
  225. len(self.DATADICT["stacks"]), len(self.DATADICT[pulse]["TIMES"]) ))
  226. for ipm in range(self.DATADICT["nPulseMoments"]):
  227. istack = 0
  228. for sstack in self.DATADICT["stacks"]:
  229. if self.pulseType == "FID" or pulse == "Pulse 2":
  230. mod = (-1)**(ipm%2) * (-1)**(sstack%2)
  231. elif self.pulseType == "4PhaseT1":
  232. mod = (-1)**(ipm%2) * (-1)**(((sstack-1)/2)%2)
  233. ichan = 0
  234. for chan in self.DATADICT[pulse]["chan"]:
  235. stack[ichan,ipm,istack,:] += mod*self.DATADICT[pulse][chan][ipm][sstack]
  236. ichan += 1
  237. istack += 1
  238. Stack[pulse] = stack
  239. #########################################
  240. # simple stack and plot of simple stack #
  241. #########################################
  242. canvas.reAxH2(np.shape(stack)[0], False, False)
  243. axes = canvas.fig.axes
  244. SimpleStack = {}
  245. VarStack = {}
  246. for pulse in self.DATADICT["PULSES"]:
  247. SimpleStack[pulse] = {}
  248. VarStack[pulse] = {}
  249. ichan = 0
  250. for chan in self.DATADICT[pulse]["chan"]:
  251. SimpleStack[pulse][chan] = 1e9*np.average( Stack[pulse][ichan], 1 )
  252. VarStack[pulse][chan] = 1e9*np.std( Stack[pulse][ichan], 1 )
  253. ax1 = axes[ 2*ichan ]
  254. #ax1.get_yaxis().get_major_formatter().set_useOffset(False)
  255. y_formatter = matplotlib.ticker.ScalarFormatter(useOffset=False)
  256. ax1.yaxis.set_major_formatter(y_formatter)
  257. ax1.plot( 1e3*self.DATADICT[pulse]["TIMES"], np.average( SimpleStack[pulse][chan], 0 ), color='darkblue' )
  258. ax1.set_title("Ch." + str(chan) + ": avg FID", fontsize=8)
  259. ax1.set_xlabel(r"time (ms)", fontsize=8)
  260. if ichan == 0:
  261. ax1.set_ylabel(r"signal [nV]", fontsize=8)
  262. else:
  263. plt.setp(ax1.get_yticklabels(), visible=False)
  264. plt.setp(ax1.get_yaxis().get_offset_text(), visible=False)
  265. # if ichan == 1:
  266. # canvas.ax2.plot( 1e3*self.DATADICT[pulse]["TIMES"], np.average( SimpleStack[pulse][chan], 0 ), color='darkblue' )
  267. # canvas.ax2.set_title("Ch." + str(chan) + ": total average FID", fontsize=8)
  268. # canvas.ax2.set_xlabel(r"time [ms]", fontsize=8)
  269. # if ichan == 2:
  270. # canvas.ax3.plot( 1e3*self.DATADICT[pulse]["TIMES"], np.average( SimpleStack[pulse][chan], 0 ), color='darkblue' )
  271. # canvas.ax3.set_title("Ch." + str(chan) + ": total average FID", fontsize=8)
  272. # canvas.ax3.set_xlabel(r"time [ms]", fontsize=8)
  273. # if ichan == 3:
  274. # canvas.ax4.plot( 1e3*self.DATADICT[pulse]["TIMES"], np.average( SimpleStack[pulse][chan], 0 ), color='darkblue' )
  275. # canvas.ax4.set_title("Ch." + str(chan) + ": total average FID", fontsize=8)
  276. # canvas.ax4.set_xlabel(r"time [ms]", fontsize=8)
  277. ichan += 1
  278. #########################
  279. # Oulier rejectig stack #
  280. #########################
  281. if outlierTest == "MAD":
  282. MADStack = {}
  283. VarStack = {}
  284. #1.4826 is assumption of gaussian noise
  285. madstack = np.zeros(( len(self.DATADICT[pulse]["chan"]),\
  286. self.DATADICT["nPulseMoments"], len(self.DATADICT[pulse]["TIMES"]) ))
  287. varstack = np.zeros(( len(self.DATADICT[pulse]["chan"]),\
  288. self.DATADICT["nPulseMoments"], len(self.DATADICT[pulse]["TIMES"]) ))
  289. for pulse in self.DATADICT["PULSES"]:
  290. MADStack[pulse] = {}
  291. VarStack[pulse] = {}
  292. ichan = 0
  293. for chan in self.DATADICT[pulse]["chan"]:
  294. ax1 = axes[ 2*ichan ]
  295. for ipm in range(self.DATADICT["nPulseMoments"]):
  296. # # brutal loop over time, can this be vectorized?
  297. # for it in range(len(self.DATADICT[pulse]["TIMES"])):
  298. # x = 1e9 *Stack[pulse][ichan,ipm,:,it]
  299. # MAD = 1.4826 * np.median( np.abs(x-np.median(x)) )
  300. # good = 0
  301. # for istack in self.DATADICT["stacks"]:
  302. # if (np.abs(x[istack-1]-np.median(x))) / MAD < 2:
  303. # good += 1
  304. # madstack[ ichan, ipm, it ] += x[istack-1]
  305. # else:
  306. # pass
  307. # madstack[ichan, ipm, it] /= good
  308. # percent = int(1e2* (float)(ipm) / (float)(self.DATADICT["nPulseMoments"]) )
  309. # self.progressTrigger.emit(percent)
  310. # Vectorized version of above...much, much faster
  311. x = 1e9*copy.deepcopy(Stack[pulse][ichan][ipm,:,:]) # stack and time indices
  312. tile_med = np.tile( np.median(x, axis=0), (np.shape(x)[0],1))
  313. MAD = MADcutoff * np.median(np.abs(x - tile_med), axis=0)
  314. tile_MAD = np.tile( MAD, (np.shape(x)[0],1))
  315. good = np.abs(x-tile_med)/tile_MAD < 2. # 1.4826 # 2
  316. madstack[ichan][ipm] = copy.deepcopy( np.ma.masked_array(x, good != True).mean(axis=0) )
  317. varstack[ichan][ipm] = copy.deepcopy( np.ma.masked_array(x, good != True).std(axis=0) )
  318. # reporting
  319. percent = int(1e2* (float)((ipm)+ichan*self.DATADICT["nPulseMoments"]) /
  320. (float)(self.DATADICT["nPulseMoments"] * len(self.DATADICT[pulse]["chan"])))
  321. self.progressTrigger.emit(percent)
  322. ax1.plot( 1e3*self.DATADICT[pulse]["TIMES"], np.average( madstack[ichan], 0 ) , color='darkred')
  323. MADStack[pulse][chan] = madstack[ichan]
  324. VarStack[pulse][chan] = varstack[ichan]
  325. ichan += 1
  326. self.DATADICT["stack"] = MADStack
  327. else:
  328. self.DATADICT["stack"] = SimpleStack
  329. #########################################
  330. # Plot Fourier Transform representation #
  331. #########################################
  332. # canvas.fig.subplots_adjust(right=0.8)
  333. # cbar_ax = canvas.fig.add_axes([0.85, 0.1, 0.015, 0.355])
  334. # cbar_ax.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  335. im2 = []
  336. for pulse in self.DATADICT["PULSES"]:
  337. ichan = 0
  338. axes = canvas.fig.axes
  339. vvmin = 1e10
  340. vvmax = 0
  341. for chan in self.DATADICT[pulse]["chan"]:
  342. ax1 = axes[2*ichan ]
  343. ax2 = axes[2*ichan+1] # TODO fix hard coded number
  344. if outlierTest == "MAD":
  345. X = np.fft.rfft( MADStack[pulse][chan][0,:] )
  346. nu = np.fft.fftfreq(len( MADStack[pulse][chan][0,:]), d=self.dt)
  347. else:
  348. X = np.fft.rfft( SimpleStack[pulse][chan][0,:] )
  349. nu = np.fft.fftfreq(len( SimpleStack[pulse][chan][0,:]), d=self.dt)
  350. nu = nu[0:len(X)]
  351. nu[-1] = np.abs(nu[-1])
  352. df = nu[1] - nu[0]
  353. of = 0
  354. istart = int((self.transFreq-150.)/df)
  355. iend = int((self.transFreq+150.)/df)
  356. of = nu[istart]
  357. def freqlabel(xxx, pos):
  358. return '%1.0f' %(of + xxx*df)
  359. formatter = FuncFormatter(freqlabel)
  360. SFFT = np.zeros( (self.DATADICT["nPulseMoments"], len(X)), dtype=np.complex64 )
  361. SFFT[0,:] = X
  362. for ipm in range(1, self.DATADICT["nPulseMoments"]):
  363. if outlierTest == "MAD":
  364. SFFT[ipm,:] = np.fft.rfft( MADStack[pulse][chan][ipm,:] )
  365. else:
  366. SFFT[ipm,:] = np.fft.rfft( SimpleStack[pulse][chan][ipm,:] )
  367. # convert to dB and add colorbars
  368. #db = 20.*np.log10(np.abs(SFFT[:,istart:iend]))
  369. db = (np.abs(SFFT[:,istart:iend]))
  370. vvmin = min(vvmin, np.min (db))
  371. vvmax = max(vvmax, np.max (db))
  372. im2.append(ax2.matshow( db, aspect='auto', cmap=cmocean.cm.ice, vmin=vvmin, vmax=vvmax))
  373. #im2 = ax2.matshow( db, aspect='auto', cmap=cmocean.cm.ice, vmin=vvmin, vmax=vvmax)
  374. if ichan == 0:
  375. ax2.set_ylabel(r"$q$ (A $\cdot$ s)", fontsize=8)
  376. else:
  377. #ax2.yaxis.set_ticklabels([])
  378. plt.setp(ax2.get_yticklabels(), visible=False)
  379. ax2.xaxis.set_major_formatter(formatter)
  380. ax2.xaxis.set_ticks_position('bottom')
  381. ax2.xaxis.set_major_locator(MaxNLocator(3))
  382. y_formatter = matplotlib.ticker.ScalarFormatter(useOffset=False)
  383. ax2.yaxis.set_major_formatter(y_formatter)
  384. #if chan == self.DATADICT[pulse]["chan"][-1]:
  385. #cb2 = canvas.fig.colorbar(im2, cax=cbar_ax, format='%1.0e')
  386. #cb2 = canvas.fig.colorbar(im2[0], ax=ax2, format='%1.0e', orientation='horizontal')
  387. #cb2 = canvas.fig.colorbar(im2, ax=ax2, format='%1.0e', orientation='horizontal')
  388. #cb2.ax.tick_params(axis='both', which='major', labelsize=8)
  389. #cb2.set_label("signal (dB)", fontsize=8)
  390. ichan += 1
  391. canvas.fig.subplots_adjust(hspace=.1, wspace=.05, left=.075, right=.95 )#left=None, bottom=None, right=None, top=None, wspace=None, hspace=None)
  392. #cb1 = canvas.fig.colorbar(im, ax=axes[0::2], format='%1.0e', orientation='horizontal', shrink=.35, aspect=30)
  393. #cb1.ax.tick_params(axis='both', which='major', labelsize=8)
  394. #cb1.set_label("$\mathcal{V}_N$ (nV)", fontsize=8)
  395. cb2 = canvas.fig.colorbar(im2[-1], ax=axes[1::2], format='%1.0e', orientation='horizontal', shrink=.35, aspect=30)
  396. cb2.ax.tick_params(axis='both', which='major', labelsize=8)
  397. cb2.set_label("$\mathcal{V}_N$ (nV)", fontsize=8)
  398. #canvas.fig.tight_layout()
  399. canvas.draw()
  400. self.doneTrigger.emit()
  401. def sumData(self, canvas, fred):
  402. chans = copy.deepcopy(self.DATADICT[self.DATADICT["PULSES"][0]]["chan"]) #= np.array( ( self.DATADICT[pulse]["chan"][0], ) )
  403. nchan = len(chans)
  404. # Sum permutations of two channel combos
  405. for ich in range(nchan-1):
  406. for ch in chans[ich+1:]:
  407. chsum = chans[ich] + "+" + ch
  408. for pulse in self.DATADICT["PULSES"]:
  409. self.DATADICT[pulse][chsum] = {}
  410. for ipm in range(self.DATADICT["nPulseMoments"]):
  411. self.DATADICT[pulse][chsum][ipm] = {}
  412. for istack in self.DATADICT["stacks"]:
  413. self.DATADICT[pulse][chsum][ipm][istack] = self.DATADICT[pulse][chans[ich]][ipm][istack] + self.DATADICT[pulse][ch][ipm][istack]
  414. self.DATADICT[pulse]["chan"].append(chsum)
  415. # Sum all channels
  416. chsum = ""
  417. for ch in chans:
  418. chsum += ch + "+"
  419. chsum = chsum[0:-1] # remove last "+"
  420. for pulse in self.DATADICT["PULSES"]:
  421. self.DATADICT[pulse][chsum] = {}
  422. for ipm in range(self.DATADICT["nPulseMoments"]):
  423. self.DATADICT[pulse][chsum][ipm] = {}
  424. for istack in self.DATADICT["stacks"]:
  425. self.DATADICT[pulse][chsum][ipm][istack] = copy.deepcopy(self.DATADICT[pulse][chans[0]][ipm][istack])
  426. for ch in chans[1:]:
  427. self.DATADICT[pulse][chsum][ipm][istack] += self.DATADICT[pulse][ch][ipm][istack]
  428. self.DATADICT[pulse]["chan"].append(chsum)
  429. # if nchan > 2:
  430. # for ch in chans:
  431. # chsum += ch
  432. # for ch2 in chans[1::]:
  433. # for pulse in self.DATADICT["PULSES"]:
  434. # self.DATADICT[pulse][chsum] = {}
  435. # for istack in self.DATADICT["stacks"]:
  436. # self.DATADICT[pulse][chsum][ipm][istack] = self.DATADICT[pulse][chans[ich]][ipm][istack] + self.DATADICT[pulse][ch][ipm][istack]
  437. self.doneTrigger.emit()
  438. def quadDet(self, clip, phase, canvas):
  439. from scipy import signal
  440. self.RotatedAmplitude = True
  441. wL = self.transFreq * 2*np.pi
  442. vL = self.transFreq
  443. #T = 50
  444. dt = self.dt
  445. #DT = 0.01
  446. CA = {} # corrected amplitude
  447. IP = {} # instantaneous phase
  448. NR = {} # Noise residual
  449. RE = {} # Real channel
  450. IM = {} # Imaginary channel
  451. # global maximums for plotting
  452. CAmax = {}
  453. NRmax = {}
  454. REmax = {}
  455. IMmax = {}
  456. first = False
  457. self.sigma = {}
  458. for pulse in self.DATADICT["PULSES"]:
  459. CA[pulse] = {}
  460. IP[pulse] = {}
  461. NR[pulse] = {}
  462. RE[pulse] = {}
  463. IM[pulse] = {}
  464. CAmax[pulse] = 0
  465. NRmax[pulse] = 0
  466. REmax[pulse] = 0
  467. IMmax[pulse] = 0
  468. ichan = 0
  469. self.sigma[pulse] = {}
  470. for chan in self.DATADICT[pulse]["chan"]:
  471. CA[pulse][chan] = np.zeros( (self.DATADICT["nPulseMoments"], len(self.DATADICT[pulse]["TIMES"])-clip ) )
  472. IP[pulse][chan] = np.zeros( (self.DATADICT["nPulseMoments"], len(self.DATADICT[pulse]["TIMES"])-clip ) )
  473. NR[pulse][chan] = np.zeros( (self.DATADICT["nPulseMoments"], len(self.DATADICT[pulse]["TIMES"])-clip ) )
  474. RE[pulse][chan] = np.zeros( (self.DATADICT["nPulseMoments"], len(self.DATADICT[pulse]["TIMES"])-clip ) )
  475. IM[pulse][chan] = np.zeros( (self.DATADICT["nPulseMoments"], len(self.DATADICT[pulse]["TIMES"])-clip ) )
  476. for ipm in range(0, self.DATADICT["nPulseMoments"]):
  477. #t = self.DATADICT[pulse]["TIMES"] - self.DATADICT[pulse]["PULSE_TIMES"][-1]
  478. xn = self.DATADICT["stack"][pulse][chan][ipm,:]
  479. ht = signal.hilbert(xn)*np.exp(-1j*wL*self.DATADICT[pulse]["TIMES"])
  480. #############################################################
  481. # Quadrature signal
  482. RE[pulse][chan][ipm,:] = np.real(ht[clip::])
  483. IM[pulse][chan][ipm,:] = np.imag(ht[clip::])
  484. REmax[pulse] = max(REmax[pulse], np.max(np.real(ht[clip::])))
  485. IMmax[pulse] = max(IMmax[pulse], np.max(np.imag(ht[clip::])))
  486. #############################################################
  487. # Instantaneous phase
  488. IP[pulse][chan][ipm,:] = np.angle(ht)[clip::]
  489. #############################################################
  490. # Rotated amplitude
  491. [E0, df, phi, T2] = decay.quadratureDetect( ht.real, ht.imag, self.DATADICT[pulse]["TIMES"] )
  492. D = self.RotateAmplitude( ht.real, ht.imag, phi, df, self.DATADICT[pulse]["TIMES"] )
  493. CA[pulse][chan][ipm,:] = D.imag[clip::] # amplitude data
  494. NR[pulse][chan][ipm,:] = D.real[clip::] # noise data
  495. CAmax[pulse] = max(CAmax[pulse], np.max(D.imag[clip::]) )
  496. NRmax[pulse] = max(NRmax[pulse], np.max(D.real[clip::]) )
  497. self.sigma[pulse][chan] = np.std(NR[pulse][chan])
  498. # reporting
  499. percent = int(1e2* (float)((ipm)+ichan*self.DATADICT["nPulseMoments"]) /
  500. (float)(self.DATADICT["nPulseMoments"] * len(self.DATADICT[pulse]["chan"])))
  501. self.progressTrigger.emit(percent)
  502. ichan += 1
  503. self.DATADICT["CA"] = CA
  504. self.DATADICT["IP"] = IP
  505. self.DATADICT["NR"] = NR
  506. self.DATADICT["RE"] = RE
  507. self.DATADICT["IM"] = IM
  508. self.DATADICT["CAmax"] = CAmax
  509. self.DATADICT["NRmax"] = NRmax
  510. self.DATADICT["REmax"] = REmax
  511. self.DATADICT["IMmax"] = IMmax
  512. self.doneTrigger.emit()
  513. def plotQuadDet(self, clip, phase, canvas):
  514. canvas.reAxH2( len(self.DATADICT[ self.DATADICT["PULSES"][0] ]["chan"] ), False, False)
  515. ###############
  516. # Plot on GUI #
  517. ###############
  518. dcmap = cmocean.cm.balance_r #"RdBu" #YlGn" # "coolwarm_r" # diverging
  519. canvas.reAxH2( len(self.DATADICT[ self.DATADICT["PULSES"][0] ]["chan"] ), False, False)
  520. for pulse in self.DATADICT["PULSES"]:
  521. ichan = 0
  522. axes = canvas.fig.axes
  523. mmaxr = 0.
  524. mmaxi = 0.
  525. if clip > 0:
  526. time_sp = 1e3 * (self.DATADICT[pulse]["TIMES"][clip-1::] - self.DATADICT[pulse]["PULSE_TIMES"][-1] )
  527. else:
  528. time_sp = 1e3 * (self.DATADICT[pulse]["TIMES"] - self.DATADICT[pulse]["PULSE_TIMES"][-1] )
  529. QQ = np.average(self.DATADICT[pulse]["Q"], axis=1 )
  530. for chan in self.DATADICT[pulse]["chan"]:
  531. ax1 = axes[2*ichan ]
  532. ax2 = axes[2*ichan+1] # TODO fix hard coded number
  533. if phase == 0: # Re Im
  534. im1 = ax1.pcolormesh( time_sp, QQ, self.DATADICT["RE"][pulse][chan], cmap=dcmap, rasterized=True,\
  535. vmin=-self.DATADICT["REmax"][pulse] , vmax=self.DATADICT["REmax"][pulse] )
  536. im2 = ax2.pcolormesh( time_sp, QQ, self.DATADICT["IM"][pulse][chan], cmap=dcmap, rasterized=True,\
  537. vmin=-self.DATADICT["IMmax"][pulse], vmax=self.DATADICT["IMmax"][pulse] )
  538. if phase == 1: # Amp phase
  539. im1 = ax1.pcolormesh( time_sp, QQ, self.DATADICT["CA"][pulse][chan], cmap=dcmap, rasterized=True,
  540. vmin=-self.DATADICT["CAmax"][pulse] , vmax=self.DATADICT["CAmax"][pulse] )
  541. im2 = ax2.pcolormesh( time_sp, QQ, self.DATADICT["IP"][pulse][chan], cmap=cmocean.cm.phase, rasterized=True,\
  542. vmin=-np.pi, vmax=np.pi)
  543. if phase == 2: # CA NR
  544. im1 = ax1.pcolormesh( time_sp, QQ, self.DATADICT["CA"][pulse][chan], cmap=dcmap, rasterized=True,\
  545. vmin=-self.DATADICT["CAmax"][pulse] , vmax=self.DATADICT["CAmax"][pulse] )
  546. im2 = ax2.pcolormesh( time_sp, QQ, self.DATADICT["NR"][pulse][chan], cmap=dcmap, rasterized=True,
  547. vmin=-self.DATADICT["NRmax"][pulse] , vmax=self.DATADICT["NRmax"][pulse] )
  548. # cb2 = canvas.fig.colorbar(im2, ax=ax2, format='%1.0e')
  549. # cb2.set_label("Noise residual (nV)", fontsize=8)
  550. # cb2.ax.tick_params(axis='both', which='major', labelsize=8)
  551. # cb1 = canvas.fig.colorbar(im1, ax=ax1, format='%1.0e')
  552. # cb1.set_label("Phased amplitude (nV)", fontsize=8)
  553. # cb1.ax.tick_params(axis='both', which='major', labelsize=8)
  554. # cb2 = canvas.fig.colorbar(im2, ax=ax2, format="%1.0e")
  555. # cb2.set_label("Phase (rad)", fontsize=8)
  556. # cb2.ax.tick_params(axis='both', which='major', labelsize=8)
  557. # cb1 = canvas.fig.colorbar(im1, ax=ax1, format="%1.0e")
  558. # cb1.set_label("FID (nV)", fontsize=8)
  559. # cb1.ax.tick_params(axis='both', which='major', labelsize=8)
  560. # if you save these as pdf or eps, there are artefacts
  561. # for cbar in [cb1,cb2]:
  562. # #cbar.solids.set_rasterized(True)
  563. # cbar.solids.set_edgecolor("face")
  564. # reporting
  565. percent = int(1e2* (float)(ichan)/len(self.DATADICT[pulse]["chan"]))
  566. self.progressTrigger.emit(percent)
  567. if ichan == 0:
  568. ax1.set_ylabel(r"$q$ ( $\mathrm{A}\cdot\mathrm{s}$)", fontsize=8)
  569. ax2.set_ylabel(r"$q$ ( $\mathrm{A}\cdot\mathrm{s}$)", fontsize=8)
  570. else:
  571. #ax2.yaxis.set_ticklabels([])
  572. #ax1.yaxis.set_ticklabels([])
  573. plt.setp(ax1.get_yticklabels(), visible=False)
  574. plt.setp(ax2.get_yticklabels(), visible=False)
  575. ichan += 1
  576. ax1.set_yscale('log')
  577. ax2.set_yscale('log')
  578. plt.setp(ax1.get_xticklabels(), visible=False)
  579. ax1.set_ylim( np.min(QQ), np.max(QQ) )
  580. ax2.set_ylim( np.min(QQ), np.max(QQ) )
  581. ax1.set_xlim( np.min(time_sp), np.max(time_sp) )
  582. ax2.set_xlim( np.min(time_sp), np.max(time_sp) )
  583. #ax2.set_xlabel(r"Time since end of pulse (ms)", fontsize=8)
  584. ax2.set_xlabel(r"Time (ms)", fontsize=8)
  585. canvas.fig.subplots_adjust(hspace=.15, wspace=.05, left=.075, right=.95, bottom=.1, top=.95 )#left=None, bottom=None, right=None, top=None, wspace=None, hspace=None)
  586. tick_locator = MaxNLocator(nbins=5)
  587. cb1 = canvas.fig.colorbar(im1, ax=axes[0::2], format='%1.0e', orientation='horizontal', shrink=.35, aspect=30)
  588. cb1.ax.tick_params(axis='both', which='major', labelsize=8)
  589. cb1.set_label("$\mathcal{V}_N$ (nV)", fontsize=8)
  590. cb1.locator = tick_locator
  591. cb1.update_ticks()
  592. cb2 = canvas.fig.colorbar(im2, ax=axes[1::2], format='%1.0e', orientation='horizontal', shrink=.35, aspect=30, pad=.2)
  593. cb2.ax.tick_params(axis='both', which='major', labelsize=8)
  594. cb2.set_label("$\mathcal{V}_N$ (nV)", fontsize=8)
  595. cb2.locator = tick_locator
  596. cb2.update_ticks()
  597. canvas.draw()
  598. self.doneTrigger.emit()
  599. def RotateAmplitude(self, X, Y, zeta, df, t):
  600. V = X + 1j*Y
  601. return np.abs(V) * np.exp( 1j * ( np.angle(V) - zeta - 2.*np.pi*df*t ) )
  602. def gateIntegrate( self, gpd, clip, canvas ):
  603. """ Gate integrate the real, imaginary, phased, and noise residual channels
  604. """
  605. self.gated = True
  606. self.GATED = {}
  607. for pulse in self.DATADICT["PULSES"]:
  608. QQ = np.average(self.DATADICT[pulse]["Q"], axis=1 )
  609. ichan = 0
  610. for chan in self.DATADICT[pulse]["chan"]:
  611. self.GATED[chan] = {}
  612. for ipm in range(0, self.DATADICT["nPulseMoments"]):
  613. # Time since pulse rather than since record starts...
  614. if clip > 0:
  615. time_sp = 1e3 * (self.DATADICT[pulse]["TIMES"][clip-1::] - self.DATADICT[pulse]["PULSE_TIMES"][-1] )
  616. else:
  617. time_sp = 1e3 * (self.DATADICT[pulse]["TIMES"] - self.DATADICT[pulse]["PULSE_TIMES"][-1] )
  618. #GT, GD, GTT, sig_stack, isum = rotate.gateIntegrate( self.DATADICT["CA"][pulse][chan][ipm,:], time_sp, gpd, self.sigma[pulse][chan], 1.5 )
  619. #GT2, GP, GTT, sig_stack_err, isum = rotate.gateIntegrate( self.DATADICT["NR"][pulse][chan][ipm,:], time_sp, gpd, self.sigma[pulse][chan], 1.5 )
  620. GT, GCA, GTT, sig_stack, isum = rotate.gateIntegrate( self.DATADICT["CA"][pulse][chan][ipm,:], time_sp, gpd, self.sigma[pulse][chan], 1.5 )
  621. GT, GNR, GTT, sig_stack, isum = rotate.gateIntegrate( self.DATADICT["NR"][pulse][chan][ipm,:], time_sp, gpd, self.sigma[pulse][chan], 1.5 )
  622. GT, GRE, GTT, sig_stack, isum = rotate.gateIntegrate( self.DATADICT["RE"][pulse][chan][ipm,:], time_sp, gpd, self.sigma[pulse][chan], 1.5 )
  623. GT, GIM, GTT, sig_stack, isum = rotate.gateIntegrate( self.DATADICT["IM"][pulse][chan][ipm,:], time_sp, gpd, self.sigma[pulse][chan], 1.5 )
  624. GT, GIP, GTT, sig_stack, isum = rotate.gateIntegrate( self.DATADICT["IP"][pulse][chan][ipm,:], time_sp, gpd, self.sigma[pulse][chan], 1.5 )
  625. if ipm == 0:
  626. # self.GATED[chan]["DATA"] = np.zeros( ( self.DATADICT["nPulseMoments"], len(GT)) )
  627. # self.GATED[chan]["ERR"] = np.zeros( ( self.DATADICT["nPulseMoments"], len(GT)) )
  628. # self.GATED[chan]["SIGMA"] = np.zeros( ( self.DATADICT["nPulseMoments"], len(GT)) )
  629. self.GATED[chan]["CA"] = np.zeros( ( self.DATADICT["nPulseMoments"], len(GT)) )
  630. self.GATED[chan]["NR"] = np.zeros( ( self.DATADICT["nPulseMoments"], len(GT)) )
  631. self.GATED[chan]["RE"] = np.zeros( ( self.DATADICT["nPulseMoments"], len(GT)) )
  632. self.GATED[chan]["IM"] = np.zeros( ( self.DATADICT["nPulseMoments"], len(GT)) )
  633. self.GATED[chan]["IP"] = np.zeros( ( self.DATADICT["nPulseMoments"], len(GT)) )
  634. self.GATED[chan]["isum"] = isum
  635. #self.GATED[chan]["DATA"][ipm] = GD.real
  636. self.GATEDABSCISSA = GT
  637. self.GATEDWINDOW = GTT
  638. #self.GATED[chan]["SIGMA"][ipm] = sig_stack #_err # GP.real
  639. #self.GATED[chan]["ERR"][ipm] = GP.real
  640. self.GATED[chan]["CA"][ipm] = GCA.real
  641. self.GATED[chan]["NR"][ipm] = GNR.real
  642. self.GATED[chan]["RE"][ipm] = GRE.real
  643. self.GATED[chan]["IM"][ipm] = GIM.real
  644. self.GATED[chan]["IP"][ipm] = GIP.real
  645. percent = int(1e2* (float)((ipm)+ichan*self.DATADICT["nPulseMoments"]) /
  646. (float)(self.DATADICT["nPulseMoments"] * len(self.DATADICT[pulse]["chan"])))
  647. self.progressTrigger.emit(percent)
  648. self.GATED[chan]["GTT"] = GTT
  649. self.GATED[chan]["GT"] = GT
  650. self.GATED[chan]["QQ"] = QQ
  651. ichan += 1
  652. self.doneTrigger.emit()
  653. def bootstrap_resample(self, X, n=None):
  654. # from http://nbviewer.jupyter.org/gist/aflaxman/6871948
  655. """ Bootstrap resample an array_like
  656. Parameters
  657. ----------
  658. X : array_like
  659. data to resample
  660. n : int, optional
  661. length of resampled array, equal to len(X) if n==None
  662. Results
  663. -------
  664. returns X_resamples
  665. """
  666. if n == None:
  667. n = len(X)
  668. resample_i = np.floor(np.random.rand(n)*len(X)).astype(int)
  669. return X[resample_i]
  670. def bootstrap_sigma(self, pulse, chan):
  671. # bootstrap resample
  672. nt = len(self.GATED[chan]["GT"])
  673. nb = 5000
  674. XS = np.zeros( (nb, nt) )
  675. for ii in range(nb):
  676. for it in range(nt):
  677. if self.GATED[chan]["isum"][it] < 8:
  678. XS[ii, it] = self.sigma[pulse][chan]
  679. else:
  680. if it == 0:
  681. X = self.bootstrap_resample( np.concatenate( (self.GATED[chan]["NR"][:,it], self.GATED[chan]["NR"][:,it+1], \
  682. self.GATED[chan]["NR"][:,it+2], self.GATED[chan]["NR"][:,it+3] ) ), n=nt )
  683. elif it == 1:
  684. X = self.bootstrap_resample( np.concatenate( (self.GATED[chan]["NR"][:,it-1], self.GATED[chan]["NR"][:,it], \
  685. self.GATED[chan]["NR"][:,it+1], self.GATED[chan]["NR"][:,it+2] ) ), n=nt )
  686. elif it == nt-2:
  687. X = self.bootstrap_resample( np.concatenate( (self.GATED[chan]["NR"][:,it+1], self.GATED[chan]["NR"][:,it], \
  688. self.GATED[chan]["NR"][:,it-1], self.GATED[chan]["NR"][:,it-2] ) ), n=nt )
  689. elif it == nt-1:
  690. X = self.bootstrap_resample( np.concatenate( (self.GATED[chan]["NR"][:,it], self.GATED[chan]["NR"][:,it-1], \
  691. self.GATED[chan]["NR"][:,it-2], self.GATED[chan]["NR"][:,it-3] ) ), n=nt )
  692. else:
  693. X = self.bootstrap_resample( np.concatenate( (self.GATED[chan]["NR"][:,it-2] , self.GATED[chan]["NR"][:,it-1], \
  694. self.GATED[chan]["NR"][:,it], self.GATED[chan]["NR"][:,it+1], self.GATED[chan]["NR"][:,it+2] )), n=nt )
  695. XS[ii,it] = np.std(X)
  696. return XS
  697. def plotGateIntegrate( self, gpd, clip, phase, canvas ):
  698. """ Plot the gate integration
  699. """
  700. canvas.reAxH2( len(self.DATADICT[ self.DATADICT["PULSES"][0] ]["chan"] ), False, False)
  701. axes = canvas.fig.axes
  702. cmap = cmocean.cm.balance_r
  703. # Calculate maximum for plotting...TODO move into loop above
  704. vmax1 = 0
  705. vmax2 = 0
  706. for pulse in self.DATADICT["PULSES"]:
  707. for chan in self.DATADICT[pulse]["chan"]:
  708. if phase == 0:
  709. vmax1 = max(vmax1, np.max(np.abs(self.GATED[chan]["RE"])))
  710. vmax2 = max(vmax2, np.max(np.abs(self.GATED[chan]["IM"])))
  711. elif phase == 1:
  712. vmax1 = max(vmax1, np.max(np.abs(self.GATED[chan]["CA"])))
  713. vmax2 = np.pi
  714. elif phase == 2:
  715. vmax1 = max(vmax1, np.max(np.abs(self.GATED[chan]["CA"])))
  716. vmax2 = max(vmax2, np.max(np.abs(self.GATED[chan]["NR"])))
  717. for pulse in self.DATADICT["PULSES"]:
  718. ichan = 0
  719. for chan in self.DATADICT[pulse]["chan"]:
  720. ax1 = axes[2*ichan ]
  721. ax2 = axes[2*ichan+1]
  722. if phase == 0:
  723. im1 = ax1.pcolormesh(self.GATED[chan]["GTT"], self.GATED[chan]["QQ"], self.GATED[chan]["RE"], cmap=cmap, vmin=-vmax1, vmax=vmax1)
  724. im2 = ax2.pcolormesh(self.GATED[chan]["GTT"], self.GATED[chan]["QQ"], self.GATED[chan]["IM"], cmap=cmap, vmin=-vmax2, vmax=vmax2)
  725. elif phase == 1:
  726. im1 = ax1.pcolormesh(self.GATED[chan]["GTT"], self.GATED[chan]["QQ"], self.GATED[chan]["CA"], cmap=cmap, vmin=-vmax1, vmax=vmax1)
  727. im2 = ax2.pcolormesh(self.GATED[chan]["GTT"], self.GATED[chan]["QQ"], self.GATED[chan]["IP"], cmap=cmocean.cm.phase, vmin=-vmax2, vmax=vmax2)
  728. elif phase == 2:
  729. im1 = ax1.pcolormesh(self.GATED[chan]["GTT"], self.GATED[chan]["QQ"], self.GATED[chan]["CA"], cmap=cmap, vmin=-vmax1, vmax=vmax1)
  730. XS = self.bootstrap_sigma(pulse, chan)
  731. #im2 = ax2.pcolormesh(self.GATED[chan]["GTT"], self.GATED[chan]["QQ"], self.GATED[chan]["NR"], cmap=cmap, vmin=-vmax2, vmax=vmax2)
  732. # bootstrap resample
  733. # nt = len(self.GATED[chan]["GT"])
  734. # nb = 5000
  735. # XS = np.zeros( (nb, nt) )
  736. # for ii in range(nb):
  737. # #XS = []
  738. # for it in range(nt):
  739. # if self.GATED[chan]["isum"][it] < 8:
  740. # XS[ii, it] = self.sigma[pulse][chan]
  741. # else:
  742. # if it == 0:
  743. # X = self.bootstrap_resample( np.concatenate( (self.GATED[chan]["NR"][:,it], self.GATED[chan]["NR"][:,it+1], \
  744. # self.GATED[chan]["NR"][:,it+2], self.GATED[chan]["NR"][:,it+3] ) ), n=nt )
  745. # if it == 1:
  746. # X = self.bootstrap_resample( np.concatenate( (self.GATED[chan]["NR"][:,it-1], self.GATED[chan]["NR"][:,it], \
  747. # self.GATED[chan]["NR"][:,it+1], self.GATED[chan]["NR"][:,it+2] ) ), n=nt )
  748. # elif it == nt-2:
  749. # X = self.bootstrap_resample( np.concatenate( (self.GATED[chan]["NR"][:,it+1], self.GATED[chan]["NR"][:,it], \
  750. # self.GATED[chan]["NR"][:,it-1], self.GATED[chan]["NR"][:,it-2] ) ), n=nt )
  751. # elif it == nt-1:
  752. # X = self.bootstrap_resample( np.concatenate( (self.GATED[chan]["NR"][:,it], self.GATED[chan]["NR"][:,it-1], \
  753. # self.GATED[chan]["NR"][:,it-2], self.GATED[chan]["NR"][:,it-3] ) ), n=nt )
  754. # else:
  755. # X = self.bootstrap_resample( np.concatenate( (self.GATED[chan]["NR"][:,it-2] , self.GATED[chan]["NR"][:,it-1], \
  756. # self.GATED[chan]["NR"][:,it], self.GATED[chan]["NR"][:,it+1], self.GATED[chan]["NR"][:,it+2] )), n=nt )
  757. # XS[ii,it] = np.std(X)
  758. #if ii == 0:
  759. # ax2.plot( self.GATED[chan]["GT"], XS[ii], '-', linewidth=1, markersize=4, alpha=.5, color='lightgrey', label = "bootstrap sim" )
  760. #else:
  761. # ax2.plot( self.GATED[chan]["GT"], XS[ii], '-', linewidth=1, markersize=4, alpha=.5, color='lightgrey' )
  762. ax2.plot( self.GATED[chan]["GT"], np.std(self.GATED[chan]["NR"], axis=0), color='darkgrey', linewidth=2, label="gate std" )
  763. ax2.plot( np.tile(self.GATED[chan]["GT"], (5000,1) ), XS, '.', color='lightgrey', linewidth=1, alpha=.5 )
  764. ax2.plot( self.GATED[chan]["GT"], np.average(XS, axis=0), color='black', linewidth=2, label="bootstrap avg." )
  765. ax2.plot( self.GATED[chan]["GT"], np.median(XS, axis=0), color='black', linewidth=2, label="bootstrap med." )
  766. ax2.legend()
  767. im1.set_edgecolor('face')
  768. if phase != 2:
  769. im2.set_edgecolor('face')
  770. plt.setp(ax1.get_xticklabels(), visible=False)
  771. ax1.set_ylim( np.min(self.GATED[chan]["QQ"]), np.max(self.GATED[chan]["QQ"]) )
  772. if phase != 2:
  773. ax2.set_ylim( np.min(self.GATED[chan]["QQ"]), np.max(self.GATED[chan]["QQ"]) )
  774. ax1.set_xlim( np.min(self.GATED[chan]["GTT"]), np.max(self.GATED[chan]["GTT"]) )
  775. ax2.set_xlim( np.min(self.GATED[chan]["GTT"]), np.max(self.GATED[chan]["GTT"]) )
  776. ax1.set_yscale('log')
  777. ax2.set_yscale('log')
  778. qlabs = np.append(np.concatenate( ( self.GATED[chan]["QQ"][0:1], self.GATED[chan]["QQ"][9::10] )), self.GATED[chan]["QQ"][-1] )
  779. ax1.yaxis.set_ticks( qlabs ) # np.append(np.concatenate( (QQ[0:1],QQ[9::10] )), QQ[-1] ) )
  780. if phase != 2:
  781. ax2.yaxis.set_ticks( qlabs ) #np.append(np.concatenate( (QQ[0:1],QQ[9::10] )), QQ[-1] ) )
  782. #formatter = matplotlib.ticker.LogFormatter(10, labelOnlyBase=False)
  783. formatter = matplotlib.ticker.FuncFormatter(lambda x, pos: str((round(x,1))))
  784. ax1.yaxis.set_major_formatter(formatter) #matplotlib.ticker.FormatStrFormatter('%d.1'))
  785. ax2.yaxis.set_major_formatter(formatter) #matplotlib.ticker.FormatStrFormatter('%d.1'))
  786. ax1.xaxis.set_major_formatter(formatter) #matplotlib.ticker.FormatStrFormatter('%d.1'))
  787. ax2.xaxis.set_major_formatter(formatter) #matplotlib.ticker.FormatStrFormatter('%d.1'))
  788. ax1.set_xscale('log')
  789. ax2.set_xscale('log')
  790. if ichan == 0:
  791. ax1.set_ylabel(r"$q$ ( $\mathrm{A}\cdot\mathrm{s}$)", fontsize=8)
  792. if phase == 2:
  793. ax2.set_ylabel(r"noise est. (nV)", fontsize=8)
  794. else:
  795. ax2.set_ylabel(r"$q$ ( $\mathrm{A}\cdot\mathrm{s}$)", fontsize=8)
  796. else:
  797. plt.setp(ax1.get_yticklabels(), visible=False)
  798. plt.setp(ax2.get_yticklabels(), visible=False)
  799. ax2.set_xlabel(r"$t-\tau_p$ (ms)", fontsize=8)
  800. ichan += 1
  801. #canvas.fig.subplots_adjust(hspace=.1, wspace=.05, left=.075, right=.925 )#left=None, bottom=None, right=None, top=None, wspace=None, hspace=None)
  802. #canvas.fig.tight_layout()
  803. #canvas.draw()
  804. canvas.fig.subplots_adjust(hspace=.15, wspace=.05, left=.075, right=.95, bottom=.1, top=.95 )#left=None, bottom=None, right=None, top=None, wspace=None, hspace=None)
  805. tick_locator = MaxNLocator(nbins=5)
  806. cb1 = canvas.fig.colorbar(im1, ax=axes[0::2], format='%1.0e', orientation='horizontal', shrink=.35, aspect=30)
  807. cb1.ax.tick_params(axis='both', which='major', labelsize=8)
  808. cb1.set_label("$\mathcal{V}_N$ (nV)", fontsize=8)
  809. cb1.locator = tick_locator
  810. cb1.update_ticks()
  811. if phase != 2:
  812. cb2 = canvas.fig.colorbar(im2, ax=axes[1::2], format='%1.0e', orientation='horizontal', shrink=.35, aspect=30, pad=.2)
  813. cb2.ax.tick_params(axis='both', which='major', labelsize=8)
  814. cb2.set_label("$\mathcal{V}_N$ (nV)", fontsize=8)
  815. cb2.locator = tick_locator
  816. cb2.update_ticks()
  817. canvas.draw()
  818. self.doneTrigger.emit()
  819. def FDSmartStack(self, cv, canvas):
  820. from matplotlib.colors import LogNorm
  821. from matplotlib.ticker import MaxNLocator
  822. """
  823. Currently this stacks 4-phase second pulse data only, we need to generalise
  824. """
  825. try:
  826. canvas.fig.clear()
  827. except:
  828. pass
  829. self.dataCubeFFT( )
  830. # canvas.ax1 = canvas.fig.add_axes([.1, .1, .8, .8])
  831. canvas.ax1 = canvas.fig.add_axes([.1, .1, .2, .8])
  832. canvas.ax2 = canvas.fig.add_axes([.325, .1, .2, .8])
  833. canvas.ax3 = canvas.fig.add_axes([.55, .1, .2, .8])
  834. canvas.ax4 = canvas.fig.add_axes([.815, .1, .05, .8]) #cb
  835. canvas.ax1.tick_params(axis='both', which='major', labelsize=8)
  836. canvas.ax2.tick_params(axis='both', which='major', labelsize=8)
  837. canvas.ax3.tick_params(axis='both', which='major', labelsize=8)
  838. canvas.ax4.tick_params(axis='both', which='major', labelsize=8)
  839. canvas.ax1.set_ylabel("pulse index", fontsize=8)
  840. canvas.ax1.set_xlabel(r"$\omega$ bin", fontsize=8)
  841. canvas.ax2.set_xlabel(r"$\omega$ bin", fontsize=8)
  842. canvas.ax3.set_xlabel(r"$\omega$ bin", fontsize=8)
  843. canvas.ax2.yaxis.set_ticklabels([])
  844. canvas.ax3.yaxis.set_ticklabels([])
  845. #canvas.ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  846. # # Look at pulses
  847. # for pulse in self.DATADICT["PULSES"]:
  848. # for istack in self.DATADICT["stacks"]:
  849. # for ipm in range(0,3):
  850. # canvas.ax1.plot( self.DATADICT[pulse]["CURRENT"][ipm][istack] , label="istack "+str(istack) + " ipm=" + str(ipm) + pulse )
  851. # canvas.draw()
  852. # Create Container for stacks
  853. # sandbox determine pulse sequence again
  854. for pulse in self.DATADICT["PULSES"]:
  855. for ichan in self.DATADICT[pulse]["chan"]:
  856. #for ipm in range(10,11):
  857. CONTAINER = {}
  858. CONTAINER["Cycle 1"] = [] # These are actually subtracted cycles... v+ - v
  859. CONTAINER["Cycle 2"] = []
  860. for istack in self.DATADICT["stacks"]:
  861. #canvas.ax1.clear()
  862. ipm = 8
  863. #for ipm in range(self.DATADICT["nPulseMoments"]):
  864. #canvas.ax1.matshow( np.real(self.DATADICT[pulse][ichan]["FFT"][istack]), aspect='auto' )
  865. #canvas.draw()
  866. if not istack%4%4:
  867. # phase cycle 4, aligned with 1 after sub
  868. CONTAINER["Cycle 1"].append(-self.DATADICT[pulse][ichan]["FFT"][istack])
  869. #canvas.ax1.plot( self.DATADICT[pulse]["TIMES"], -self.DATADICT[pulse][ichan][ipm][istack], label="istack "+str(istack)+ " " + pulse )
  870. elif not istack%4%3:
  871. # phase cycle 3, aligned with 2 after sub
  872. CONTAINER["Cycle 2"].append(-self.DATADICT[pulse][ichan]["FFT"][istack])
  873. #canvas.ax1.plot( self.DATADICT[pulse]["TIMES"], -self.DATADICT[pulse][ichan][ipm][istack], label="istack "+str(istack)+ " " + pulse )
  874. elif not istack%4%2:
  875. # phase cycle 2
  876. CONTAINER["Cycle 2"].append( self.DATADICT[pulse][ichan]["FFT"][istack])
  877. #canvas.ax1.plot( self.DATADICT[pulse]["TIMES"], self.DATADICT[pulse][ichan][ipm][istack], label="istack "+str(istack)+ " " + pulse )
  878. else:
  879. # phase cycle 1
  880. CONTAINER["Cycle 1"].append( self.DATADICT[pulse][ichan]["FFT"][istack])
  881. #canvas.ax1.plot( self.DATADICT[pulse]["TIMES"], self.DATADICT[pulse][ichan][ipm][istack], label="istack "+str(istack)+ " " + pulse )
  882. #canvas.ax1.matshow(np.array(np.average(self.DATADICT[pulse][ichan]["FFT"]), axis=2), aspect='auto' )
  883. #canvas.ax1.plot( self.DATADICT[pulse]["PULSE_TIMES"], self.DATADICT[pulse]["CURRENT"][ipm][istack] , color='black', label="istack "+str(istack) )
  884. #canvas.ax1.plot( self.DATADICT[pulse]["CURRENT"][ipm][istack] , label="istack "+str(istack) + " iFID" + str(iFID) )
  885. #canvas.ax1.plot( self.DATADICT[pulse]["TIMES"], self.DATADICT[pulse][ichan][ipm][istack], label="istack "+str(istack)+ " " + pulse )
  886. #canvas.ax1.legend(prop={'size':6})
  887. #canvas.draw()
  888. # Boostrap
  889. # stack.
  890. #scipy.random.shuffle(x)
  891. # Stack and calculate the pooled variance (http://en.wikipedia.org/wiki/Pooled_variance)
  892. """ All this phase cycling wreaks havoc on a normal calculation of std. and variance. Instead, we resort to calculating
  893. a pooled variance. In this assumption is that the precision of the measurment is constant. This is a poor choice for
  894. any type of moving sensor.
  895. """
  896. # if a window filter has been applied
  897. #self.WINDOW
  898. #self.IWindowStart
  899. #self.iWindowEnd
  900. #self.FFTtimes
  901. CONTAINER = .5*(np.array(CONTAINER["Cycle 2"]) - np.array(CONTAINER["Cycle 1"]))
  902. print ("container shape", np.shape( CONTAINER), self.iWindowStart+1, self.iWindowEnd-1)
  903. dmin = np.min(np.abs(np.average(np.array(CONTAINER)[:,:,self.iWindowStart+1:self.iWindowEnd-1], axis=0)))
  904. dmax = np.max(np.abs(np.average(np.array(CONTAINER)[:,:,self.iWindowStart+1:self.iWindowEnd-1], axis=0)))
  905. mn = canvas.ax1.matshow( 20.*np.log10(np.abs(np.average(np.array(CONTAINER)[:,:, self.iWindowStart+1:self.iWindowEnd-1], axis=0))), aspect='auto', vmin=-120, vmax=-40)
  906. #mn = canvas.ax1.matshow(20.*np.log10(XA[:,istart:iend+1]), aspect='auto', vmax=-40, vmin=-120) #, norm=LogNorm())
  907. canvas.ax2.matshow( 20*np.log10(np.std(np.real(np.array(CONTAINER)[:,:,self.iWindowStart+1:self.iWindowEnd-1]), axis=0)), aspect='auto', vmin=-120, vmax=-40)
  908. canvas.ax3.matshow( 20*np.log10(np.std(np.imag(np.array(CONTAINER)[:,:,self.iWindowStart+1:self.iWindowEnd-1]), axis=0)), aspect='auto', vmin=-120, vmax=-40)
  909. #canvas.ax1.legend(prop={'size':6})
  910. cb1 = mpl.colorbar.Colorbar(canvas.ax4, mn)
  911. cb1.ax.tick_params(labelsize=8)
  912. cb1.set_label("power [dB]", fontsize=8)
  913. canvas.ax1.xaxis.set_major_locator(MaxNLocator(4))
  914. canvas.ax2.xaxis.set_major_locator(MaxNLocator(4))
  915. canvas.ax3.xaxis.set_major_locator(MaxNLocator(4))
  916. canvas.draw()
  917. self.doneTrigger.emit()
  918. def effectivePulseMoment(self, cv, canvas):
  919. canvas.reAxH(2)
  920. nstack = len(self.DATADICT["stacks"])
  921. #canvas.ax1.set_yscale('log')
  922. for pulse in self.DATADICT["PULSES"]:
  923. self.DATADICT[pulse]["qeff"] = {}
  924. self.DATADICT[pulse]["q_nu"] = {}
  925. for ipm in range(self.DATADICT["nPulseMoments"]):
  926. self.DATADICT[pulse]["qeff"][ipm] = {}
  927. self.DATADICT[pulse]["q_nu"][ipm] = {}
  928. #canvas.ax1.clear()
  929. scolours = np.array([0.,0.,1.])
  930. for istack in self.DATADICT["stacks"]:
  931. #self.DATADICT[pulse]["PULSE_TIMES"]
  932. x = self.DATADICT[pulse]["CURRENT"][ipm][istack]
  933. X = np.fft.rfft(x)
  934. v = np.fft.fftfreq(len(x), self.dt)
  935. v = v[0:len(X)]
  936. v[-1] = np.abs(v[-1])
  937. # calculate effective current/moment
  938. I0 = np.abs(X)/len(X)
  939. qeff = I0 * (self.DATADICT[pulse]["PULSE_TIMES"][-1]-self.DATADICT[pulse]["PULSE_TIMES"][0])
  940. canvas.ax1.set_title(r"pulse moment index " +str(ipm), fontsize=10)
  941. canvas.ax1.set_xlabel(r"$\nu$ [Hz]", fontsize=8)
  942. canvas.ax1.set_ylabel(r"$q_{eff}$ [A$\cdot$sec]", fontsize=8)
  943. if nstack > 1:
  944. canvas.ax1.plot(v, qeff, color=scolours+istack*np.array((0.,1./(nstack+1.),-1./(nstack+1.)) )) # eff current
  945. else:
  946. canvas.ax1.plot(v, qeff, color=scolours) # eff current
  947. self.DATADICT[pulse]["qeff"][ipm][istack] = qeff
  948. self.DATADICT[pulse]["q_nu"][ipm][istack] = v
  949. canvas.draw()
  950. percent = int(1e2* (float)((istack)+ipm*self.DATADICT["nPulseMoments"]) /
  951. (float)(len(self.DATADICT["PULSES"])*self.DATADICT["nPulseMoments"]*nstack))
  952. self.progressTrigger.emit(percent)
  953. canvas.draw()
  954. self.plotQeffNu(cv, canvas.ax2)
  955. canvas.draw()
  956. self.doneTrigger.emit()
  957. def plotQeffNu(self, cv, ax):
  958. ####################################
  959. # TODO label fid1 and fid2, and make a legend, and colour by pulse
  960. scolours = ['blue','green']
  961. nstack = len(self.DATADICT["stacks"])
  962. iFID = 0
  963. for pulse in self.DATADICT["PULSES"]:
  964. self.DATADICT[pulse]["Q"] = np.zeros( (self.DATADICT["nPulseMoments"], len(self.DATADICT["stacks"])) )
  965. ilabel = True
  966. for ipm in range(self.DATADICT["nPulseMoments"]):
  967. scolours = np.array([0.,0.,1.])
  968. istack = 0
  969. for stack in self.DATADICT["stacks"]:
  970. # find index
  971. icv = int (round(cv / self.DATADICT[pulse]["q_nu"][ipm][stack][1]))
  972. self.DATADICT[pulse]["Q"][ipm,istack] = self.DATADICT[pulse]["qeff"][ipm][stack][icv]
  973. if ilabel:
  974. ax.scatter(ipm, self.DATADICT[pulse]["qeff"][ipm][stack][icv], facecolors='none', edgecolors=scolours, label=(str(pulse)))
  975. ilabel = False
  976. else:
  977. ax.scatter(ipm, self.DATADICT[pulse]["qeff"][ipm][stack][icv], facecolors='none', edgecolors=scolours)
  978. scolours += np.array((0,1./(nstack+1),-1/(nstack+1.)))
  979. percent = int(1e2* (float)((istack)+ipm*self.DATADICT["nPulseMoments"]) /
  980. (float)(len(self.DATADICT["PULSES"])*self.DATADICT["nPulseMoments"]*nstack))
  981. self.progressTrigger.emit(percent)
  982. istack += 1
  983. iFID += 1
  984. ax.set_xlabel(r"pulse moment index", fontsize=8)
  985. ax.set_ylabel(r"$q_{eff}$ [A$\cdot$sec]", fontsize=8)
  986. ax.set_yscale('log')
  987. ax.set_xlim(0, ax.get_xlim()[1])
  988. ax.legend(loc='upper right', scatterpoints = 1, prop={'size':6})
  989. def enableDSP(self):
  990. self.enableDSPTrigger.emit()
  991. def adaptiveFilter(self, M, flambda, truncate, mu, PCA, canvas):
  992. canvas.reAx2(shx=False, shy=False)
  993. canvas.ax2.tick_params(axis='both', which='major', labelsize=8)
  994. canvas.ax2.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  995. canvas.ax2.tick_params(axis='both', which='major', labelsize=8)
  996. canvas.ax2.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  997. if truncate:
  998. itrunc =(int) ( round( 1e-3*truncate*self.samp ) )
  999. print( "adaptive filter size", 1e3*self.dt*M, " [ms]" )
  1000. Filt = adapt.AdaptiveFilter(flambda)
  1001. H = {}
  1002. for pulse in self.DATADICT["PULSES"]:
  1003. H[pulse] = {}
  1004. for ichan in self.DATADICT[pulse]["chan"]:
  1005. print("setting H", pulse, ichan)
  1006. H[pulse][ichan] = np.zeros(M)
  1007. iFID = 0
  1008. # original ordering...
  1009. #for pulse in self.DATADICT["PULSES"]:
  1010. # for ipm in range(self.DATADICT["nPulseMoments"]):
  1011. # for istack in self.DATADICT["stacks"]:
  1012. # This order makes more sense, same as data collection, verify
  1013. for istack in self.DATADICT["stacks"]:
  1014. for ipm in range(self.DATADICT["nPulseMoments"]):
  1015. for pulse in self.DATADICT["PULSES"]:
  1016. canvas.ax1.clear()
  1017. canvas.ax2.clear()
  1018. for ichan in self.DATADICT[pulse]["chan"]:
  1019. #H = np.zeros(M)
  1020. canvas.ax1.plot( self.DATADICT[pulse]["TIMES"], 1e9* self.DATADICT[pulse][ichan][ipm][istack],\
  1021. label = "noisy")
  1022. RX = []
  1023. for irchan in self.DATADICT[pulse]["rchan"]:
  1024. RX.append(self.DATADICT[pulse][irchan][ipm][istack][::-1])
  1025. if all(H[pulse][ichan]) == 0:
  1026. # call twice to reconcile filter wind-up
  1027. [e,H[pulse][ichan]] = Filt.adapt_filt_Ref( self.DATADICT[pulse][ichan][ipm][istack][::-1],\
  1028. RX,\
  1029. M, mu, PCA, flambda, H[pulse][ichan])
  1030. [e,H[pulse][ichan]] = Filt.adapt_filt_Ref( self.DATADICT[pulse][ichan][ipm][istack][::-1],\
  1031. RX,\
  1032. M, mu, PCA, flambda, H[pulse][ichan])
  1033. else:
  1034. [e,H[pulse][ichan]] = Filt.adapt_filt_Ref( self.DATADICT[pulse][ichan][ipm][istack][::-1],\
  1035. RX,\
  1036. M, mu, PCA, flambda, H[pulse][ichan])
  1037. # replace
  1038. if truncate:
  1039. canvas.ax1.plot( self.DATADICT[pulse]["TIMES"][0:itrunc], 1e9* e[::-1][0:itrunc],\
  1040. label = pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " ichan=" + str(ichan))
  1041. self.DATADICT[pulse][ichan][ipm][istack] = e[::-1][0:itrunc]
  1042. else:
  1043. canvas.ax1.plot( self.DATADICT[pulse]["TIMES"], 1e9* e[::-1],\
  1044. label = pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " ichan=" + str(ichan))
  1045. self.DATADICT[pulse][ichan][ipm][istack] = e[::-1]
  1046. canvas.ax2.plot( H[pulse][ichan] , label="taps")
  1047. canvas.ax1.legend(prop={'size':6})
  1048. canvas.ax2.legend(prop={'size':6})
  1049. canvas.ax1.set_xlabel(r"time [s]", fontsize=8)
  1050. canvas.ax1.set_ylabel(r"signal [nV]", fontsize=8)
  1051. canvas.ax2.set_xlabel(r"filter index", fontsize=8)
  1052. canvas.ax2.set_ylabel(r"scale factor", fontsize=8)
  1053. canvas.draw()
  1054. # truncate the reference channels too, in case you still need them for something.
  1055. # Otherwise they are no longer aligned with the data
  1056. for rchan in self.DATADICT[pulse]["rchan"]:
  1057. if truncate:
  1058. self.DATADICT[pulse][rchan][ipm][istack] = self.DATADICT[pulse][rchan][ipm][istack][0:itrunc]
  1059. #percent = (int)(1e2*((float)(iFID*self.DATADICT["nPulseMoments"]+(ipm))/( len(self.DATADICT["PULSES"])*self.nPulseMoments)))
  1060. percent = (int)(1e2*((float)(istack*self.DATADICT["nPulseMoments"]+(ipm))/( len(self.DATADICT["PULSES"])*self.nPulseMoments*(len(self.DATADICT["stacks"])+1) )))
  1061. self.progressTrigger.emit(percent)
  1062. # # why is this loop here, istack is not part of rest?
  1063. # for istack in self.DATADICT["stacks"]:
  1064. # if truncate:
  1065. # self.DATADICT[pulse]["TIMES"] = self.DATADICT[pulse]["TIMES"][0:itrunc]
  1066. # percent = (int)(1e2*((float)(iFID*self.DATADICT["nPulseMoments"]+(ipm))/( len(self.DATADICT["PULSES"])*self.nPulseMoments)))
  1067. # self.progressTrigger.emit(percent)
  1068. # iFID += 1
  1069. if truncate:
  1070. self.DATADICT[pulse]["TIMES"] = self.DATADICT[pulse]["TIMES"][0:itrunc]
  1071. self.doneTrigger.emit()
  1072. self.updateProcTrigger.emit()
  1073. #self.plotFT(canvas)
  1074. def plotFT(self, canvas, istart=0, iend=0):
  1075. try:
  1076. canvas.fig.clear()
  1077. except:
  1078. pass
  1079. canvas.ax1 = canvas.fig.add_axes([.1, .1, .65, .8])
  1080. canvas.ax1c = canvas.fig.add_axes([.8, .1, .05, .8])
  1081. canvas.ax1.tick_params(axis='both', which='major', labelsize=8)
  1082. for pulse in self.DATADICT["PULSES"]:
  1083. for istack in self.DATADICT["stacks"]:
  1084. for ichan in self.DATADICT[pulse]["chan"]:
  1085. # FFT of stack
  1086. XA = np.zeros((self.DATADICT["nPulseMoments"] , len(self.DATADICT[pulse][ichan][0][istack])/2+1))
  1087. nu = np.fft.fftfreq(self.DATADICT[pulse][ichan][0][istack].size, d=self.dt)
  1088. nu[-1] *= -1
  1089. df = nu[1]
  1090. of = 0
  1091. if istart:
  1092. of = nu[istart]
  1093. def freqlabel(x, pos):
  1094. return '%1.0f' %(of + x*df)
  1095. formatter = FuncFormatter(freqlabel)
  1096. canvas.ax1.clear()
  1097. for ipm in range(self.DATADICT["nPulseMoments"]):
  1098. X = np.fft.rfft(self.DATADICT[pulse][ichan][ipm][istack])
  1099. XA[ipm,:] = np.abs(X)
  1100. if istart:
  1101. mn = canvas.ax1.matshow(20.*np.log10(XA[:,istart:iend+1]), aspect='auto', vmax=-40, vmin=-120) #, norm=LogNorm())
  1102. else:
  1103. mn = canvas.ax1.matshow(20.*np.log10(XA), aspect='auto', vmax=-40, vmin=-120) #, norm=LogNorm())
  1104. smin = np.min(20.*np.log10(XA))
  1105. smax = np.max(20.*np.log10(XA))
  1106. canvas.ax1.xaxis.set_major_formatter(formatter)
  1107. cb1 = mpl.colorbar.Colorbar(canvas.ax1c, mn)
  1108. cb1.ax.tick_params(labelsize=8)
  1109. cb1.set_label("signal [dB]", fontsize=8)
  1110. canvas.ax1.set_xlabel(r"$\nu$ [Hz]", fontsize=10)
  1111. canvas.ax1.set_ylabel(r"$q_{index}$", fontsize=10)
  1112. canvas.draw()
  1113. def plotFT(self, canvas, istart=0, iend=0):
  1114. try:
  1115. canvas.fig.clear()
  1116. except:
  1117. pass
  1118. canvas.ax1 = canvas.fig.add_axes([.1, .1, .65, .8])
  1119. canvas.ax1c = canvas.fig.add_axes([.8, .1, .05, .8])
  1120. canvas.ax1.tick_params(axis='both', which='major', labelsize=8)
  1121. for pulse in self.DATADICT["PULSES"]:
  1122. for istack in self.DATADICT["stacks"]:
  1123. for ichan in self.DATADICT[pulse]["chan"]:
  1124. # FFT of stack
  1125. XA = np.zeros((self.DATADICT["nPulseMoments"] , len(self.DATADICT[pulse][ichan][0][istack])//2+1))
  1126. nu = np.fft.fftfreq(self.DATADICT[pulse][ichan][0][istack].size, d=self.dt)
  1127. nu[-1] *= -1
  1128. df = nu[1]
  1129. of = 0
  1130. if istart:
  1131. of = nu[istart]
  1132. def freqlabel(x, pos):
  1133. return '%1.0f' %(of + x*df)
  1134. formatter = FuncFormatter(freqlabel)
  1135. canvas.ax1.clear()
  1136. for ipm in range(self.DATADICT["nPulseMoments"]):
  1137. X = np.fft.rfft(self.DATADICT[pulse][ichan][ipm][istack])
  1138. XA[ipm,:] = np.abs(X)
  1139. if istart:
  1140. mn = canvas.ax1.matshow(20.*np.log10(XA[:,istart:iend+1]), aspect='auto', vmax=-40, vmin=-120, cmap='viridis') #, norm=LogNorm())
  1141. else:
  1142. mn = canvas.ax1.matshow(20.*np.log10(XA), aspect='auto', vmax=-40, vmin=-120, cmap='viridis') #, norm=LogNorm())
  1143. canvas.ax1.xaxis.set_major_formatter(formatter)
  1144. cb1 = mpl.colorbar.Colorbar(canvas.ax1c, mn)
  1145. cb1.ax.tick_params(labelsize=8)
  1146. cb1.set_label("signal [dB]", fontsize=8)
  1147. canvas.ax1.set_xlabel(r"$\nu$ [Hz]", fontsize=10)
  1148. canvas.ax1.set_ylabel(r"$q_{index}$", fontsize=10)
  1149. canvas.draw()
  1150. def dataCubeFFT(self):
  1151. """
  1152. Performs FFT on entire cube of DATA, and REFERENCE channels, but not pulse currents,
  1153. Results are saved to a new field in the data structure
  1154. The GMR varies phase as a function of pulse moment index, so that the first pusle moment is zero phase,
  1155. the second is pi/2 the third is zero. This method corrects for this, so that all pulse moments are in phase.
  1156. Technically we may not want to do this, if there is some system response that this cycles away, and we lose track of
  1157. how many of each cycle we have, could this be problomatic? I think it will come out in the wash as we keep track of the
  1158. rest of the phase cycles. Holy phase cycling batman.
  1159. """
  1160. for pulse in self.DATADICT["PULSES"]:
  1161. for ichan in np.append(self.DATADICT[pulse]["chan"], self.DATADICT[pulse]["rchan"]):
  1162. # FFT of stack
  1163. self.DATADICT[pulse][ichan]["FFT"] = {}
  1164. self.DATADICT[pulse][ichan]["FFT"]["nu"] = np.fft.fftfreq(self.DATADICT[pulse][ichan][0][self.DATADICT["stacks"][0]].size, d=self.dt)
  1165. self.DATADICT[pulse][ichan]["FFT"]["nu"][-1] *= -1
  1166. for istack in self.DATADICT["stacks"]:
  1167. self.DATADICT[pulse][ichan]["FFT"][istack] = np.zeros((self.DATADICT["nPulseMoments"] , len(self.DATADICT[pulse][ichan][0][istack])//2+1), dtype=complex)
  1168. for ipm in range(self.DATADICT["nPulseMoments"]):
  1169. # Mod works for FID pulse sequences, TODO generalize this for 4 phase T1, etc..
  1170. #mod = (-1)**(ipm%2) * (-1)**(istack%2)
  1171. self.DATADICT[pulse][ichan]["FFT"][istack][ipm,:] = np.fft.rfft( self.DATADICT[pulse][ichan][ipm][istack] )
  1172. #if ipm%2:
  1173. # odd, phase cycled from previous
  1174. # self.DATADICT[pulse][ichan]["FFT"][istack][ipm,:] = np.fft.rfft(-self.DATADICT[pulse][ichan][ipm][istack])
  1175. #else:
  1176. # even, we define as zero phase, first pulse moment has this
  1177. # self.DATADICT[pulse][ichan]["FFT"][istack][ipm,:] = np.fft.rfft(self.DATADICT[pulse][ichan][ipm][istack])
  1178. def adaptiveFilterFD(self, ftype, band, centre, canvas):
  1179. try:
  1180. canvas.fig.clear()
  1181. except:
  1182. pass
  1183. canvas.ax1 = canvas.fig.add_axes([.1, .5, .7, .4])
  1184. canvas.ax1c = canvas.fig.add_axes([.85, .5, .05, .4])
  1185. canvas.ax1.tick_params(axis='both', which='major', labelsize=8)
  1186. #canvas.ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  1187. canvas.ax2 = canvas.fig.add_axes([.1, .05, .7, .4])
  1188. canvas.ax2c = canvas.fig.add_axes([.85, .05, .05, .4])
  1189. canvas.ax2.tick_params(axis='both', which='major', labelsize=8)
  1190. #canvas.ax2.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  1191. self.dataCubeFFT()
  1192. Filt = adapt.AdaptiveFilter(0.)
  1193. for pulse in self.DATADICT["PULSES"]:
  1194. # Compute window function and dimensions
  1195. [WINDOW, nd, wstart, wend, dead] = self.computeWindow(pulse, band, centre, ftype)
  1196. for istack in self.DATADICT["stacks"]:
  1197. for ichan in self.DATADICT[pulse]["chan"]:
  1198. # FFT of stack
  1199. nd = len(self.DATADICT[pulse][ichan][0][istack])
  1200. XX = np.zeros((self.DATADICT["nPulseMoments"] , len(self.DATADICT[pulse][ichan][0][istack])//2+1), dtype=complex)
  1201. nu = np.fft.fftfreq(self.DATADICT[pulse][ichan][0][istack].size, d=self.dt)
  1202. nu[-1] *= -1
  1203. #nu = self.DATADICT[pulse][ichan]["FFT"]["nu"]
  1204. def freqlabel(x, pos):
  1205. return '%1.0f' %((wstart)*nu[1] + x*nu[1])
  1206. formatter = FuncFormatter(freqlabel)
  1207. canvas.ax1.clear()
  1208. for ipm in range(self.DATADICT["nPulseMoments"]):
  1209. X = np.fft.rfft(self.DATADICT[pulse][ichan][ipm][istack])
  1210. XX[ipm,:] = X
  1211. XX = XX*WINDOW
  1212. XX = XX[:,wstart:wend]
  1213. smin = np.min(20.*np.log10(np.abs(XX)))
  1214. smax = np.max(20.*np.log10(np.abs(XX)))
  1215. #if smin != smin:
  1216. smax = -40
  1217. smin = -120
  1218. mn = canvas.ax1.matshow(20.*np.log10(np.abs(XX)), aspect='auto', vmin=smin, vmax=smax) #, norm=LogNorm())
  1219. canvas.ax1.xaxis.set_major_formatter(formatter)
  1220. cb1 = mpl.colorbar.Colorbar(canvas.ax1c, mn)
  1221. RX = []
  1222. for ichan in self.DATADICT[pulse]["rchan"]:
  1223. R = np.zeros((self.DATADICT["nPulseMoments"] , len(self.DATADICT[pulse][ichan][0][istack])//2+1), dtype=complex)
  1224. for ipm in range(self.DATADICT["nPulseMoments"]):
  1225. R[ipm,:] = np.fft.rfft(self.DATADICT[pulse][ichan][ipm][istack])
  1226. RX.append(R[:,wstart:wend])
  1227. XC = Filt.transferFunctionFFT(XX, RX)
  1228. # TODO inverse FFT, but we need to map back to origional matrix size
  1229. #for ichan in self.DATADICT[pulse]["chan"]:
  1230. # for ipm in range(self.DATADICT["nPulseMoments"]):
  1231. # self.DATADICT[pulse][ichan][ipm][istack] = np.fft.irfft(XC[] , nd)
  1232. mc = canvas.ax2.matshow(20.*np.log10(np.abs(XC)), aspect='auto', vmin=smin, vmax=smax) #, norm=LogNorm())
  1233. cb2 = mpl.colorbar.Colorbar(canvas.ax2c, mc)
  1234. cmin = np.min(20.*np.log10(np.abs(XC)))
  1235. cmax = np.max(20.*np.log10(np.abs(XC)))
  1236. canvas.ax2.xaxis.set_major_formatter(formatter)
  1237. #canvas.ax2.colorbar(mn)
  1238. canvas.draw()
  1239. ##############################3
  1240. # TODO inverse FFT to get the damn data back!!!
  1241. # self.progressTrigger.emit(percent)
  1242. # #label = "iFID="+str(iFID) + " ipm=" + str(ipm) + " istack=" + str(istack) + " ichan=" + str(ichan))
  1243. self.doneTrigger.emit()
  1244. def findSpikes(self, x, width, threshold, rollOn):
  1245. import scipy.ndimage as im
  1246. spikes = np.zeros( len(x) )
  1247. med = im.median_filter(x, width,mode='nearest')
  1248. std = np.std(x)
  1249. spikes = (np.abs(x-med) > threshold * std)
  1250. return np.array(np.where(spikes[rollOn::])) + rollOn
  1251. # def despike(self, width, threshold, itype, rollOn, win, canvas):
  1252. # from scipy import interpolate
  1253. # """ This was a stab at a despike filter. Better results were achieved using the SmartStack approach
  1254. # """
  1255. # try:
  1256. # canvas.fig.clear()
  1257. # except:
  1258. # pass
  1259. #
  1260. # canvas.ax1 = canvas.fig.add_axes([.125,.1,.725,.8])
  1261. # canvas.ax1.tick_params(axis='both', which='major', labelsize=8)
  1262. # canvas.ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  1263. # iFID = 0
  1264. # for pulse in self.DATADICT["PULSES"]:
  1265. # for ipm in range(self.DATADICT["nPulseMoments"]):
  1266. # for istack in self.DATADICT["stacks"]:
  1267. # canvas.ax1.clear()
  1268. # for ichan in np.append(self.DATADICT[pulse]["chan"], self.DATADICT[pulse]["rchan"]):
  1269. # x = self.findSpikes(self.DATADICT[pulse][ichan][ipm][istack], width, threshold, rollOn)
  1270. # canvas.ax1.plot( self.DATADICT[pulse]["TIMES"], self.DATADICT[pulse][ichan][ipm][istack],
  1271. # label = pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " ichan=" + str(ichan))
  1272. # canvas.ax1.plot( self.DATADICT[pulse]["TIMES"][x], self.DATADICT[pulse][ichan][ipm][istack][x], '.', color='red' , markersize=6 )
  1273. #
  1274. # FIXED = np.zeros(len(x[0]))
  1275. # ii = 0
  1276. # for spike in np.array(x[0]).tolist():
  1277. # f = interpolate.interp1d(np.delete(self.DATADICT[pulse]["TIMES"][spike-win/2:spike+win/2], x[0]-(spike-win/2)), \
  1278. # np.delete(self.DATADICT[pulse][ichan][ipm][istack][spike-win/2:spike+win/2], x[0]-(spike-win/2)), itype)
  1279. # FIXED[ii] = f(self.DATADICT[pulse]["TIMES"][spike])
  1280. # ii += 1
  1281. # canvas.ax1.plot( self.DATADICT[pulse]["TIMES"][x[0]] , FIXED, '.', color='black' , markersize=4 )
  1282. # self.DATADICT[pulse][ichan][ipm][istack][x[0]] = FIXED
  1283. #
  1284. # canvas.ax1.legend(prop={'size':6})
  1285. # canvas.draw()
  1286. # percent = (int)(1e2*((float)(iFID*self.DATADICT["nPulseMoments"]+(ipm))/( len(self.DATADICT["PULSES"])*self.nPulseMoments)))
  1287. # self.progressTrigger.emit(percent)
  1288. # iFID += 1
  1289. # self.doneTrigger.emit()
  1290. def designFilter(self, cf, PB, SB, gpass, gstop, ftype, canvas):
  1291. ''' cf is central frequency
  1292. pb is pass band
  1293. sb is stop band
  1294. '''
  1295. TS = (cf) / (.5/self.dt)
  1296. PB = PB / (.5/self.dt) # 1/2 width pass band Muddy Creek
  1297. SB = SB / (.5/self.dt) # 1/2 width stop band Muddy Creek
  1298. # if butterworth
  1299. #[bord, wn] = signal.buttord([TS-PB,TS+PB], [TS-SB,TS+SB], 1e-1, 5.)
  1300. if ftype=="Butterworth":
  1301. [bord, wn] = signal.buttord([TS-PB,TS+PB], [TS-SB,TS+SB], gpass, gstop)
  1302. [self.filt_b, self.filt_a] = signal.butter(bord, wn, btype='bandpass', output='ba')
  1303. [self.filt_z, self.filt_p, self.filt_k] = signal.butter(bord, wn, btype='band', output='zpk')
  1304. elif ftype == "Chebychev Type II":
  1305. [bord, wn] = signal.cheb2ord([TS-PB,TS+PB], [TS-SB,TS+SB], gpass, gstop)
  1306. [self.filt_b, self.filt_a] = signal.cheby2(bord, gstop, wn, btype='bandpass', output='ba')
  1307. [self.filt_z, self.filt_p, self.filt_k] = signal.cheby2(bord, gstop, wn, btype='band', output='zpk')
  1308. elif ftype == "Elliptic":
  1309. [bord, wn] = signal.ellipord([TS-PB,TS+PB], [TS-SB,TS+SB], gpass, gstop)
  1310. [self.filt_b, self.filt_a] = signal.ellip(bord, gpass, gstop, wn, btype='bandpass', output='ba')
  1311. [self.filt_z, self.filt_p, self.filt_k] = signal.ellip(bord, gpass, gstop, wn, btype='band', output='zpk')
  1312. # if cheby2
  1313. impulse = self.mfreqz2(self.filt_b, self.filt_a, canvas)
  1314. self.fe = -5
  1315. for it in range(len(impulse[0])):
  1316. if abs(impulse[1][0][it][0]) >= .1 * gpass:# gpass:
  1317. self.fe = impulse[0][it]
  1318. return [bord, self.fe]
  1319. def downsample(self, truncate, dec, canvas):
  1320. """ Downsamples and truncates the raw signal.
  1321. """
  1322. canvas.reAx2()
  1323. canvas.ax1.tick_params(axis='both', which='major', labelsize=8)
  1324. canvas.ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  1325. canvas.ax2.tick_params(axis='both', which='major', labelsize=8)
  1326. canvas.ax2.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  1327. self.samp /= dec
  1328. self.dt = 1./self.samp
  1329. if truncate:
  1330. itrunc = (int)( 1e-3*truncate*self.samp )
  1331. iFID = 0
  1332. for pulse in self.DATADICT["PULSES"]:
  1333. for ipm in range(self.DATADICT["nPulseMoments"]):
  1334. for istack in self.DATADICT["stacks"]:
  1335. canvas.ax1.clear()
  1336. canvas.ax2.clear()
  1337. for ichan in np.append(self.DATADICT[pulse]["chan"], self.DATADICT[pulse]["rchan"]):
  1338. # trim off indices that don't divide evenly
  1339. ndi = np.shape(self.DATADICT[pulse][ichan][ipm][istack])[0]%dec
  1340. if ndi:
  1341. [self.DATADICT[pulse][ichan][ipm][istack], RSTIMES] = signal.resample(self.DATADICT[pulse][ichan][ipm][istack][0:-ndi],\
  1342. len(self.DATADICT[pulse][ichan][ipm][istack][0:-ndi])//dec,\
  1343. self.DATADICT[pulse]["TIMES"][0:-ndi], window='hamm')
  1344. else:
  1345. [self.DATADICT[pulse][ichan][ipm][istack], RSTIMES] = signal.resample(self.DATADICT[pulse][ichan][ipm][istack],\
  1346. len(self.DATADICT[pulse][ichan][ipm][istack])//dec,\
  1347. self.DATADICT[pulse]["TIMES"], window='hamm')
  1348. if truncate:
  1349. self.DATADICT[pulse][ichan][ipm][istack] = self.DATADICT[pulse][ichan][ipm][istack][0:itrunc]
  1350. RSTIMES = RSTIMES[0:itrunc]
  1351. for ichan in self.DATADICT[pulse]["chan"]:
  1352. canvas.ax2.plot( RSTIMES, 1e9*self.DATADICT[pulse][ichan][ipm][istack], \
  1353. label = pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " ichan=" + str(ichan))
  1354. for ichan in self.DATADICT[pulse]["rchan"]:
  1355. canvas.ax1.plot( RSTIMES, 1e9*self.DATADICT[pulse][ichan][ipm][istack], \
  1356. label = pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " ichan=" + str(ichan))
  1357. canvas.ax1.set_xlabel(r"time [s]", fontsize=8)
  1358. canvas.ax1.set_ylabel(r"signal [nV]", fontsize=8)
  1359. canvas.ax2.set_xlabel(r"time [s]", fontsize=8)
  1360. canvas.ax2.set_ylabel(r"signal [nV]", fontsize=8)
  1361. canvas.ax1.legend(prop={'size':6})
  1362. canvas.ax2.legend(prop={'size':6})
  1363. canvas.ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  1364. canvas.ax2.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  1365. canvas.draw()
  1366. percent = (int)(1e2*((float)(iFID*self.DATADICT["nPulseMoments"]+(ipm))/( len(self.DATADICT["PULSES"])*self.nPulseMoments)))
  1367. self.progressTrigger.emit(percent)
  1368. iFID += 1
  1369. self.DATADICT[pulse]["TIMES"] = RSTIMES
  1370. #####################################
  1371. # resample pulse data
  1372. for pulse in self.DATADICT["PULSES"]:
  1373. for ipm in range(self.DATADICT["nPulseMoments"]):
  1374. for istack in self.DATADICT["stacks"]:
  1375. ndi = np.shape(self.DATADICT[pulse]["CURRENT"][ipm][istack])[0]%dec
  1376. if ndi:
  1377. [self.DATADICT[pulse]["CURRENT"][ipm][istack], RSPTIMES] = signal.resample(self.DATADICT[pulse]["CURRENT"][ipm][istack][0:-ndi],\
  1378. len(self.DATADICT[pulse]["CURRENT"][ipm][istack][0:-ndi])//dec,\
  1379. self.DATADICT[pulse]["PULSE_TIMES"][0:-ndi], window='hamm')
  1380. else:
  1381. [self.DATADICT[pulse]["CURRENT"][ipm][istack], RSPTIMES] = signal.resample(self.DATADICT[pulse]["CURRENT"][ipm][istack],\
  1382. len(self.DATADICT[pulse]["CURRENT"][ipm][istack])//dec,\
  1383. self.DATADICT[pulse]["PULSE_TIMES"], window='hamm')
  1384. self.DATADICT[pulse]["PULSE_TIMES"] = RSPTIMES
  1385. self.doneTrigger.emit()
  1386. self.updateProcTrigger.emit()
  1387. def computeWindow(self, pulse, band, centre, ftype, canvas=None):
  1388. # Compute window
  1389. nd = len(self.DATADICT[pulse][self.DATADICT[pulse]["chan"][0]][0][self.DATADICT["stacks"][0]])
  1390. fft1 = np.fft.rfft(self.DATADICT[pulse][self.DATADICT[pulse]["chan"][0]][0][self.DATADICT["stacks"][0]])
  1391. freqs = np.fft.fftfreq(nd, self.dt)
  1392. df = freqs[1] - freqs[0]
  1393. N = int((round)(band/df))
  1394. if ftype == "Hamming":
  1395. window = np.hamming(N)
  1396. elif ftype == "Hanning":
  1397. window = np.hanning(N)
  1398. elif ftype == "Rectangular":
  1399. window = np.ones(N)
  1400. elif ftype == "Flat top":
  1401. window = signal.flattop(N)
  1402. else:
  1403. print ("in windowFilter, window type undefined")
  1404. WINDOW = np.zeros(len(fft1))
  1405. ifreq = int(round(centre/df))
  1406. istart = ifreq-len(window)//2
  1407. iend = 0
  1408. if N%2:
  1409. WINDOW[ifreq-N//2:ifreq+N//2+1] = window
  1410. iend = ifreq+N//2+1
  1411. else:
  1412. WINDOW[ifreq-N//2:ifreq+N//2] = window
  1413. iend = ifreq+N//2
  1414. self.WINDOW = WINDOW
  1415. self.iWindowStart = istart
  1416. self.iWindowEnd = iend
  1417. self.FFTtimes = nd
  1418. fft1 = np.fft.irfft(WINDOW)
  1419. # calculate dead time
  1420. self.windead = 0.
  1421. for ift in np.arange(100,0,-1):
  1422. #print( ift, fft1[ift] )
  1423. if fft1[ift] > .001:
  1424. #print ("DEAD TIME", 1e3*self.DATADICT[pulse]["TIMES"][ift] - 1e3*self.DATADICT[pulse]["TIMES"][0] )
  1425. dead = 1e3*self.DATADICT[pulse]["TIMES"][ift] - 1e3*self.DATADICT[pulse]["TIMES"][0]
  1426. self.windead = self.DATADICT[pulse]["TIMES"][ift] - self.DATADICT[pulse]["TIMES"][0]
  1427. break
  1428. if canvas != None:
  1429. canvas.fig.clear()
  1430. canvas.ax1 = canvas.fig.add_axes([.1, .6, .75, .35])
  1431. canvas.ax2 = canvas.fig.add_axes([.1, .1, .75, .35])
  1432. canvas.ax1.plot(WINDOW)
  1433. canvas.ax2.plot( 1e3* self.DATADICT[pulse]["TIMES"][0:100] - 1e3*self.DATADICT[pulse]["TIMES"][0], fft1[0:100] )
  1434. canvas.ax2.set_xlabel("time (ms)")
  1435. canvas.ax2.set_title("IFFT")
  1436. canvas.draw()
  1437. return [WINDOW, nd, istart, iend, dead]
  1438. def windowFilter(self, ftype, band, centre, canvas):
  1439. ###############################
  1440. # Window Filter (Ormsby filter http://www.xsgeo.com/course/filt.htm)
  1441. # apply window
  1442. iFID = 0
  1443. for pulse in self.DATADICT["PULSES"]:
  1444. [WINDOW, nd, istart, iend,dead] = self.computeWindow(pulse, band, centre, ftype)
  1445. for istack in self.DATADICT["stacks"]:
  1446. for ipm in range(self.DATADICT["nPulseMoments"]):
  1447. for ichan in np.append(self.DATADICT[pulse]["chan"], self.DATADICT[pulse]["rchan"]):
  1448. fft = np.fft.rfft( self.DATADICT[pulse][ichan][ipm][istack] )
  1449. fft *= WINDOW
  1450. self.DATADICT[pulse][ichan][ipm][istack] = np.fft.irfft(fft , nd)
  1451. percent = (int)(1e2*((float)(iFID*self.DATADICT["nPulseMoments"]+(ipm))/(len(self.DATADICT["PULSES"])*self.nPulseMoments)))
  1452. self.progressTrigger.emit(percent)
  1453. iFID += 1
  1454. self.plotFT(canvas, istart, iend)
  1455. self.doneTrigger.emit()
  1456. def bandpassFilter(self, canvas, blank, plot=True):
  1457. if plot:
  1458. canvas.reAx2()
  1459. canvas.ax1.tick_params(axis='both', which='major', labelsize=8)
  1460. canvas.ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  1461. canvas.ax2.tick_params(axis='both', which='major', labelsize=8)
  1462. canvas.ax2.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  1463. ife = (int)( max(self.fe, self.windead) * self.samp )
  1464. # Data
  1465. iFID = 0
  1466. for pulse in self.DATADICT["PULSES"]:
  1467. self.DATADICT[pulse]["TIMES"] = self.DATADICT[pulse]["TIMES"][ife:-ife]
  1468. for ipm in range(self.DATADICT["nPulseMoments"]):
  1469. for istack in self.DATADICT["stacks"]:
  1470. canvas.ax1.clear()
  1471. canvas.ax2.clear()
  1472. #for ichan in np.append(self.DATADICT[pulse]["chan"], self.DATADICT[pulse]["rchan"]):
  1473. for ichan in self.DATADICT[pulse]["rchan"]:
  1474. # reflect signal back on itself to reduce gibbs effects on early times
  1475. #nr = len( self.DATADICT[pulse][ichan][ipm][istack] ) - 1 + ife
  1476. #refl = np.append( -1*self.DATADICT[pulse][ichan][ipm][istack][::-1][0:-1], self.DATADICT[pulse][ichan][ipm][istack] )
  1477. #reflfilt = signal.filtfilt( self.filt_b, self.filt_a, refl )
  1478. #self.DATADICT[pulse][ichan][ipm][istack] = reflfilt[nr:-ife]
  1479. # don't reflect
  1480. self.DATADICT[pulse][ichan][ipm][istack] = \
  1481. signal.filtfilt(self.filt_b, self.filt_a, self.DATADICT[pulse][ichan][ipm][istack])[ife:-ife]
  1482. # plot
  1483. if plot:
  1484. canvas.ax1.plot( self.DATADICT[pulse]["TIMES"], 1e9*self.DATADICT[pulse][ichan][ipm][istack], \
  1485. label = pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " rchan=" + str(ichan))
  1486. for ichan in self.DATADICT[pulse]["chan"]:
  1487. # reflect signal back on itself to reduce gibbs effects on early times
  1488. #nr = len( self.DATADICT[pulse][ichan][ipm][istack] ) - 1 + ife
  1489. #refl = np.append( -1*self.DATADICT[pulse][ichan][ipm][istack][::-1][0:-1], self.DATADICT[pulse][ichan][ipm][istack] )
  1490. #reflfilt = signal.filtfilt( self.filt_b, self.filt_a, refl )
  1491. #self.DATADICT[pulse][ichan][ipm][istack] = reflfilt[nr:-ife]
  1492. # don't reflect
  1493. self.DATADICT[pulse][ichan][ipm][istack] = \
  1494. scipy.signal.filtfilt(self.filt_b, self.filt_a, self.DATADICT[pulse][ichan][ipm][istack])[ife:-ife]
  1495. # plot
  1496. if plot:
  1497. canvas.ax2.plot( self.DATADICT[pulse]["TIMES"], 1e9*self.DATADICT[pulse][ichan][ipm][istack], \
  1498. label = "data " + pulse + " ipm=" + str(ipm) + " istack=" + str(istack) + " chan=" + str(ichan))
  1499. if plot:
  1500. canvas.ax1.set_xlabel(r"time [s]", fontsize=8)
  1501. canvas.ax1.set_ylabel(r"signal [nV]", fontsize=8)
  1502. canvas.ax2.set_xlabel(r"time [s]", fontsize=8)
  1503. canvas.ax2.set_ylabel(r"signal [nV]", fontsize=8)
  1504. canvas.ax1.legend(prop={'size':6})
  1505. canvas.ax2.legend(prop={'size':6})
  1506. canvas.draw()
  1507. percent = (int)(1e2*((float)(iFID*self.DATADICT["nPulseMoments"]+(ipm))/(len(self.DATADICT["PULSES"])*self.nPulseMoments)))
  1508. self.progressTrigger.emit(percent)
  1509. iFID += 1
  1510. self.doneTrigger.emit()
  1511. self.updateProcTrigger.emit()
  1512. def loadFIDData(self, base, procStacks, chanin, rchanin, FIDProc, canvas, deadTime, plot):
  1513. '''
  1514. Loads a GMR FID dataset, reads OLD ASCII files
  1515. '''
  1516. import struct
  1517. canvas.reAx3(True,False)
  1518. chan = []
  1519. for ch in chanin:
  1520. chan.append(str(ch))
  1521. rchan = []
  1522. for ch in rchanin:
  1523. rchan.append(str(ch))
  1524. # not in any headers but this has changed, NOT the place to do this. MOVE
  1525. #self.prePulseDelay = 0.01 # delay before pulse
  1526. self.deadTime = deadTime # instrument dead time before measurement
  1527. self.samp = 50000. # in case this is a reproc, these might have
  1528. self.dt = 1./self.samp # changed
  1529. invGain = 1./self.RxGain
  1530. invCGain = self.CurrentGain
  1531. #################################################################################
  1532. # figure out key data indices
  1533. # Pulse
  1534. nps = (int)((self.prePulseDelay)*self.samp)
  1535. npul = (int)(self.pulseLength[0]*self.samp) #+ 100
  1536. # Data
  1537. nds = nps+npul+(int)((self.deadTime)*self.samp); # indice pulse 1 data starts
  1538. nd1 = (int)(1.*self.samp) # samples in first pulse
  1539. #################################################################################
  1540. # Data structures
  1541. PULSES = [FIDProc]
  1542. PULSES = ["Pulse 1"]
  1543. self.DATADICT = {}
  1544. self.DATADICT["nPulseMoments"] = self.nPulseMoments
  1545. self.DATADICT["stacks"] = procStacks
  1546. self.DATADICT["PULSES"] = PULSES
  1547. for pulse in PULSES:
  1548. self.DATADICT[pulse] = {}
  1549. self.DATADICT[pulse]["chan"] = chan # TODO these should not be a subet of pulse! for GMR all
  1550. self.DATADICT[pulse]["rchan"] = rchan # data are consistent
  1551. self.DATADICT[pulse]["CURRENT"] = {}
  1552. for ichan in np.append(chan,rchan):
  1553. self.DATADICT[pulse][ichan] = {}
  1554. for ipm in range(self.nPulseMoments):
  1555. self.DATADICT[pulse][ichan][ipm] = {}
  1556. self.DATADICT[pulse]["CURRENT"][ipm] = {}
  1557. for istack in procStacks:
  1558. self.DATADICT[pulse][ichan][ipm][istack] = np.zeros(3)
  1559. self.DATADICT[pulse]["CURRENT"][ipm][istack] = np.zeros(3)
  1560. ##############################################
  1561. # Read in binary (.lvm) data
  1562. iistack = 0
  1563. hack = False
  1564. for istack in procStacks:
  1565. rawFile = open(base + "_" + str(istack) + ".lvm", 'rb')
  1566. if hack:
  1567. subFile = open(base + "sub" + "_" + str(istack) + ".lvm", 'wb')
  1568. for ipm in range(self.nPulseMoments):
  1569. # frequency cycing modulation
  1570. #mod = (-1)**(ipm%2) * (-1)**(istack%2)
  1571. buf1 = rawFile.read(4)
  1572. buf2 = rawFile.read(4)
  1573. if hack:
  1574. # hack to do some data analysis
  1575. subFile.write(buf1)
  1576. subFile.write(buf2)
  1577. # end hack
  1578. #N_chan = struct.unpack('>i', rawFile.read(4))[0]
  1579. #N_samp = struct.unpack('>i', rawFile.read(4))[0]
  1580. N_chan = struct.unpack('>i', buf1 )[0]
  1581. N_samp = struct.unpack('>i', buf2 )[0]
  1582. T = N_samp * self.dt
  1583. TIMES = np.arange(0, T, self.dt) - .0002 # small offset in GMR DAQ?
  1584. DATA = np.zeros([N_samp, N_chan+1])
  1585. for ichan in range(N_chan):
  1586. DATADUMP = rawFile.read(4*N_samp)
  1587. #subFile.write(DATADUMP)
  1588. for irec in range(N_samp):
  1589. DATA[irec,ichan] = struct.unpack('>f', DATADUMP[irec*4:irec*4+4])[0]
  1590. if hack:
  1591. # hack to do some data analysis (array survey)
  1592. for ichan in range(5):
  1593. for irec in range(N_samp):
  1594. bdata = struct.pack( '>f', DATA[irec,ichan] )
  1595. subFile.write(bdata)
  1596. for ichan in range(5,6):
  1597. for irec in range(N_samp):
  1598. bdata = struct.pack( '>f', DATA[irec,ichan] + DATA[irec,ichan+2] )
  1599. #bdata = struct.pack( '>f', .23 )
  1600. subFile.write(bdata)
  1601. for ichan in range(6,N_chan):
  1602. for irec in range(N_samp):
  1603. bdata = struct.pack( '>f', DATA[irec,ichan] )
  1604. subFile.write(bdata)
  1605. if plot:
  1606. canvas.ax1.clear()
  1607. canvas.ax2.clear()
  1608. canvas.ax3.clear()
  1609. # Save into Data Cube
  1610. for ichan in chan:
  1611. self.DATADICT["Pulse 1"][ichan][ipm][istack] = DATA[:,eval(ichan)+3][nds:nds+nd1] * invGain
  1612. self.DATADICT["Pulse 1"]["TIMES"] = TIMES[nds:nds+nd1]
  1613. self.DATADICT["Pulse 1"]["CURRENT"][ipm][istack] = DATA[:,1][nps:nps+npul] * invCGain
  1614. self.DATADICT["Pulse 1"]["PULSE_TIMES"] = TIMES[nps:nps+npul]
  1615. if plot:
  1616. #canvas.ax2.plot(self.DATADICT["Pulse 1"]["TIMES"], self.DATADICT["Pulse 1"][ichan][ipm][istack], label="Pulse 1 FID data ch. "+str(ichan)) #, color='blue')
  1617. canvas.ax1.plot(self.DATADICT["Pulse 1"]["PULSE_TIMES"], self.DATADICT["Pulse 1"]["CURRENT"][ipm][istack] , color='black')
  1618. canvas.ax3.plot(self.DATADICT["Pulse 1"]["TIMES"], self.DATADICT["Pulse 1"][ichan][ipm][istack], label="Pulse 1 FID data ch. "+str(ichan)) #, color='blue')
  1619. #canvas.draw()
  1620. # plot reference channels?
  1621. for ichan in rchan:
  1622. self.DATADICT["Pulse 1"][ichan][ipm][istack] = DATA[:,eval(ichan)+3][nds:nds+nd1] * invGain
  1623. self.DATADICT["Pulse 1"]["TIMES"] = TIMES[nds:nds+nd1]
  1624. if plot:
  1625. canvas.ax2.plot(self.DATADICT["Pulse 1"]["TIMES"], self.DATADICT["Pulse 1"][ichan][ipm][istack], label="Pulse 1 FID ref ch. "+str(ichan)) #, color='blue')
  1626. if plot:
  1627. canvas.ax3.legend(prop={'size':6})
  1628. canvas.ax2.legend(prop={'size':6})
  1629. canvas.ax1.set_title("stack "+str(istack)+" pulse index " + str(ipm), fontsize=8)
  1630. canvas.ax1.set_xlabel("time [s]", fontsize=8)
  1631. canvas.ax1.set_ylabel("Current [A]", fontsize=8)
  1632. canvas.ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  1633. canvas.ax2.set_ylabel("RAW signal [V]", fontsize=8)
  1634. canvas.ax2.tick_params(axis='both', which='major', labelsize=8)
  1635. canvas.ax2.tick_params(axis='both', which='minor', labelsize=6)
  1636. canvas.ax2.set_xlabel("time [s]", fontsize=8)
  1637. canvas.ax2.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  1638. canvas.ax3.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  1639. canvas.draw()
  1640. percent = (int) (1e2*((float)((iistack*self.nPulseMoments+ipm+1)) / (len(procStacks)*self.nPulseMoments)))
  1641. #self.progressTrigger.emit(percent)
  1642. self.progressTrigger.emit(percent)
  1643. iistack += 1
  1644. self.enableDSP()
  1645. self.doneTrigger.emit()
  1646. def load4PhaseT1Data(self, base, procStacks, chan, rchan, FIDProc, canvas, deadTime, plot):
  1647. """
  1648. Designed to load GMR 4-phase data which use the following convention for phase cycles
  1649. P1 P2
  1650. Stack 1 -> 0 0 <-- <--
  1651. Stack 2 -> 0 pi/2 | <-- <--
  1652. Stack 3 -> pi/2 0 <-- | <--
  1653. Stack 4 -> pi/2 pi/2 <-- <--
  1654. The cycle is determined by stack indice. Walbrecker proposes for pulse2 data (Stack2 - Stack1) / 2
  1655. equivalently (Stack 4 - Stack3) will yield the same voltage response wrt. the second pulse.
  1656. Alternatively Stack 4 can be converted to be aligned with Stack 1 by negating, and Stack 3 Can be aligned with Stack 2 by negating
  1657. Then there are just the two phase cycles that can be stacked like normal.
  1658. Unfortunately, we need to stack each cycle first, then perform corrections for phase cycling. The reason for this is that otherwise,
  1659. the entire point is lost, as the signal that is desired to be cancelled out may not be balanced evenly across the stacks. That is to say,
  1660. if there is an uneven number of a certain phase cycle.
  1661. We could, I suppose impose this condition, but I think I would rather not?
  1662. + more samples for std. deviation calculation
  1663. + single spikes will have less residual effect
  1664. - can no longer do normality tests etc. and remove data that are suspect.
  1665. - requires a dumb stack, and may also require removal of entire stacks of data
  1666. Additonally, the GMR varies phase as a function of pulse moment index, so that the first pusle moment is zero phase, the second is pi/2 the third is zero ...
  1667. This however, is altered by the above convention. It gets a little complicated...
  1668. """
  1669. import struct
  1670. canvas.reAx2()
  1671. # not in any headers but this has changed, NOT the place to do this. MOVE
  1672. self.prePulseDelay = 0.01 # delay before pulse
  1673. self.deadTime = deadTime # instrument dead time before measurement
  1674. self.samp = 50000. # in case this is a reproc, these might have
  1675. self.dt = 1./self.samp # changed
  1676. invGain = 1./self.RxGain
  1677. invCGain = self.CurrentGain
  1678. #################################################################################
  1679. # figure out key data indices
  1680. # Pulse
  1681. nps = (int)((self.prePulseDelay)*self.samp)
  1682. nps2 = (int)((self.prePulseDelay+self.interpulseDelay)*self.samp)
  1683. npul = (int)(self.pulseLength[0]*self.samp) #+ 100
  1684. np2 = (int)(self.pulseLength[1]*self.samp) #+ 100
  1685. # Data
  1686. nds = nps+npul+(int)((self.deadTime)*self.samp); # indice pulse 1 data starts
  1687. nd1 = (int)((self.interpulseDelay)*self.samp) # samples in first pulse
  1688. nd2s = nps+npul+nd1+(int)((self.deadTime)*self.samp); # indice pulse 2 data starts
  1689. nd2 = (int)((1.)*self.samp) # samples in first pulse
  1690. nd1 -= (int)((.028)*self.samp) + nps # some time to get ready for next pulse
  1691. #################################################################################
  1692. # Data structures
  1693. PULSES = [FIDProc]
  1694. if FIDProc == "Both":
  1695. PULSES = ["Pulse 1","Pulse 2"]
  1696. self.DATADICT = {}
  1697. self.DATADICT["nPulseMoments"] = self.nPulseMoments
  1698. self.DATADICT["stacks"] = procStacks
  1699. self.DATADICT["PULSES"] = PULSES
  1700. for pulse in PULSES:
  1701. self.DATADICT[pulse] = {}
  1702. self.DATADICT[pulse]["chan"] = chan
  1703. self.DATADICT[pulse]["rchan"] = rchan
  1704. self.DATADICT[pulse]["CURRENT"] = {}
  1705. for ichan in np.append(chan,rchan):
  1706. self.DATADICT[pulse][ichan] = {}
  1707. for ipm in range(self.nPulseMoments):
  1708. self.DATADICT[pulse][ichan][ipm] = {}
  1709. self.DATADICT[pulse]["CURRENT"][ipm] = {}
  1710. for istack in procStacks:
  1711. self.DATADICT[pulse][ichan][ipm][istack] = np.zeros(3)
  1712. self.DATADICT[pulse]["CURRENT"][ipm][istack] = np.zeros(3)
  1713. ##############################################
  1714. # Read in binary data
  1715. iistack = 0
  1716. for istack in procStacks:
  1717. rawFile = open(base + "_" + str(istack) + ".lvm", 'rb')
  1718. for ipm in range(self.nPulseMoments):
  1719. N_chan = struct.unpack('>i', rawFile.read(4))[0]
  1720. N_samp = struct.unpack('>i', rawFile.read(4))[0]
  1721. T = N_samp * self.dt
  1722. TIMES = np.arange(0, T, self.dt) - .0002 # small offset in GMR DAQ?
  1723. DATA = np.zeros([N_samp, N_chan+1])
  1724. for ichan in range(N_chan):
  1725. DATADUMP = rawFile.read(4*N_samp)
  1726. for irec in range(N_samp):
  1727. DATA[irec,ichan] = struct.unpack('>f', DATADUMP[irec*4:irec*4+4])[0]
  1728. if plot:
  1729. canvas.ax1.clear()
  1730. canvas.ax2.clear()
  1731. li = np.shape( DATA[:,4][nd2s:nd2s+nd2] )[0]
  1732. ######################################
  1733. # save into DATA cube
  1734. # TODO, changing iFID to 'Pulse 1' or 'Pulse 2'
  1735. for ichan in chan:
  1736. if FIDProc == "Pulse 1":
  1737. self.DATADICT["Pulse 1"][ichan][ipm][istack] = DATA[:,ichan+3][nds:nds+nd1] * invGain
  1738. self.DATADICT["Pulse 1"]["TIMES"] = TIMES[nds:nds+nd1]
  1739. self.DATADICT["Pulse 1"]["CURRENT"][ipm][istack] = DATA[:,1][nps:nps+npul] * invCGain
  1740. self.DATADICT["Pulse 1"]["PULSE_TIMES"] = TIMES[nps:nps+npul]
  1741. if plot:
  1742. canvas.ax1.plot(self.DATADICT["Pulse 1"]["TIMES"], self.DATADICT["Pulse 1"][ichan][ipm][istack], label="Pulse 1 FID data ch. "+str(ichan)) #, color='blue')
  1743. canvas.ax2.plot(self.DATADICT["Pulse 1"]["PULSE_TIMES"], self.DATADICT["Pulse 1"]["CURRENT"][ipm][istack] , color='black')
  1744. elif FIDProc == "Pulse 2":
  1745. self.DATADICT["Pulse 2"][ichan][ipm][istack] = DATA[:,ichan+3][nd2s:nd2s+nd2] *invGain
  1746. self.DATADICT["Pulse 2"]["TIMES"] = TIMES[nd2s:nd2s+nd2]
  1747. self.DATADICT["Pulse 2"]["CURRENT"][ipm][istack] = DATA[:,1][nps2:nps2+np2] * invCGain
  1748. self.DATADICT["Pulse 2"]["PULSE_TIMES"] = TIMES[nps2:nps2+np2]
  1749. if plot:
  1750. canvas.ax1.plot(self.DATADICT["Pulse 2"]["TIMES"], self.DATADICT["Pulse 2"][ichan][ipm][istack], label="Pulse 2 FID data ch. "+str(ichan)) #, color='blue')
  1751. canvas.ax2.plot( self.DATADICT["Pulse 2"]["PULSE_TIMES"], self.DATADICT["Pulse 2"]["CURRENT"][ipm][istack], color='black' )
  1752. else:
  1753. self.DATADICT["Pulse 1"][ichan][ipm][istack] = DATA[:,ichan+3][nds:nds+nd1] * invGain
  1754. self.DATADICT["Pulse 2"][ichan][ipm][istack] = DATA[:,ichan+3][nd2s:nd2s+nd2] * invGain
  1755. self.DATADICT["Pulse 1"]["TIMES"] = TIMES[nds:nds+nd1]
  1756. self.DATADICT["Pulse 2"]["TIMES"] = TIMES[nd2s:nd2s+nd2]
  1757. self.DATADICT["Pulse 1"]["CURRENT"][ipm][istack] = DATA[:,1][nps:nps+npul] * invCGain
  1758. self.DATADICT["Pulse 1"]["PULSE_TIMES"] = TIMES[nps:nps+npul]
  1759. self.DATADICT["Pulse 2"]["CURRENT"][ipm][istack] = DATA[:,1][nps2:nps2+np2] * invCGain
  1760. self.DATADICT["Pulse 2"]["PULSE_TIMES"] = TIMES[nps2:nps2+np2]
  1761. if plot:
  1762. canvas.ax1.plot(self.DATADICT["Pulse 1"]["TIMES"], self.DATADICT["Pulse 1"][ichan][ipm][istack], label="Pulse 1 FID data ch. "+str(ichan)) #, color='blue')
  1763. canvas.ax1.plot(self.DATADICT["Pulse 2"]["TIMES"], self.DATADICT["Pulse 2"][ichan][ipm][istack], label="Pulse 2 FID data ch. "+str(ichan)) #, color='blue')
  1764. canvas.ax2.plot( self.DATADICT["Pulse 1"]["PULSE_TIMES"], self.DATADICT["Pulse 1"]["CURRENT"][ipm][istack] , color='black' )
  1765. canvas.ax2.plot( self.DATADICT["Pulse 2"]["PULSE_TIMES"], self.DATADICT["Pulse 2"]["CURRENT"][ipm][istack] , color='black')
  1766. for ichan in rchan:
  1767. if FIDProc == "Pulse 1":
  1768. self.DATADICT["Pulse 1"][ichan][ipm][istack] = DATA[:,ichan+3][nds:nds+nd1] * invGain
  1769. self.DATADICT["Pulse 1"]["TIMES"] = TIMES[nds:nds+nd1]
  1770. if plot:
  1771. canvas.ax1.plot(self.DATADICT["Pulse 1"]["TIMES"], self.DATADICT["Pulse 1"][ichan][ipm][istack], label="Pulse 1 FID ref ch. "+str(ichan)) #, color='blue')
  1772. elif FIDProc == "Pulse 2":
  1773. self.DATADICT["Pulse 2"][ichan][ipm][istack] = DATA[:,ichan+3][nd2s:nd2s+nd2] * invGain
  1774. self.DATADICT["Pulse 2"]["TIMES"] = TIMES[nd2s:nd2s+nd2]
  1775. if plot:
  1776. canvas.ax1.plot(self.DATADICT["Pulse 2"]["TIMES"], self.DATADICT["Pulse 2"][ichan][ipm][istack], label="Pulse 2 FID ref ch. "+str(ichan)) #, color='blue')
  1777. else:
  1778. self.DATADICT["Pulse 1"][ichan][ipm][istack] = DATA[:,ichan+3][nds:nds+nd1] * invGain
  1779. self.DATADICT["Pulse 2"][ichan][ipm][istack] = DATA[:,ichan+3][nd2s:nd2s+nd2] * invGain
  1780. self.DATADICT["Pulse 1"]["TIMES"] = TIMES[nds:nds+nd1]
  1781. self.DATADICT["Pulse 2"]["TIMES"] = TIMES[nd2s:nd2s+nd2]
  1782. if plot:
  1783. canvas.ax1.plot(self.DATADICT["Pulse 1"]["TIMES"], self.DATADICT["Pulse 1"][ichan][ipm][istack], label="Pulse 1 FID ref ch. "+str(ichan)) #, color='blue')
  1784. canvas.ax1.plot(self.DATADICT["Pulse 2"]["TIMES"], self.DATADICT["Pulse 2"][ichan][ipm][istack], label="Pulse 2 FID ref ch. "+str(ichan)) #, color='blue')
  1785. if plot:
  1786. canvas.ax1.legend(prop={'size':6})
  1787. canvas.ax1.set_title("stack "+str(istack)+" pulse index " + str(ipm), fontsize=8)
  1788. canvas.ax1.set_xlabel("time [s]", fontsize=8)
  1789. canvas.ax1.set_ylabel("RAW signal [V]", fontsize=8)
  1790. canvas.ax2.set_ylabel("Current [A]", fontsize=8)
  1791. canvas.ax2.tick_params(axis='both', which='major', labelsize=8)
  1792. canvas.ax2.tick_params(axis='both', which='minor', labelsize=6)
  1793. canvas.ax2.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  1794. canvas.ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
  1795. canvas.draw()
  1796. # update GUI of where we are
  1797. percent = (int) (1e2*((float)((iistack*self.nPulseMoments+ipm+1)) / (len(procStacks)*self.nPulseMoments)))
  1798. self.progressTrigger.emit(percent)
  1799. iistack += 1
  1800. self.enableDSP()
  1801. self.doneTrigger.emit()
  1802. if __name__ == "__main__":
  1803. if len(sys.argv) < 4:
  1804. print( "mrsurvey path/to/header <stack1> <stackN> ")
  1805. exit()
  1806. GMR = GMRDataProcessor()
  1807. GMR.readHeaderFile(sys.argv[1])
  1808. GMR.Print()
  1809. if GMR.pulseType == "FID":
  1810. GMR.loadFIDData(sys.argv[1], sys.argv[2], sys.argv[3], 5)
  1811. if GMR.pulseType == "4PhaseT1":
  1812. GMR.load4PhaseT1Data(sys.argv[1], sys.argv[2], sys.argv[3], 5)
  1813. pylab.show()