NatComm2020/in-vitro/PhaseAmpCouplingNoise.py

# -*- coding: utf-8 -*-
"""
Created on Fri Jul  5 09:21:49 2019

@author: rgarcia1
"""

from PhyREC.NeoInterface import NeoSegment#, ReadMCSFile
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
from scipy.signal import hilbert


def ReadMCSFile(McsFile,Include, 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
            if ChName not in Include:
                continue
            
            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 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#the gain (1e6) is already applied to the saved signal  ## Check this gain
ivgainAC = 1188*pq.V
TWind = (140*pq.s, 2350*pq.s)

#DCch = ('ME5', 'ME7', 'ME29', 'ME31', 'SE5', 'SE7', 'SE29', 'SE31')

DCch = ('SE7','ME31',)#'ME1','ME2')#'SE7','SE31',)#'ME1','ME2')#'SE7','SE31',)#'SE29','SE31','ME5','ME7','ME29','ME31',)

Include = DCch#('ME1','ME2','ME4','ME6')
color = ['r','orange','g','b']


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

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

#%% Plot Sig
plt.close('all')
Slots= []
SlotsComp = []
counter = -1
for sig in Rec.Signals():
    SigProDC = [{'function': RPro.Gain, 'args': {'Gain': pq.A/(Gm[sig.name][8]*(pq.A/pq.V)*ivgainDC)}},
            {'function': RPro.SetZero, 'args' : {'TWind' :TWind}},
#            {'function': RPro.DownSampling, 'args' : {'Fact' :10}},
            {'function': RPro.Filter, 'args': {'Type':'bandpass',
                                               'Order':2,
                                               'Freqs':(0.005,0.05)}},
            ]
    if sig.name not in DCch:
        continue
    counter +=1
    sig.ProcessChain = SigProDC
    sig = sig.GetSignal(TWind)
    
    if counter==0:
        sig0 = sig.GetSignal(TWind)
    else:
        sig1 = sig.GetSignal(TWind)
        
    Slots.append(Rplt.WaveSlot(sig,
                                 Position=0,
                                 Alpha=0.5))
    SlotsComp.append(Rplt.WaveSlot(sig,
                                 Position=0,
                                 Alpha=0.5))
Splots = Rplt.PlotSlots(Slots)
Splots.PlotChannels(Time=TWind,
                    Units='V')    


fig, (ax1,ax2,ax3) = plt.subplots(3,1)
HT0 = hilbert(np.array(sig0).T[0,:])

ax1.plot([],color='w', label='[2-4 Hz]')
ax1.plot(sig0.times, sig0*1e6,'k',label = 'signal')
ax1.plot(sig0.times,np.real(HT0)*1e6,'b',label = 'Real(HT)')
ax1.plot(sig0.times,np.abs(HT0)*1e6,'r',label = 'Mag(HT)')

HT1= hilbert(np.array(sig1).T[0,:])

ax2.plot([],color='w', label='ISA [0.007-0.1 Hz]')
ax2.plot(sig1.times,sig1*1e6,'k',label = 'signal')
ax2.plot(sig1.times,np.real(HT1)*1e6,'k',label = 'Real(HT)')
ax2.plot(sig1.times,np.abs(HT1)*1e6,'r',label = 'Mag(HT)')
ax3.plot(sig1.times,np.angle(HT1,deg=True))
ax1.set_xlim([140,2350])
ax2.set_xlim([140,2350])
ax3.set_xlim([140,2350])
ax1.set_ylabel('Sig ($\mu$V)')
ax2.set_ylabel('Sig ($\mu$V)')
ax3.set_ylabel('Sig phase (degree)')

ax1.legend()
ax2.legend()
plt.figure()
plt.plot(np.angle(HT1,deg=True),np.abs(HT1),'*')
plt.figure()
plt.plot(np.angle(HT1,deg=True),np.real(HT1),'*')
plt.figure()
plt.plot(np.angle(HT1,deg=True), np.abs(HT0),'*')


fig, ax = plt.subplots()

phaseShift = np.angle(HT1,deg=True)-np.angle(HT0,deg=True)
for iSamp, samp in enumerate(phaseShift):
    if samp <= -180:
        phaseShift[iSamp] = samp+360 
    if samp >=180:
        phaseShift[iSamp] = samp-360 
        
h = plt.hist2d(phaseShift, np.log10(np.abs(HT0)), normed = True, bins=50, range= [[np.min(phaseShift), np.max(phaseShift)], [-6.0, -2.5]],vmin=0, vmax=0.04, cmap='jet')#norm = LogNorm())#
plt.colorbar(h[3])