Surface NMR processing and inversion GUI
Nevar pievienot vairāk kā 25 tēmas Tēmai ir jāsākas ar burtu vai ciparu, tā var saturēt domu zīmes ('-') un var būt līdz 35 simboliem gara.

mrsurvey.py 109KB

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