NatComm2020/in-vitro/LongRec.py

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