NatComm2020/In-vivo/Calibration and longitudinal assessment/Day1-2/Analysis scripts/PhaseAmplitude_Rec2.py

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

@author: rgarcia1
"""

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
import deepdish as dd
from matplotlib.pyplot import cohere
from scipy.signal import hilbert
from matplotlib.colors import LogNorm
import matplotlib.colors as mcolors
from matplotlib.pyplot import cohere

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

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 = '../Raw Data/'
InFileM = Path + '2019-07-31T18-00-05B12784O18-T2-Longterm-Rec2.h5' ############
InFileS = Path + '2019-07-31T18-00-05B12784O18-T2-Longterm-Rec2_2.h5' ##########

Gm = dd.io.load('../CalibrationFile/Cal-Gm20190731_033840.h5')

BW = 100
ivgainDC = 118.8*pq.V
ivgainAC = 1188*pq.V
TWind = (12850*pq.s,15450*pq.s)


DCch = ('SE7','ME31',)

Include = DCch
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')

TWind = (None,20800*pq.s)

Slots= []
SlotsComp = []
counter = -1
for sig in Rec.Signals():
    SigProDC = [{'function': RPro.Gain, 'args': {'Gain': pq.A/(Gm['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)}},
            ]

    SigProRMS = [{'function':RPro.sliding_window, 'args': {'func':RPro.rms,'timewidth': 300*pq.s}}]
    
    if sig.name not in DCch:
        continue
    counter +=1
    sig.ProcessChain = SigProDC
    sig = sig.GetSignal(TWind)
    sigRMS = sig.GetSignal(TWind)
    sigRMS.ProcessChain = SigProRMS
    sigRMS = sigRMS.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')    


HT0 = hilbert(np.array(sig0).T[0,:])

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


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)), bins=50, normed = True, range= [[np.min(phaseShift), np.max(phaseShift)], [-6, -2.5]],vmin=0, vmax=0.006,cmap='jet')#norm = LogNorm())#
plt.colorbar(h[3])