|
@@ -0,0 +1,164 @@
|
|
|
+# -*- coding: utf-8 -*-
|
|
|
+"""
|
|
|
+Created on Fri Jul 5 09:21:49 2019
|
|
|
+
|
|
|
+@author: rgarcia1
|
|
|
+"""
|
|
|
+
|
|
|
+from PhyREC.NeoInterface import NeoSegment#, ReadMCSFile
|
|
|
+import PhyREC.PlotWaves as Rplt
|
|
|
+import quantities as pq
|
|
|
+import matplotlib.pyplot as plt
|
|
|
+import numpy as np
|
|
|
+import neo
|
|
|
+import PhyREC.SignalProcess as RPro
|
|
|
+import deepdish as dd
|
|
|
+from scipy.signal import hilbert
|
|
|
+
|
|
|
+
|
|
|
+def ReadMCSFile(McsFile,Include, OutSeg=None, SigNamePrefix=''):
|
|
|
+ import McsPy.McsData as McsData
|
|
|
+
|
|
|
+ Dat = McsData.RawData(McsFile)
|
|
|
+ Rec = Dat.recordings[0]
|
|
|
+ NSamps = Rec.duration
|
|
|
+
|
|
|
+ if OutSeg is None:
|
|
|
+ OutSeg = NeoSegment()
|
|
|
+
|
|
|
+ for AnaStrn, AnaStr in Rec.analog_streams.iteritems():
|
|
|
+ if len(AnaStr.channel_infos) == 1:
|
|
|
+ continue
|
|
|
+
|
|
|
+ for Chn, Chinfo in AnaStr.channel_infos.iteritems():
|
|
|
+ print 'Analog Stream ', Chinfo.label, Chinfo.sampling_frequency
|
|
|
+ ChName = str(SigNamePrefix + Chinfo.label)
|
|
|
+ print ChName
|
|
|
+ if ChName not in Include:
|
|
|
+ continue
|
|
|
+
|
|
|
+ Fs = Chinfo.sampling_frequency
|
|
|
+ Var, Unit = AnaStr.get_channel_in_range(Chn, 0, NSamps)
|
|
|
+ sig = neo.AnalogSignal(pq.Quantity(Var, Chinfo.info['Unit']),
|
|
|
+ t_start=0*pq.s,
|
|
|
+ sampling_rate=Fs.magnitude*pq.Hz,
|
|
|
+ name=ChName)
|
|
|
+
|
|
|
+ OutSeg.AddSignal(sig)
|
|
|
+ return OutSeg
|
|
|
+
|
|
|
+def MeanStd(Data):
|
|
|
+ Arr = np.zeros([len(Data.keys()),len(Data[Data.keys()[0]]['psd'])])
|
|
|
+ for iT,TrtName in enumerate(Data.keys()):
|
|
|
+ Arr[iT,:] = Data[TrtName]['psd'][:,0]
|
|
|
+
|
|
|
+ return np.mean(Arr,0), np.std(Arr,0)
|
|
|
+
|
|
|
+
|
|
|
+Path = '23072019/B12784O18-T3/'
|
|
|
+InFileM = Path + '2019-07-23T16-59-52B12784O18-T3-ACDC-PostEth-LowFreqNoise-0.25V.h5' ############
|
|
|
+InFileS = Path + '2019-07-23T16-59-52B12784O18-T3-ACDC-PostEth-LowFreqNoise-0.25V_2.h5' ##########
|
|
|
+
|
|
|
+Gm = dd.io.load(Path + 'GM-B12784O18-T3')
|
|
|
+
|
|
|
+BW = 100
|
|
|
+ivgainDC = 118.8*pq.V#the gain (1e6) is already applied to the saved signal ## Check this gain
|
|
|
+ivgainAC = 1188*pq.V
|
|
|
+TWind = (140*pq.s, 2350*pq.s)
|
|
|
+
|
|
|
+#DCch = ('ME5', 'ME7', 'ME29', 'ME31', 'SE5', 'SE7', 'SE29', 'SE31')
|
|
|
+
|
|
|
+DCch = ('SE7','ME31',)#'ME1','ME2')#'SE7','SE31',)#'ME1','ME2')#'SE7','SE31',)#'SE29','SE31','ME5','ME7','ME29','ME31',)
|
|
|
+
|
|
|
+Include = DCch#('ME1','ME2','ME4','ME6')
|
|
|
+color = ['r','orange','g','b']
|
|
|
+
|
|
|
+
|
|
|
+Rec = ReadMCSFile(InFileM,
|
|
|
+ Include,
|
|
|
+ OutSeg=None,
|
|
|
+ SigNamePrefix='M')
|
|
|
+
|
|
|
+Rec = ReadMCSFile(InFileS,
|
|
|
+ Include,
|
|
|
+ OutSeg=Rec,
|
|
|
+ SigNamePrefix='S')
|
|
|
+
|
|
|
+#%% Plot Sig
|
|
|
+plt.close('all')
|
|
|
+Slots= []
|
|
|
+SlotsComp = []
|
|
|
+counter = -1
|
|
|
+for sig in Rec.Signals():
|
|
|
+ SigProDC = [{'function': RPro.Gain, 'args': {'Gain': pq.A/(Gm[sig.name][8]*(pq.A/pq.V)*ivgainDC)}},
|
|
|
+ {'function': RPro.SetZero, 'args' : {'TWind' :TWind}},
|
|
|
+# {'function': RPro.DownSampling, 'args' : {'Fact' :10}},
|
|
|
+ {'function': RPro.Filter, 'args': {'Type':'bandpass',
|
|
|
+ 'Order':2,
|
|
|
+ 'Freqs':(0.005,0.05)}},
|
|
|
+ ]
|
|
|
+ if sig.name not in DCch:
|
|
|
+ continue
|
|
|
+ counter +=1
|
|
|
+ sig.ProcessChain = SigProDC
|
|
|
+ sig = sig.GetSignal(TWind)
|
|
|
+
|
|
|
+ if counter==0:
|
|
|
+ sig0 = sig.GetSignal(TWind)
|
|
|
+ else:
|
|
|
+ sig1 = sig.GetSignal(TWind)
|
|
|
+
|
|
|
+ Slots.append(Rplt.WaveSlot(sig,
|
|
|
+ Position=0,
|
|
|
+ Alpha=0.5))
|
|
|
+ SlotsComp.append(Rplt.WaveSlot(sig,
|
|
|
+ Position=0,
|
|
|
+ Alpha=0.5))
|
|
|
+Splots = Rplt.PlotSlots(Slots)
|
|
|
+Splots.PlotChannels(Time=TWind,
|
|
|
+ Units='V')
|
|
|
+
|
|
|
+
|
|
|
+fig, (ax1,ax2,ax3) = plt.subplots(3,1)
|
|
|
+HT0 = hilbert(np.array(sig0).T[0,:])
|
|
|
+
|
|
|
+ax1.plot([],color='w', label='[2-4 Hz]')
|
|
|
+ax1.plot(sig0.times, sig0*1e6,'k',label = 'signal')
|
|
|
+ax1.plot(sig0.times,np.real(HT0)*1e6,'b',label = 'Real(HT)')
|
|
|
+ax1.plot(sig0.times,np.abs(HT0)*1e6,'r',label = 'Mag(HT)')
|
|
|
+
|
|
|
+HT1= hilbert(np.array(sig1).T[0,:])
|
|
|
+
|
|
|
+ax2.plot([],color='w', label='ISA [0.007-0.1 Hz]')
|
|
|
+ax2.plot(sig1.times,sig1*1e6,'k',label = 'signal')
|
|
|
+ax2.plot(sig1.times,np.real(HT1)*1e6,'k',label = 'Real(HT)')
|
|
|
+ax2.plot(sig1.times,np.abs(HT1)*1e6,'r',label = 'Mag(HT)')
|
|
|
+ax3.plot(sig1.times,np.angle(HT1,deg=True))
|
|
|
+ax1.set_xlim([140,2350])
|
|
|
+ax2.set_xlim([140,2350])
|
|
|
+ax3.set_xlim([140,2350])
|
|
|
+ax1.set_ylabel('Sig ($\mu$V)')
|
|
|
+ax2.set_ylabel('Sig ($\mu$V)')
|
|
|
+ax3.set_ylabel('Sig phase (degree)')
|
|
|
+
|
|
|
+ax1.legend()
|
|
|
+ax2.legend()
|
|
|
+plt.figure()
|
|
|
+plt.plot(np.angle(HT1,deg=True),np.abs(HT1),'*')
|
|
|
+plt.figure()
|
|
|
+plt.plot(np.angle(HT1,deg=True),np.real(HT1),'*')
|
|
|
+plt.figure()
|
|
|
+plt.plot(np.angle(HT1,deg=True), np.abs(HT0),'*')
|
|
|
+
|
|
|
+
|
|
|
+fig, ax = plt.subplots()
|
|
|
+
|
|
|
+phaseShift = np.angle(HT1,deg=True)-np.angle(HT0,deg=True)
|
|
|
+for iSamp, samp in enumerate(phaseShift):
|
|
|
+ if samp <= -180:
|
|
|
+ phaseShift[iSamp] = samp+360
|
|
|
+ if samp >=180:
|
|
|
+ phaseShift[iSamp] = samp-360
|
|
|
+
|
|
|
+h = plt.hist2d(phaseShift, np.log10(np.abs(HT0)), normed = True, bins=50, range= [[np.min(phaseShift), np.max(phaseShift)], [-6.0, -2.5]],vmin=0, vmax=0.04, cmap='jet')#norm = LogNorm())#
|
|
|
+plt.colorbar(h[3])
|