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

mrsurvey.py 131KB

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