123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201 |
- 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)
|