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