Surface NMR processing and inversion GUI
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

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