Browse Source

Subir archivos a 'in-vitro'

Ragmon Garcia Cortadella 1 year ago
parent
commit
3f872a4fe2
1 changed files with 164 additions and 0 deletions
  1. 164 0
      in-vitro/PhaseAmpCouplingNoise.py

+ 164 - 0
in-vitro/PhaseAmpCouplingNoise.py

@@ -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])