KIAPBCI/kiap_bci/analytics/lfp_analysis1.m

clear all;

file_name1 = '/home/vlachos/data_kiap/20190321-140130/20190321-140130-002.ns2';
file_name2 = '/home/vlachos/data_kiap/20190321-140130/20190321-140130-002.nev';


%% defining parameters


pars.ch_nr = 2;

pars.lfp.filter_order = 4;
pars.lfp.filter_fc_low = 30;
pars.lfp.filter_fc_high = 60;
pars.lfp.fs = 1000;                 % sampling rate

pars.artifact_thr = 1000;

% psth parameters
pars.psth_win =[-1000, 4999];

%% LOAD FILES

ns2 = openNSx('read', 'report',file_name1, 'p:int16');
nev = openNEV('read', 'report',file_name2);

%% CAR
lfp = double(ns2.Data)';
lfp = lfp - mean(lfp,2);
disp(size(lfp))

%% FILTER

[b, a] = butter(pars.lfp.filter_order, [pars.lfp.filter_fc_low, pars.lfp.filter_fc_high] / (pars.lfp.fs / 2));
% [b, a] = butter(pars.lfp.filter_order, [100] / (pars.lfp.fs / 2),'low');
% lfp_filt = filtfilt(b, a, lfp(:,1:pars.ch_nr));
lfp_filt = filtfilt(b, a, lfp(:,30:31));
disp(size(lfp_filt))

%% COMPUTE FFT

Y_lfp = fft(lfp_filt);
L = length(Y_lfp);
P2 = abs(Y_lfp/L);
P1 = P2(1:L/2+1,:);
P1(2:end-1,:) = 2*P1(2:end-1,:);
ff = pars.lfp.fs*(0:(L/2))/L;

% clf
% for ii=1:128
%     plot(ff,P1(:,ii))
%     disp(ii)
%     input('press enter')
% end

%% GET YES/NO TRIGGERS FROM NEV FILE

triggers_ts = nev.Data.Comments.TimeStamp;
triggers_txt = nev.Data.Comments.Text;

triggers_yes = [];
triggers_no = [];

for ii=1:length(triggers_txt)
    if contains(triggers_txt(ii,:), 'no, response, start')
        triggers_no = [triggers_no, triggers_ts(ii)];
    elseif contains(triggers_txt(ii,:), 'yes, response, start')
        triggers_yes = [triggers_yes, triggers_ts(ii)];
    end
end

% get timestamps for ns2
triggers_yes_ns2 = triggers_yes / 30;
triggers_no_ns2 = triggers_no / 30;

%% COMPUTE PSTH

psth_win = pars.psth_win;
psth_len = diff(psth_win)+1;
psth_yes = zeros(length(triggers_yes_ns2),  psth_len, pars.ch_nr);
psth_no = zeros(length(triggers_no_ns2), psth_len, pars.ch_nr);

for ii=1:length(triggers_yes_ns2)
    if triggers_yes_ns2(ii)+psth_win(1)>0 && triggers_yes_ns2(ii)+psth_win(2)<length(lfp_filt)
        tmp = lfp_filt((triggers_yes_ns2(ii)+psth_win(1)):(triggers_yes_ns2(ii)+psth_win(2)),:);
        if ~any(any(tmp>pars.artifact_thr))
            psth_yes(ii,:,:) = tmp;
        else 
            disp(['Trial excluded: ',num2str(ii)])
        end
    end
end

for ii=1:length(triggers_no_ns2)
    if triggers_no_ns2(ii)-psth_win(1)>0 && triggers_no_ns2(ii)+psth_win(2)<length(lfp_filt)
        tmp = lfp_filt((triggers_no_ns2(ii)+psth_win(1)):(triggers_no_ns2(ii)+psth_win(2)),:);
        if ~any(any(tmp>pars.artifact_thr))
            psth_no(ii,:,:) = tmp;
        else
            disp(['Trial excluded: ',num2str(ii)])
        end
    end
end

%% FFT OF PSTH

Y_psth_yes = fft(psth_yes,[],2);
L = length(Y_psth_yes);
P2 = abs(Y_psth_yes/L);
P1 = P2(:,1:L/2+1,:);
P1(:,2:end-1,:) = 2*P1(:,2:end-1,:);
ff = pars.lfp.fs*(0:(L/2))/L;


%% PLOT RESULTS
subplot(2,2,1)
plot(squeeze(mean(psth_no,1)));
subplot(2,2,3)
plot(squeeze(mean(psth_yes,1)),'b');

subplot(2,2,2)
plot(ff, squeeze(mean(P1,1)),'r');
xlim([0,100])

subplot(2,2,4)
plot(ff, squeeze(mean(P1,1)),'b');
xlim([0,100])