# -*- coding: utf-8 -*- """ Created on Sat Oct 13 18:50:52 2018 @author: aemdlabs """ from PhyREC.NeoInterface import NeoSegment, NeoSignal#, 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 deepdish as dd import csv from datetime import datetime import PyGFET.DataClass as Gfetdat 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 GetSwitchTimes(Sig, Thres=-1e-4, Plot=True): s = Sig.GetSignal(None) ds = np.abs(np.diff(np.array(s), axis=0)) if Plot: plt.figure() plt.plot(s.times, s) plt.plot(s.times[1:], ds) ds = Sig.duplicate_with_new_array(signal=ds) Times = Ran.threshold_detection(ds, threshold=Thres, RelaxTime=5*pq.s) return Times def CalcVgeffAC(Sig, TACchar, Vgs, vgs): i = list(np.around(Vgs,3)).index(vgs) gm = TACchar[i]*pq.A/pq.V return neo.AnalogSignal(Sig/gm, units='V', t_start=Sig.t_start, sampling_rate=Sig.sampling_rate, name=str(Sig.name), file_origin=Sig.file_origin) if __name__ == '__main__': #### Test autosweep Path = '../../Day0/IV charact/' PBSInFileM = Path + '2019-07-30T15-39-02B12784O18-T3-ACDC-PreImplant-NoMemb.h5' ############ PBSInFileS = Path + '2019-07-30T15-39-02B12784O18-T3-ACDC-PreImplant-NoMemb_2.h5' ########## PBSLogFile = Path + 'B12784O18-T3-ACDC-PBS-PreImplant-20s-NoMemb.txt' StartCycle = 0 Directory = '../IV charact/' File = '20190731_033840 Current vs Vgs_5sStep.txt' Fin = open(Directory+File) LogVals = ReadLogFile(PBSLogFile) delta = np.mean([t-LogVals['Time'][it] for it, t in enumerate(LogVals['Time'][1:])])*pq.s Delay = delta * StartCycle DCch = ('ME5', 'ME7', 'ME29', 'ME31', 'SE5', 'SE7', 'SE29', 'SE31') TrigChannel = 'ME29' TrigThres = 1e-5 Vgs = np.array(LogVals['Vgs']) VgsInt = np.linspace(Vgs[0],Vgs[-1],len(Vgs)*10) Vds = [LogVals['Vds'][0]] ivgain1 = 12e3*pq.V ivgain2 = 101 ACgain = 10*1/pq.V DCgain = 1*pq.V Fsig = 10 StabTime = 2*pq.s GuardTime = 1*pq.s BW = 100 ivgainDC = 118.8*pq.V ivgainAC = 1188*pq.V Rec = ReadMCSFile(PBSInFileM, OutSeg=None, SigNamePrefix='M') Rec = ReadMCSFile(PBSInFileS, OutSeg=Rec, SigNamePrefix='S') # %% In vivo Ids plt.close('all') ChKo = {} reader = csv.reader(Fin, delimiter='\t') lenVgs = 1000 Vgs = np.zeros(lenVgs) Ids = np.zeros([lenVgs,8]) reader = csv.reader(Fin, delimiter='\t') for il, e in enumerate(reader): if il == 0: continue print(e) vg = e[0] vg = vg.replace(',','.') Vgs[il-1] = float(vg) for iids, ids in enumerate(e[1:]): ids = ids.replace(',','.') Ids[il-1,iids] = float(ids) color = ['m','r','orange','y','g','c','b','violet'] plt.figure() Ud0 = [] for i in range(8): Vgs = np.append(Vgs[0:2],np.trim_zeros(Vgs[2:])) ids = np.trim_zeros(Ids[:,i]) plt.plot(Vgs,ids,color = color[i]) if np.any(ids>=5.0): index = np.where(ids == np.amin(ids)) print(ids, np.amin(ids), index) Ud0.append(Vgs[index[0][0]]) Ud0meanInVivo = np.mean(np.array(Ud0)) #%% SwTimes = GetSwitchTimes(Sig=Rec.GetSignal(TrigChannel), Thres=TrigThres, Plot=True) SlotsDC = [] SlotsAC = [] for sig in Rec.Signals(): sig.__class__ = NeoSignal if sig.name in DCch: if sig.name.startswith('M'): col = 'r' else: col = 'g' SlotsDC.append(Rplt.WaveSlot(sig, Position=0, Color=col, Alpha=0.5)) if sig.name.startswith('M'): col = 'k' else: col = 'b' SlotsAC.append(Rplt.WaveSlot(sig, Position=1, Color=col, Alpha=0.5)) Splots = Rplt.PlotSlots(SlotsDC) Splots.PlotChannels(Time=None, Units='mV') Splots.PlotEvents(SwTimes, color='k') Splots.PlotEvents(LogVals['Time']*pq.s+SwTimes[0]+Delay, color='r') SwTimes = LogVals['Time']*pq.s+SwTimes[0]+Delay figId, Axiv = plt.subplots() fig, Axt = plt.subplots() Ids = {} for sl in SlotsDC: Ids[sl.name] = [] for isw, (t, vg) in enumerate(zip(SwTimes, Vgs)): if isw == len(SwTimes)-1: ts = sl.Signal.t_stop else: ts = SwTimes[isw+1] TWind = (t+StabTime, ts-GuardTime) s = sl.GetSignal(TWind, Units='V') Axt.plot(s) vio = np.mean(s).magnitude ids = (vio*101-(-vg+Vds))/12e3 if ids <= 7e-6: ChKo[sl.name] = True else: ChKo[sl.name] = False Ids[sl.name].append(ids) if ChKo[sl.name] == False: Axiv.plot(Vgs,Ids[sl.name], label=sl.name) DevACvals = {} Ud0mean = 0 dev = Ids Ud0 = np.zeros(len(dev)) for ich, (ch, ids) in enumerate(dev.iteritems()): if ChKo[ch] == False: print(ich) DevACvals[ch] = { 'Ids': np.array(ids[:-1]), 'Vgs': Vgs[:-1], 'Vds': Vds} DataDC = Gfetdat.DataCharDC(DevACvals[ch]) Ud0[ich] = DataDC.GetUd0() print(Ud0[ich]) Ud0mean=Ud0mean+Ud0[ich] Axiv.plot(VgsInt, DataDC.GetIds(VgsInt),'--') Ud0meanPBS = Ud0mean/(ich+1) #%% Calc GM GM = {} CalFile = {} figgm, axgm = plt.subplots() for sig in Rec.Signals(): GM[sig.name] = [] fig, (AxPsd, Axt) = plt.subplots(2,1) for isw, (t, vg) in enumerate(zip(SwTimes, Vgs)): if isw == len(SwTimes)-1: ts = sl.Signal.t_stop else: ts = SwTimes[isw+1] TWind = (t+StabTime, ts-GuardTime) s = sig.GetSignal(TWind, Units ='V') if s.name in DCch: s = (s*ivgain2-(-vg+Vds)*pq.V)/ivgain1 else: s = s/(ivgain1*ACgain/ivgain2) Axt.plot(s.times, s, label=vg, alpha=0.5) PSD = Ran.PlotPSD((s,), Time = TWind, Ax=AxPsd, FMin=1, Label=str(vg), scaling='spectrum') psd = PSD[sig.name]['psd'] Fpsd = PSD[sig.name]['ff'] indicesPeak = np.where( ((Fpsd >= Fsig-4) & (Fpsd<=Fsig+4))) IDSpeak = np.sqrt(psd[np.argmax(psd[indicesPeak])+indicesPeak[0][0]]+ psd[np.argmax(psd[indicesPeak])+indicesPeak[0][0]+1]+ psd[np.argmax(psd[indicesPeak])+indicesPeak[0][0]-1]) gm = IDSpeak*1000/0.707 GM[sig.name]=np.append(GM[sig.name],gm) if np.any(GM[sig.name]>= 2e-5) : ChKo[sig.name] = False color = 'k' else: ChKo[sig.name] = True color = 'r' if sig.name in DCch: color = 'g' axgm.plot(Vgs, GM[sig.name], label=sig.name, color = color) plt.close(fig) axgm.legend() CalFile = {'GM': GM, 'Ud0PBS':Ud0meanPBS, 'Ud0InVivo':Ud0meanInVivo, 'DeltaUd0':Ud0meanPBS-Ud0meanInVivo, 'Vgs': Vgs-(Ud0meanPBS-Ud0meanInVivo)} CalFile['ChKo'] = ChKo dd.io.save('../Cal Data/'+'Cal-Gm'+File[:15]+'.h5', CalFile)