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

mrsurvey.py 153KB

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