123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229 |
- # -*- coding: utf-8 -*-
- """
- Created on Fri Jul 5 09:21:49 2019
- @author: rgarcia1
- """
- from PhyREC.NeoInterface import NeoSegment#, ReadMCSFile
- import PhyREC.SignalAnalysis as Ran
- 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
- import csv
- from datetime import datetime
- def ReadMCSFile(McsFile, 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
- 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 ReadLogFile(File):
- Fin = open(File)
-
- reader = csv.reader(Fin, delimiter='\t')
-
- LogVals = {}
- ValsPos = {}
- for il, e in enumerate(reader):
- if il == 0:
- for ih, v in enumerate(e):
- ValsPos[ih] = v
- LogVals[v] = []
- else:
- for ih, v in enumerate(e):
- par = ValsPos[ih]
- if (par=='Vgs') or (par=='Vds') or (par=='Vref'):
- LogVals[par].append(float(v.replace(',','.')))
- elif par == 'Date/Time':
- LogVals[par].append(datetime.strptime(v, '%d/%m/%Y %H:%M:%S'))
- else:
- LogVals[par].append(v)
-
- deltas = np.array(LogVals['Date/Time'])[:]-LogVals['Date/Time'][0]
- LogVals['Time'] = []
- for d in deltas:
- LogVals['Time'].append(d.total_seconds())
-
- Fin.close()
-
- return LogVals
- 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
- ivgainAC = 1188*pq.V
- TWind = (140*pq.s, 2350*pq.s)
- DCch = ('ME5', 'ME7', 'ME29', 'ME31', 'SE5', 'SE7', 'SE29', 'SE31')
- Rec = ReadMCSFile(InFileM,
- OutSeg=None,
- SigNamePrefix='M')
- Rec = ReadMCSFile(InFileS,
- OutSeg=Rec,
- SigNamePrefix='S')
- #plt.close('all')
- #%% Plot Sig High Freq
- plt.close('all')
- SlotsAC= []
- for sig in Rec.Signals():
- SigProDC = [{'function': RPro.Gain, 'args': {'Gain': pq.A/(Gm[sig.name][8]*(pq.A/pq.V)*ivgainAC)}},
- {'function': RPro.SetZero, 'args' : {'TWind' :TWind}},
- {'function': RPro.Filter, 'args': {'Type':'bandpass',
- 'Order':2,
- 'Freqs':(20,200)}},
- ]
- if sig.name in DCch+('SE30',):
- continue
- sig.ProcessChain = SigProDC
- sig = sig.GetSignal(TWind)
- SlotsAC.append(Rplt.WaveSlot(sig,
- Position=0,
- Alpha=0.5))
- Splots = Rplt.PlotSlots(SlotsAC)
- Splots.PlotChannels(Time=TWind,
- Units='V')
- plt.xlabel('time (s)')
- plt.ylabel('Sig (uV)')
- #%% Plot PSD
- fig, AxPs = plt.subplots()
- PSD = Ran.PlotPSD((s for s in SlotsAC),
- Time = TWind,
- Ax=AxPs,
- FMin=1,
- scaling='density',
- Units = 'V',
- label = '20-200Hz',
- color='orange')
- MeanPSDhigh, StdPSDhigh = MeanStd(PSD)
- FpsdHigh = PSD[PSD.keys()[0]]['ff']
- #%%Plot hystogram
- HistData = NeoSegment()
- plt.figure(5)
- signal = np.array([])
- for iSl, sl in enumerate(SlotsAC):
- signal = np.append(signal,np.array(sl.GetSignal(TWind))[:,0]*1000000)
- hist, vect = np.histogram(signal, bins=1000 , range = (np.min(signal),np.max(signal)))
- plt.hist(signal, bins =1000, density = True, range = (np.min(signal),np.max(signal)),alpha=0.3,color='orange',label = '1-100Hz')
- plt.legend()
- #%% Plot Sig
- Slots= []
- for sig in Rec.Signals():
- SigProAC = [{'function': RPro.Gain, 'args': {'Gain': pq.A/(Gm[sig.name][8]*(pq.A/pq.V)*ivgainDC)}},
- {'function': RPro.SetZero, 'args' : {'TWind' :TWind}},
- {'function': RPro.Filter, 'args': {'Type':'bandpass',
- 'Order':2,
- 'Freqs':(0.05,0.5)}},
- ]
- if sig.name not in DCch:
- continue
- sig.ProcessChain = SigProAC
- sig = sig.GetSignal(TWind)
- Slots.append(Rplt.WaveSlot(sig,
- Position=0,
- Alpha=0.5))
- Splots = Rplt.PlotSlots(Slots)
- Splots.PlotChannels(Time=TWind,
- Units='V')
- plt.xlabel('time (s)')
- plt.ylabel('Sig (uV)')
- #%% Plot PSD
- #fig, (AxPs, Axt) = plt.subplots(2,1)
- PSD = Ran.PlotPSD((s for s in Slots),
- Time = TWind,
- Ax=AxPs,
- FMin=0.01,
- scaling='density',
- Units = 'V',
- label = '0.05-0.5Hz',
- color= 'k')
- MeanPSDlow, StdPSDlow = MeanStd(PSD)
- FpsdLow = PSD[PSD.keys()[0]]['ff']
- #%%Plot hystogram
- HistData = NeoSegment()
- plt.figure(5)
- signal = np.array([])
- for iSl, sl in enumerate(Slots):
- signal = np.append(signal,np.array(sl.GetSignal(TWind))[:,0]*1000000)
- hist, vect = np.histogram(signal, bins=1000 , range = (np.min(signal),np.max(signal)))
- plt.hist(signal, bins =1000, density = True, range = (np.min(signal),np.max(signal)),alpha=0.3,color='k', label = '0.01-1Hz')
- plt.xlabel('Signal amplitude ($\mu$V)')
- plt.ylabel('Probability density')
- #%%Plot PSD average
- plt.figure()
- plt.loglog(FpsdLow,MeanPSDlow,'k')
- plt.fill_between(FpsdLow,MeanPSDlow+StdPSDlow,MeanPSDlow-StdPSDlow, color ='k',alpha=0.3)
- plt.loglog(FpsdHigh,MeanPSDhigh,'orange')
- plt.fill_between(FpsdHigh,MeanPSDhigh+StdPSDhigh,MeanPSDhigh-StdPSDhigh, color ='orange',alpha=0.3)
|