KIAPBCI/kiap_bci/analytics/lfp_analysis1_2.py

import numpy as np
from scipy import signal
import aux
import importlib
import matplotlib.pyplot as plt
import sys
from analytics import lfp_analytics as lfpan
importlib.reload(lfpan)


interactive_off = int(sys.argv[1])
if interactive_off:
    plt.ioff()
else:
    plt.ion()

# plt.close('all')

params = aux.load_config()
# path = '/home/vlachos/devel/vg/kiapvmdev/data/clinical/neural_new/nsx/'
path = '/media/vlachos/bck_disk1/kiap/recordings/20190705/20190705-143437/'
# file_name0 = '20190607-123832-001.ns2'     # motor mapping 1
file_name0 = '20190705-143437-001.ns2'     # motor mapping 1


file_name = path + file_name0    # motor session 5
# params.lfp.array1 = np.delete(params.lfp.array1, params.lfp.array1_exclude)
# params.lfp.array2 = np.delete(params.lfp.array2, params.lfp.array2_exclude)
psth_win = params.lfp.psth_win

# file_name = '/home/vlachos/data_kiap/20190321-135311/20190321-135311-001'
# path = '/media/kiap/kiap_backup/Recordings/K01/Recordings/20190326-160445/'
# path = '/kiap/data/clinical/neural_new/nsx/'
# file_name = path + '20190326-160445-001.ns2'

array1 = params.lfp.array1
array2 = params.lfp.array2

color = ['C0', 'C1', 'C2', 'C3', 'C4']
args = {'arr_nr': 2, 'erase': True, 'log': True, 'xmin': 0.01, 'xmax': 20, 'ymin': -1e2, 'ymax': 5e3,
    'save_fig': False}

arr_id = 1

if 'an' not in locals():
    an = lfpan.lfp_analytics(params)
an.read_data(file_name)
an.preprocess()
an.get_triggers(speller='feedback')
an.compute_psth(arr_nr=arr_id)
an.calc_stft_psth(arr_nr=arr_id, nperseg=2000)
an.plot_psth_stft(arr_nr=arr_id, ch_ids=range(32), fmax=200, save_fig=True)  
xx
an.plot_channels(ch_ids=[], arr_id=args['arr_nr'], save_fig=args['save_fig'], i_stop=-1)
for ii in range(5):

    markerline, stemlines, baseline = plt.stem(an.triggers_mapping[ii]/1000, np.arange(an.triggers_mapping[ii].size) * 0+10, color[ii])
    plt.setp(stemlines, 'linewidth', 3)
xx
an.plot_psth_time(ch_ids=params.lfp.array21, save_fig=True)
xx

an.calc_stft(arr_nr=2, nperseg=2000)
an.plot_stft(arr_nr=2, save_fig=True)
# an.calc_spectra(arr_nr=args['arr_nr'], nperseg=20000, axis=0, force=True, onesided=True, detrend='constant')
# an.plot_spectra_bands(**args)
# an.calc_cc_time(save_fig=args['save_fig'])

sys.exit(1)
# an.get_triggers(n_triggers=20)
# trials_msk = get_valid_trials(triggers_no_ns2, triggers_yes_ns2)
# data1[data1 > params.lfp.artifact_thr] = 0
# data1[data1 < -params.lfp.artifact_thr] = 0
# data0 = data0 - np.mean(data0, 1)[:, np.newaxis]
# reader.nev_data['Comments'][0]

# ff, tt, S = signal.spectrogram(an.data2[:, :, 0], params.lfp.fs, axis=0, nperseg=nperseg, noverlap=int(nperseg * 0.9)) 
fidx = (ff > 0) & (ff < 2)
plt.figure(1)
plt.clf()
# plt.plot(ff[fidx], np.log10(S[fidx, 4, :]), 'C1')
# plt.pcolormesh(tt, ff[fidx], np.log10(S[fidx, 4, :]))
plt.plot(tt, np.sum(S[fidx, :, :], axis=0))
plt.draw()
plt.show()
# nperseg = 20000
# ff2, Y2 = signal.welch(an.data02, fs=an.params.lfp.fs, axis=0, nperseg=nperseg, detrend=False, scaling='spectrum', return_onesided=False)
# from scipy.fftpack import fftshift
# Y2 = fftshift(Y2)
# ff2 = fftshift(ff2)
plt.xlim(-10,10)




# tidx = tt < 1000
plt.figure(1)
plt.clf()
plt.subplot(121)
band_id = 2
# fidx = (an.ff2 >= 0) & (an.ff2 < 20)
fidx = (an.ff2 >= 0) & (an.ff2 < 200)
plt.plot(an.ff2[fidx], np.log10(an.Y2[fidx, :, band_id]), 'C0')
plt.plot(an.ff2[fidx], np.mean(np.log10(an.Y2[fidx, :, band_id]), axis=1), 'k')
plt.plot(an.ff2[fidx], np.median(np.log10(an.Y2[fidx, :, band_id]), axis=1), 'k--')
plt.ylim(-2, 4)
plt.xlabel('Hz')
plt.ylabel('log spectrum V**2')
plt.title(f'Sensor mapping, 20190606-214850-001.ns2')
plt.draw()
plt.show()

# plt.subplot(122)
# plt.plot(ff2[fidx], np.log10(Y2[fidx, :, 0]), 'C1')
# plt.plot(ff2[fidx], np.mean(np.log10(Y2[fidx, :, 0]), axis=1), 'k')
# plt.plot(ff2[fidx], np.median(np.log10(Y2[fidx, :, 0]), axis=1), 'k--')
# plt.ylim(-2, 4)
# plt.xlabel('Hz')
# plt.ylabel('log spectrum V**2')
# plt.title(f'Motor mapping, 20190607-134213-001.ns2')
# plt.plot(an.ff2[fidx], np.log10(an.Y2[fidx, :, 0].mean(1)), 'k')
# plt.pcolormesh(tt, ff[fidx], np.log10(S[fidx, :, :].mean(1)))
# plt.plot(ff, np.log10(Y[:, :, 2]), 'C0')
# plt.pcolormesh(ff, range(Y.shape[1]), np.log10(Y[:, :, 2].T))
# plt.pcolormesh(ff, tt, np.log10(S[:, 4, :].T))
# plt.colorbar()
# plt.xlim(0,20)
# plt.xlim(, 10)
plt.draw()
plt.show()



xxx
# COMPUTE PSTH
# ----------------------------------

# psth shape: (trials, time, channels)
offset = 600 * params.lfp.fs  # offset in sec to get data for normalization
triggers_yes_norm = np.hstack((triggers_yes_ns2 + offset, triggers_yes_ns2 + offset + 10))
triggers_no_norm = np.hstack((triggers_yes_ns2 + offset, triggers_yes_ns2 + offset + 10))
if np.any(triggers_yes_norm > data01.shape[1]):
    log.error('triggers_norm outside imported data')
psth_no, psth_yes = compute_psth(triggers_yes_ns2, triggers_no_ns2, arr_nr=1, sub_band=params.lfp.sub_band)
psth_norm1, psth_norm2 = compute_psth(triggers_yes_norm, triggers_no_norm, arr_nr=1, sub_band=params.lfp.sub_band)


# COMPUTE PSTH SPECTROGRAMS
# ----------------------------------

nperseg = params.lfp.spectra.spgr_len
nperseg = min(nperseg, psth_no.shape[1])
# ff, tt, Sxx1 = signal.spectrogram(psth_no.mean(0), params.lfp.fs, axis=0, nperseg=nperseg, noverlap=int(nperseg * 0.9))
# ff, tt, Sxx2 = signal.spectrogram(psth_yes.mean(0), params.lfp.fs, axis=0, nperseg=nperseg, noverlap=int(nperseg * 0.9))

# Sxx1.shape = (trials, ff, channels, time)
# Sxx1_norm.shape = (ff,time, channels)
ff, tt, Sxx1 = signal.spectrogram(psth_no, params.lfp.fs, axis=1, nperseg=nperseg, noverlap=int(nperseg * 0.9))
ff, tt, Sxx2 = signal.spectrogram(psth_yes, params.lfp.fs, axis=1, nperseg=nperseg, noverlap=int(nperseg * 0.9))

ff, tt, Sxx1_norm = signal.spectrogram(psth_norm1.mean(0), params.lfp.fs, axis=0, nperseg=nperseg, noverlap=int(nperseg * 0.9))
ff, tt, Sxx2_norm = signal.spectrogram(psth_norm2.mean(0), params.lfp.fs, axis=0, nperseg=nperseg, noverlap=int(nperseg * 0.9))


# Sxx1 = Sxx1.mean(1)
# Sxx2 = Sxx2.mean(1)
Sxx1_norm = Sxx1_norm.mean(2)
Sxx2_norm = Sxx2_norm.mean(2)

tt = tt + psth_win[0] / params.lfp.fs


t_idx = tt < 1

Sxx1_norm = Sxx1_norm[:, np.newaxis].repeat(len(tt), 1)
Sxx2_norm = Sxx2_norm[:, np.newaxis].repeat(len(tt), 1)


# if params.lfp.normalize:
    # Sxx1 = Sxx1 / Sxx1_norm
    # Sxx2 = Sxx2 / Sxx2_norm

# 11. PLOT RESULTS
# ----------------------------------

if params.lfp.plot.general:
    # plot_spectra()
    plot_psth(ff=ff, sub_band=params.lfp.sub_band, ch_id=0, trial_id=0, step=100)
    pass

# xxx
plt.ioff()

# for trial_id in range(4, 5):
#     for ch_id in range(ch_nr1 - 1, ch_nr1):
        
#         plot_psth(ff=ff, sub_band=0, ch_id=ch_id, trial_id=0, step=100)
#         fig2 = plt.gcf()
#         fname = path + f'results/trial_{trial_id}_ch_{ch_id}.png'
#         # fig2.savefig(fname)
#         print(ch_id)