123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164 |
- # -*- 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])
|