Surface NMR processing and inversion GUI
Nelze vybrat více než 25 témat Téma musí začínat písmenem nebo číslem, může obsahovat pomlčky („-“) a může být dlouhé až 35 znaků.

mrsurvey.py 152KB

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