NatComm2020/in-vitro/CalcGM_Noise_LMU.py

# -*- coding: utf-8 -*-
"""
Created on Fri Oct 05 15:49:46 2018

@author: aemdlabs
"""

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Wed Sep  5 12:05:21 2018

@author: aguimera
"""

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 csv
from datetime import datetime
from scipy import integrate


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 MeanStd(Data, var):
    Arr = np.zeros([len(Data.keys()),len(Data[Data.keys()[0]][var])])
    for iT,TrtName in enumerate(Data.keys()):
        Arr[iT,:] = Data[TrtName][var][:,0]
    
    return np.mean(Arr,0), np.std(Arr,0)

def MeanStdGM(Data):
    Arr = np.zeros([len(Data.keys()),len(Data[Data.keys()[0]])])
    for iT,TrtName in enumerate(Data.keys()):
        Arr[iT,:] = Data[TrtName]
    
    return np.mean(Arr,0), np.std(Arr,0)

def Integrate(PSD, Freqs, Fmin, Fmax):
    indices = np.where((Freqs >= Fmin) & (Freqs<=Fmax))
    print( Freqs[indices])
    Irms = np.sqrt(integrate.trapz(PSD[indices], Freqs[indices]))      
    return Irms 


MCSMapI={'SE1':'Ch03',
         'SE2':'Ch05',
         'SE3':'Ch01',
         'SE4':'Ch02',
         'SE5':'Ch22',
         'SE6':'Ch06',
         'SE7':'Ch16',
         'SE8':'Ch37',
         'SE9':'Ch20',
         'SE10':'Ch10',
         'SE11':'Ch24',
         'SE12':'Ch08',
         'SE13':'Ch14',
         'SE14':'Ch04',
         'SE15':'Ch18',
         'SE16':'Ch33',
         'SE17':'Ch34',
         'SE18':'Ch60',
         'SE19':'Ch38',
         'SE20':'Ch64',
         'SE21':'Ch40',
         'SE22':'Ch56',
         'SE23':'Ch42',
         'SE24':'Ch70',
         'SE25':'Ch66',
         'SE26':'Ch65',
         'SE27':'Ch68',
         'SE28':'Ch67',
         'SE29':'Ch55',
         'SE30':'Ch62',
         'SE31':'Ch58',
         'SE32':'Ch69',
         'ME1':'Ch57',
         'ME2':'Ch61',
         'ME3':'Ch53',
         'ME4':'Ch63',
         'ME5':'Ch52',
         'ME6':'Ch41',
         'ME7':'Ch49',
         'ME8':'Ch51',
         'ME9':'Ch46',
         'ME10':'Ch45',
         'ME11':'Ch44',
         'ME12':'Ch39',
         'ME13':'Ch54',
         'ME14':'Ch43',
         'ME15':'Ch50',
         'ME16':'Ch47',
         'ME17':'Ch32',
         'ME18':'Ch27',
         'ME19':'Ch30',
         'ME20':'Ch29',
         'ME21':'Ch28',
         'ME22':'Ch25',
         'ME23':'Ch26',
         'ME24':'Ch07',
         'ME25':'Ch21',
         'ME26':'Ch11',
         'ME27':'Ch17',
         'ME28':'Ch15',
         'ME29':'Ch13',
         'ME30':'Ch31',
         'ME31':'Ch19',
         'ME32':'Ch09'}

                          #Col, Row  
MCSMapFacingDown={'Ch58':(0,1),
                  'Ch57':(0,2),
                  'Ch56':(0,3),
                  'Ch55':(0,4),
                  'Ch54':(0,5),
                  'Ch53':(0,6),
                  'Ch52':(0,7),
                  'Ch51':(0,8),
                  'Ch50':(0,9),
                  'Ch49':(0,10),
                  'Ch60':(1,0),
                  'Ch61':(1,1),
                  'Ch62':(1,2),
                  'Ch63':(1,3),
                  'Ch64':(1,4),
                  'Ch65':(1,5),
                  'Ch43':(1,6),
                  'Ch44':(1,7),
                  'Ch45':(1,8),
                  'Ch46':(1,9),
                  'Ch47':(1,10),
                  'Ch70':(2,0),
                  'Ch69':(2,1),
                  'Ch68':(2,2),
                  'Ch67':(2,3),
                  'Ch66':(2,4),
                  'Ch42':(2,5),
                  'Ch41':(2,6),
                  'Ch40':(2,7),
                  'Ch39':(2,8),
                  'Ch38':(2,9),
                  'Ch37':(2,10),
                  'Ch01':(3,0),
                  'Ch02':(3,1),
                  'Ch03':(3,2),
                  'Ch04':(3,3),
                  'Ch05':(3,4),
                  'Ch06':(3,5),
                  'Ch30':(3,6),
                  'Ch31':(3,7),
                  'Ch32':(3,8),
                  'Ch33':(3,9),
                  'Ch34':(3,10),
                  'Ch11':(4,0),
                  'Ch10':(4,1),
                  'Ch09':(4,2),
                  'Ch08':(4,3),
                  'Ch07':(4,4),
                  'Ch29':(4,5),
                  'Ch28':(4,6),
                  'Ch27':(4,7),
                  'Ch26':(4,8),
                  'Ch25':(4,9),
                  'Ch24':(4,10),
                  'Ch12':None,
                  'Ch59':None,
                  'Ch13':(5,1),
                  'Ch14':(5,2),
                  'Ch15':(5,3),
                  'Ch16':(5,4),
                  'Ch17':(5,5),
                  'Ch18':(5,6),
                  'Ch19':(5,7),
                  'Ch20':(5,8),
                  'Ch21':(5,9),
                  'Ch22':(5,10)}

if __name__ == '__main__':


    Path = '23072019/B12784O18-T3/'
    InFileM = Path + '2019-07-23T18-55-56B12784O18-T3-ACDC-PostEth-PostLong.h5' ############
    InFileS = Path + '2019-07-23T18-55-56B12784O18-T3-ACDC-PostEth-PostLong_2.h5' ##########
    LogFile = Path + 'B12784O18-T3-ACDC-PostEth-PostLong.txt'
    
    StartCycle = -1

    LogVals = ReadLogFile(LogFile)
    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 = 'SE31'
    TrigThres = 5e-5
    Vgs = np.array(LogVals['Vgs'])
    Vds = LogVals['Vds'][0]
    
    ivgain1 = 12e3*pq.V  
    ivgain2 = 101
    ACgain = 10*1
    DCgain = 1*pq.V
    Fsig = 10
    StabTime = 10*pq.s
    GuardTime = 1*pq.s
    BW = 100
    ivgainDC = 118.8*pq.V
    ivgainAC = 1188*pq.V 


    Rec = ReadMCSFile(InFileM,
                      OutSeg=None,
                      SigNamePrefix='M')

    Rec = ReadMCSFile(InFileS,
                      OutSeg=Rec,
                      SigNamePrefix='S')

# %%
    plt.close('all')  
    plt.ion() 
    
    SwTimes = GetSwitchTimes(Sig=Rec.GetSignal(TrigChannel),
                             Thres=TrigThres,
                             Plot=True)
   
    SlotsDC = []
    for sig in Rec.Signals():
        if sig.name not in DCch:
            continue

        if sig.name.startswith('M'):
            col = 'r'
        else:
            col = 'g'
        SlotsDC.append(Rplt.WaveSlot(sig,
                                     Position=0,
                                     Color=col,
                                     Alpha=0.5))
    
    SwTimes = LogVals['Time']*pq.s+SwTimes[0]+Delay
#%% calc IV  DC
    DevDCVals = {}
    
    fig, Axt = plt.subplots()
    Ids = {}
    for sl in SlotsDC:
        Ids[sl.name] = []
        for isw, (t, vg) in enumerate(zip(SwTimes, Vgs)):
            ts = SwTimes[isw]+delta
            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
            Ids[sl.name].append(ids)
        

    
#%%  Calc GM
    GM = {}
    Irms = {}
    Urms = {}

    fig, (AxPsd, Axt) = plt.subplots(2,1)     
    fig2, (AxPs, Axt) = plt.subplots(2,1)  
    for sig in Rec.Signals():
        if sig.name[0:3] == 'SEn':
            continue

#        
        GM[sig.name] = []
        Irms[sig.name] = []
        Urms[sig.name] = []

        for isw, (t, vg) in enumerate(zip(SwTimes, Vgs)):
            if isw == len(SwTimes)-1:
                ts = SwTimes[isw]-Delay
            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)     

            PS = Ran.PlotPSD((s,),
                              Time = TWind,
                              Ax=AxPs,
                              FMin=1,
                              Label=str(vg),
                              scaling='spectrum')
            
            ps = PS[sig.name]['psd']
            Fps = PS[sig.name]['ff']
           
            indicesPeak = np.where( ((Fps >= Fsig-4) & (Fps<=Fsig+4)))   

            IDSpeak = np.sqrt(ps[np.argmax(ps[indicesPeak])+indicesPeak[0][0]]+
                             ps[np.argmax(ps[indicesPeak])+indicesPeak[0][0]+1]+
                             ps[np.argmax(ps[indicesPeak])+indicesPeak[0][0]-1])
            
            gm = IDSpeak*1000/0.707
            GM[sig.name] = np.append(GM[sig.name],gm)
            
            PSD = Ran.PlotPSD((s,),
                  Time = TWind,
                  Ax=AxPsd,
                  FMin=1,
                  Label=str(vg),
                  scaling='density')

            psd = PSD[sig.name]['psd'][:,0]
            Fpsd = PSD[sig.name]['ff']
            
#            irms = Integrate(psd, Fpsd, 1.9, 1.9*3.2)
#            Irms[sig.name] = np.append(Irms[sig.name],irms*np.sqrt(2))
            irms = Integrate(psd, Fpsd, 63, 200)
            Irms[sig.name] = np.append(Irms[sig.name],irms*np.sqrt(2))
            
            Irms[sig.name] = Irms[sig.name]
            
            Urms[sig.name] = Irms[sig.name]/GM[sig.name]


    
    plt.figure(12)
    GMmean, GMstd = MeanStdGM(GM)
    plt.plot(Vgs, GMmean*1000/0.1,'k',label = '1 metal layer')
    plt.fill_between(Vgs, GMmean*1000/0.1-GMstd*1000/0.1, GMmean*1000/0.1+GMstd*1000/0.1,color = 'k',alpha =0.3)
    plt.xlabel('V$_{gs}$ - V$_{CNP}$ (V)')
    plt.ylabel('G$_m$ (mS/V)')
    plt.legend()
    
    plt.figure(7)
    IrmsMean, IrmsStd = MeanStdGM(Irms)
    plt.semilogy(Vgs, IrmsMean,'k',label = 'rms')
    plt.fill_between(Vgs, IrmsMean-IrmsStd, IrmsMean+IrmsStd ,color = 'k',alpha =0.3)

    plt.xlabel('V$_{gs}$ - V$_{CNP}$ (V)')
    plt.ylabel('I$_{rms}$ (A)')
    plt.legend()
    
    plt.figure(8)
    UrmsMean, UrmsStd = MeanStdGM(Urms)
    plt.semilogy(Vgs, UrmsMean,'k',label = '1 metal layer')
    plt.fill_between(Vgs, UrmsMean-UrmsStd, UrmsMean+UrmsStd ,color = 'k',alpha =0.3)
    plt.xlabel('V$_{gs}$ - V$_{CNP}$ (V)')
    plt.ylabel('U$_{rms}$ (A)')
    plt.legend()
    


#%% plot map Urms
    
plt.figure()
A=np.log10(np.ones((11,6))*5e-17)

import matplotlib.colors as colors
for Trt in Urms.keys():
    ch = MCSMapI[Trt]
#    if Trt in DCch:
#        continue
    
    A[MCSMapFacingDown[ch][1],MCSMapFacingDown[ch][0]] = (Urms[Trt][9])*1e6

plt.imshow(A, interpolation='nearest', vmin=1, vmax=15, norm=colors.LogNorm(vmin=1, vmax=15))
plt.grid(True)
cbar=plt.colorbar()
plt.xlabel('column',fontsize=12)
plt.ylabel('row',fontsize=12)
cbar.set_label('U$_{gs-rms}$ ($\mu$V)', rotation=270, labelpad=15,fontsize=13)