NatComm2020/in-vitro/CalcGM_Noise_LMU_multipleFiles.py

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

@author: aemdlabs
"""


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 deepdish as dd
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__':

    
    Files = [
             'B12142O30-T4/2019-06-29T17-07-11B12142O30-T4-1mVVgsSweep-ICN2-2',
             'B12142O37-T2/2019-06-29T18-00-36B12142O37-T2-1mVVgsSweep-ICN2-PostEthx2-Pt',
             'B12142O37-T5/2019-07-03T09-55-56B12142O37-T5-ACDC-1mVsine-Pt-PostEth',
             'B12142O37-T6/2019-07-03T10-18-36B12142O37-T6-ACDC-1mVsine-Pt-PostEth',
             'B12142O30-T6/2019-07-03T08-31-00B12142O30-T6-ACDC-1mVsine-Pt-PostEth',
             'B12142O30-T9/2019-07-03T10-43-12B12142O30-T9-ACDC-1mVsine-Pt-PostEth',
             'B12784O18-T1/2019-07-05T08-04-13B12784O18-T1-1mVVgsSweep-PostEth-ICN2',
             'B12784O18-T2/2019-07-23T16-18-20B12784O18-T2-ACDC-PostEth-1mV_2',
             'B12784O18-T3/2019-07-23T18-55-56B12784O18-T3-ACDC-PostEth-PostLong',]############
    LogFiles = [
                'B12142O30-T4/B12142O30-T4-PreEth-ACDC1mV-AgAgCl-2.txt',
                'B12142O37-T2/B12142O37-T2-ACDC-PostEthx2-Pt.txt',
                'B12142O37-T5/B12142O37-T5-ACDC-Pt-PostEth.txt',
                'B12142O37-T6/B12142O37-T6-ACDC-Pt-PostEth.txt',
                'B12142O30-T6/B12142O30-T6_Pt_PostEth-ACDC.txt',
                'B12142O30-T9/B12142o30-T9-ACDC-Pt-PostEth.txt',
                'B12784O18-T1/B12784O18-T1-ACdc-postTH.txt',
                'B12784O18-T2/B12784O18-T2-ACDC-PostEth_Cy2.txt',
                'B12784O18-T3/B12784O18-T3-ACDC-PostEth-PostLong.txt']

    UrmsTot = {}
    IrmsTot = {}
    GmTot = {}
    IdsTot = {}
    VGS = {}
    for iFile, File in enumerate(Files):
        StartCycle = -1

        LogFile = LogFiles[iFile]
        InFileM = File + '.h5'
        InFileS = File + '_2.h5'
        
        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 = 'SE29'
        TrigThres = 5e-5
        Vgs = np.array(LogVals['Vgs'])
        Vds = LogVals['Vds'][0]
        VGS[File] = Vgs
        
        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 
    

    
    
        Rec = ReadMCSFile(InFileM,
                          OutSeg=None,
                          SigNamePrefix='M')
    
        Rec = ReadMCSFile(InFileS,
                          OutSeg=Rec,
                          SigNamePrefix='S')
    
    # %% Detect Vgs switch time and arrange DC signals into SlotsDC 
        plt.close('all')  
        plt.ion() 
        
        SwTimes = GetSwitchTimes(Sig=Rec.GetSignal(TrigChannel),
                                 Thres=TrigThres,
                                 Plot=False)
       
        SlotsDC = []
        for sig in Rec.Signals():
            if sig.name not in DCch:
                continue
    
            SlotsDC.append(Rplt.WaveSlot(sig,))
    
        
        SwTimes = LogVals['Time']*pq.s+SwTimes[0]+Delay
        
    #%% calc IV  DC
        DevDCVals = {}
        
        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')

                vio = np.mean(s).magnitude
                ids = (vio*101-(-vg+Vds))/12e3 #apply hardware gain calibration
                Ids[sl.name].append(ids)
            
        IdsTot[File.split('/')[1]] = Ids    


        
    #%%  Calc GM and integrated noise Irms
    
        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': #discard encoder channels
                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 = 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)
    
                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.3)
                Irms[sig.name] = np.append(Irms[sig.name],irms*2) #
                # square root of 2 adjusts the 1/f integrated noise to a 1 order of magnitude frequency range
                Irms[sig.name] = Irms[sig.name]
                Urms[sig.name] = Irms[sig.name]/GM[sig.name]
    
            plt.close('all')   #close psd figures created

        
        UrmsTot[File.split('/')[1]] = Urms
        IrmsTot[File.split('/')[1]] = Irms
        GmTot[File.split('/')[1]] = GM
        
    
    plt.figure()
    GMmean, GMstd = MeanStdGM(GM)
    plt.plot(Vgs-0.70, 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()
    IrmsMean, IrmsStd = MeanStdGM(Irms)
    plt.semilogy(Vgs-0.70, IrmsMean,'k',label = 'rms')
    plt.fill_between(Vgs-0.70, IrmsMean-IrmsStd, IrmsMean+IrmsStd ,color = 'k',alpha =0.3)


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


#    dd.io.save('GmIrmsUrmsIds10Probe_2.h5',(GmTot, IrmsTot, UrmsTot, IdsTot))

    #%% 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]

        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)
    
    #%% plot map Gm
        
    plt.figure()
    A=np.log10(np.ones((11,6))*5e-17)
    
    import matplotlib.colors as colors
    for Trt in Urms.keys():
        ch = MCSMapI[Trt]
        
        A[MCSMapFacingDown[ch][1], MCSMapFacingDown[ch][0]] = (GM[Trt][9]*1000/0.1)
    
    plt.imshow(A, interpolation='nearest', vmin=1, vmax=3, norm=colors.LogNorm(vmin=1, vmax=3))
    plt.grid(True)
    cbar=plt.colorbar()
    plt.xlabel('column',fontsize=12)
    plt.ylabel('row',fontsize=12)
    cbar.set_label('G$_{m}$ (mS/V)', rotation=270, labelpad=15,fontsize=13)