NatComm2020/in-vitro/B12784O18-T2/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, NeoSignal#, ReadMCSFile
import PhyREC.SignalAnalysis as Ran
import PhyREC.PlotWaves as Rplt
import PhyREC.SignalAnalysis as Ran
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
import os
from scipy import integrate
from scipy.interpolate import interp1d

import PyGFET.AnalyzeData as FETana
import PyGFET.PlotDataClass as FETplt

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 = '../Characterization/28062019/B12142O37-T2/'
#    InFileM = Path + '2019-06-29T18-00-36B12142O37-T2-1mVVgsSweep-ICN2-PostEthx2-Pt.h5' ############
#    InFileS = Path + '2019-06-29T18-00-36B12142O37-T2-1mVVgsSweep-ICN2-PostEthx2-Pt_2.h5' ##########
#    LogFile = Path + 'B12142O37-T2-ACDC-PostEthx2-Pt.txt'
#    IVFigFile = 'Figs/InVivo-B12142O37-T2-2.png'
#    GMFigFile = 'Figs/InVivo-B12142O37-T2-2.png'
#    FoutDCvals = 'InVivo-B12142O37-T2-2.h5'
    
    
#
#    Path = '../Characterization/28062019/B12142O30-T4/'
#    InFileM = Path + '2019-06-29T17-31-17B12142O30-T4-1mVVgsSweep-ICN2-Pt-2.h5' ############
#    InFileS = Path + '2019-06-29T17-31-17B12142O30-T4-1mVVgsSweep-ICN2-Pt-2_2.h5' ##########
#    LogFile = Path + 'B12142O30-T4-PreEth-ACDC-1mV-Pt-2.txt'    


#    Path = '../Characterization/03072019/B12142O37-T4/'
#    InFileM = Path + '2019-07-03T09-23-45B12142O37-T4-ACDC-1mVsine-Pt-PostEthx2.h5' ############
#    InFileS = Path + '2019-07-03T09-23-45B12142O37-T4-ACDC-1mVsine-Pt-PostEthx2_2.h5' ##########
#    LogFile = Path + 'B12142O37-T4-ACDC-Pt-PostEthx2.txt' 
#   
#    Path = '../Characterization/03072019/B12142O30-T9/'
#    InFileM = Path + '2019-07-03T10-43-12B12142O30-T9-ACDC-1mVsine-Pt-PostEth.h5' ############
#    InFileS = Path + '2019-07-03T10-43-12B12142O30-T9-ACDC-1mVsine-Pt-PostEth_2.h5' ##########
#    LogFile = Path + 'B12142o30-T9-ACDC-Pt-PostEth.txt' 
    

    Path = ''
    InFileM = Path + '2019-07-23T16-18-20B12784O18-T2-ACDC-PostEth-1mV_2.h5' ############
    InFileS = Path + '2019-07-23T16-18-20B12784O18-T2-ACDC-PostEth-1mV_2_2.h5' ##########
    LogFile = Path + 'B12784O18-T2-ACDC-PostEth_Cy2.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 #the gain (1e6) is already applied to the saved signal  ## Check this gain
    ivgainAC = 1188*pq.V 

#    
#    SigProAC = [{'function': RPro.Gain, 'args': {'Gain': pq.A/ACgain}},
#                {'function': RPro.Filter, 'args': {'Type':'highpass',
#                                                   'Order':2,
#                                                   'Freqs':(1)}},
#                ]
#    
#    SigProDC = [{'function': RPro.Gain, 'args': {'Gain': pq.A/DCgain}},
#                {'function': RPro.Filter, 'args': {'Type':'highpass',
#                                                   'Order':2,
#                                                   'Freqs':(1)}},
#                ]


    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 = []
    SlotsAC = []
    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))

    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
#%% 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)
        
        DCVals = {'Ids': np.array((Ids[sl.name],)).transpose(),
                  'Vds': np.array((Vds,)),
                  'Vgs': np.array(Vgs),
                  'ChName': sl.name,
                  'Name': sl.name,
                  'DateTime': LogVals['Date/Time'][0]}
        DevDCVals[sl.name] = DCVals
   
    FETana.CheckIsOK(DevDCVals, RdsRange=[400, 40e3])
    FETana.CalcGM(DevDCVals)
    pltDC = FETplt.PyFETPlot()
    pltDC.AddAxes(('Ids', 'Gm', 'Rds'))
    pltDC.PlotDataCh(DevDCVals, PltIsOK=True)
    pltDC.AddLegend()
#    if os.path.exists(IVFigFile[0:-7]) == False:
#        os.mkdir((IVFigFile[0:-7]))
#    pltDC.Fig.savefig(IVFigFile[0:-7]+'/'+IVFigFile[5:])  
    
#%%  Calc GM
    GM = {}
    Irms = {}
    Urms = {}
    UrmsDrift = {}
    SNR = {}
    Irms2 = {}
    figgm, axgm = plt.subplots()
    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
#        if sig.name in DCch:
#            sig.SignalProcess = SigProDC
#        else:
#            sig.SignalProcess = SigProAC
#        
        GM[sig.name] = []
        Irms[sig.name] = []
        Irms2[sig.name] = []
        Urms[sig.name] = []
        SNR[sig.name] = []
        UrmsDrift[sig.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 = 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*2)
            Irms[sig.name] = Irms[sig.name]
            
            irms2 = np.sqrt(psd[2]*Fpsd[2]*np.log(100))
            Irms2[sig.name] = np.append(Irms2[sig.name],irms2)
            Irms2[sig.name] = Irms2[sig.name]
            
            
            SNR[sig.name] = 20*np.log10(GM[sig.name]*(0.707/1000)/Irms[sig.name])
#        if float(SNR[sig.name][0]) <=15:
#            continue
            Urms[sig.name] = Irms[sig.name]/GM[sig.name]
        Polate = interp1d(Vgs-0.34,Urms[sig.name])
        VgsInt=np.linspace(-0.2,-0.01,10)
        UrmsDrift[sig.name] = Polate(VgsInt)
        plt.figure(8)
        plt.plot(Vgs-0.34, Urms[sig.name],'k',alpha=0.1)
#        axgm.plot(Vgs, GM[sig.name], label=sig.name) 
            
        plt.figure(9)
        plt.plot(Vgs, SNR[sig.name])
        axgm.plot(Vgs, GM[sig.name], label=sig.name) 
        
        
#        plt.figure(8)
#        plt.plot(Vgs-0.7, Urms[sig.name])
        
        plt.figure(5)
        plt.xlabel('Ugs(V)')
        plt.ylabel('Gm(S)')
#        fig.savefig(IVFigFile[0:-7]+'/'+sig.name+'psd'+'.png')
#        plt.close(fig)
#        plt.close()
    plt.figure(6)
    GMmean, GMstd = MeanStdGM(GM)
    plt.plot(Vgs, GMmean*1000/0.1,'k',label = '1 metal layer')
    plt.fill_between(Vgs-0.70, 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)

    IrmsMean2, IrmsStd2 = MeanStdGM(Irms2)
    plt.semilogy(Vgs, IrmsMean2,'b',label = 'a param')
    plt.fill_between(Vgs, IrmsMean2-IrmsStd2, IrmsMean2+IrmsStd2 ,color = 'b',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-0.34, UrmsMean,'k')
    plt.fill_between(Vgs-0.34, UrmsMean-UrmsStd, UrmsMean+UrmsStd ,color = 'k',alpha =0.3)
    plt.xlabel('V$_{gs}$ - V$_{CNP}$ (V)')
    plt.ylabel('V$_{gs-rms}$ ($\mu$V)')
    plt.legend()
    
    plt.figure(9)
    SNRMean, SNRStd = MeanStdGM(SNR)
    plt.plot(Vgs-0.70, SNRMean,'k',label = '1 metal layer')
    plt.fill_between(Vgs-0.70, SNRMean-SNRStd, SNRMean+SNRStd ,color = 'k',alpha =0.3)
    plt.xlabel('V$_{gs}$ - V$_{CNP}$ (V)')
    plt.ylabel('SNR (dB)')
    plt.legend()
    
    fig, ax = plt.subplots()
    UrmsDriftArray = np.zeros((len(UrmsDrift.keys()),len(UrmsDrift[UrmsDrift.keys()[0]])))
    for ikeyTrt,keyTrt in enumerate(UrmsDrift.keys()):
        for ikeyVgs in np.arange(len(UrmsDrift[UrmsDrift.keys()[0]])): 
            UrmsDriftArray[ikeyTrt, ikeyVgs] = UrmsDrift[keyTrt][ikeyVgs]*1e6
#            ax.plot(np.random.normal(VgsInt[ikeyVgs], 0.003)-0.0015, UrmsDrift[keyTrt][ikeyVgs]*1e6,'*k')
    bpr = ax.boxplot(UrmsDriftArray, positions=VgsInt, widths=0.006)
    
#    axgm.legend()
#    figgm.savefig(IVFigFile[0:-7]+'/'+'Gm'+'.png')
#    dd.io.save(FoutDCvals, (DevDCVals,GM ), ('zlib', 1))
    
    dd.io.save(Path+'GM-B12784O18-T2', GM)

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


for Trt in SNR.keys():
    ch = MCSMapI[Trt]
#    if Trt in DCch:
#        continue
    
    A[MCSMapFacingDown[ch][1],MCSMapFacingDown[ch][0]] = SNR[Trt][7]

plt.imshow(A, interpolation='nearest', vmin=10, vmax=38)
plt.grid(True)
cbar=plt.colorbar()
plt.xlabel('column')
plt.ylabel('row')
cbar.set_label('SNR [dB]', rotation=270, labelpad=15)

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

import matplotlib.colors as colors
for Trt in SNR.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=3, vmax=30, norm=colors.LogNorm(vmin=3, vmax=30))
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)