TVSD/_code/extract_MUA_v2_N.m

%%% extact MUA - THINGS
% v.2-Gemini
% 2024 P. Papale fecit

clear all
addpath(genpath('\_code\code_utils_v2\'));

% constants
monkey = 'monkeyN';
datadir_gen = '\';
temp_folder = '\_temps\';
dates = {'20220111','20220112','20220113','20220114'};
instances = 1:8;
stimbit = 1;
trial_length = .3;
pre_trial = .1;
post_trial = trial_length-pre_trial;
log_prefix = 'THINGS_';
plot_avg = 1;

%
downsample_factor = 30;
Fs = 30000;
Fbp=[500,5000];
N  = 2;    % filter order
Fn = Fs/2; % Nyquist frequency
Fl=200;
% FsD = Fs/downsample_factor;
% FnD = FsD/2; % Downsampled Nyquist frequency
[B, A] = butter(N, [min(Fbp)/Fn max(Fbp)/Fn]); % compute filter coefficients
[D, C] = butter(N,Fl/Fn,'low'); % compute filter coefficients
[F, E] = butter(N,[(50-5)/Fn (50+5)/Fn],'stop'); % compute filter coefficients

for i = 1:5
    d50{i} = designfilt('bandstopiir','FilterOrder',2, ...
        'HalfPowerFrequency1',50*i-5,'HalfPowerFrequency2',50*i+5, ...
        'DesignMethod','butter','SampleRate',Fs);
end

for i = 1:5
    d60{i} = designfilt('bandstopiir','FilterOrder',2, ...
        'HalfPowerFrequency1',60*i-5,'HalfPowerFrequency2',60*i+5, ...
        'DesignMethod','butter','SampleRate',Fs);
end

tb =-pre_trial*1000:1:(post_trial)*1000;
instance_chns = 1:128:128*length(instances);
load('\\vs03\vs03-vandc-1\Latents\Passive_Fixation\monkeyN\_logs\1024chns_mapping_20220105.mat')
restarted = 0;
get_diode = 1;
for day = 1:length(dates)
    clear datadir all_folders blocks
    datadir = [datadir_gen,monkey,'\',dates{day},'\'];
    disp(['%%%%%% DAY: ',num2str(day),' ...']);
    
    all_folders = dir(datadir);
    all_folders = all_folders([all_folders(:).isdir]);
    all_folders = all_folders(~ismember({all_folders(:).name},{'.','..'}));
    blocks = length(all_folders);
    for b = 1:blocks
        block_number = str2num(all_folders(b).name(7:end));
        tic
        disp(['%%% Block: ',num2str(block_number),' ...']);
        
        % load and group the logs
        % MAT files are: [#trial #train_pic #test_pic #pic_rep #ncount #correct]
        % ALLMAT files will be: [#trial #train_pic #test_pic #pic_rep #ncount #day]
        clear MAT temp_MAT
        load([datadir_gen,monkey,'\_logs\',log_prefix,monkey,'_',dates{day},'_B',num2str(block_number),'.mat'],'MAT')
        temp_MAT = MAT(MAT(:,end)>0,1:5);
        temp_MAT = [temp_MAT ones([length(temp_MAT) 1])*day];
        save([datadir_gen,monkey,temp_folder,'TEMP_MAT_',dates{day},'_B',num2str(block_number),'.mat'],'temp_MAT');
        
        for instance_idx = instances
            clear instance_name
            instance_name = ['instance',num2str(instance_idx),'_B',sprintf('%03.0f', block_number)];
            clear partition_idx partition_name
            partition_idx = instance_idx;
            partition_name = ['partition',num2str(partition_idx),'_B',sprintf('%03.0f', block_number)];
            if ~isfile([datadir_gen,monkey,temp_folder,'TEMP_MUA_',dates{day},'_',partition_name,'.mat'])
                if instance_idx == 1 || restarted
                    get_diode = 1;
                else
                    get_diode = 0;
                end
                
                % load the raw data
                clear instance_NS6_filename RAW
                instance_NS6_filename = [datadir,'Block_',num2str(block_number),'\',instance_name,'.ns6'];
                RAW = openNSx(instance_NS6_filename);
                
                if RAW.RawData.PausedFile==1
                    RAW.Data = cell2mat(RAW.Data);
                end
                
                % load the stored events
                clear instance_NEV_filename EVENT Trials_idxs Trials
                instance_NEV_filename = [datadir,'Block_',num2str(block_number),'\',instance_name,'.nev'];
                EVENT = openNEV(instance_NEV_filename,'nosave','noread');
                Trials_idxs=find(EVENT.Data.SerialDigitalIO.UnparsedData==2^stimbit);%starts at 2^0, till 2^7
                Trials=EVENT.Data.SerialDigitalIO.TimeStamp(Trials_idxs);%time stamps corresponding to stimulus onset
                
                assert(length(Trials)==length(MAT),'There are (%d) recorded trials while (%d) trails were logged',length(Trials),length(MAT));
                
                % get the difference between the diode and the stimulus bit
                if get_diode == 1
                    clear diode dio_onset
                    if instance_idx ~=1
                        % load the raw data
                        clear temp_instance_name temp_RAW
                        temp_instance_name = ['instance',num2str(1),'_B',sprintf('%03.0f', block_number)];
                        instance_NS6_filename = [datadir,'Block_',num2str(block_number),'\',temp_instance_name,'.ns6'];
                        temp_RAW = openNSx(instance_NS6_filename);
                        if temp_RAW.RawData.PausedFile==1
                            temp_RAW.Data = cell2mat(temp_RAW.Data);
                        end
                        clear temp_instance_NEV_filename temp_EVENT temp_Trials_idxs temp_Trials
                        temp_instance_NEV_filename = [datadir,'Block_',num2str(block_number),'\',instance_name,'.nev'];
                        temp_EVENT = openNEV(temp_instance_NEV_filename,'nosave','noread');
                        temp_Trials_idxs=find(temp_EVENT.Data.SerialDigitalIO.UnparsedData==2^stimbit);%starts at 2^0, till 2^7
                        temp_Trials=temp_EVENT.Data.SerialDigitalIO.TimeStamp(temp_Trials_idxs);%time stamps corresponding to stimulus onset
                        
                        diode = temp_RAW.Data(143,:);
                        clear temp_RAW
                    else
                        diode = RAW.Data(143,:);
                        temp_Trials = Trials;
                    end
                    for trl=1:length(temp_Trials)
                        clear diode_trl
                        temp_int = temp_Trials(trl):temp_Trials(trl)+Fs*trial_length-1;
                        diode_trl=diode(temp_int);
                        f=find(abs(diode_trl)>1*10^4,1,'first');
                        if isempty(f)
                            dio_onset(trl) = 0;
                        else
                            dio_onset(trl)= f;
                        end
                    end
                    restarted = 0;
                end
                
                clear temp_block_mua
                temp_block_mua = nan([128 length(temp_MAT) length(tb)]);
                for chn = 1:128
                    clear S dum1 dum2
                    S=double(squeeze(RAW.Data(chn,:)))';
                    
                    % bandpass
                    dum1 = FiltFiltM(B, A, S); % apply filter to the data
                    % rectify
                    dum2 = abs(dum1);
                    
                    % low-pass
                    clear muafilt
                    muafilt = FiltFiltM(D, C, dum2);
                    
                    % remove 50Hz line noise (+ harmonics)
                    buttLoop = muafilt;
                    for i = 1:5
%                         d = designfilt('bandstopiir','FilterOrder',2, ...
%                             'HalfPowerFrequency1',50*i-5,'HalfPowerFrequency2',50*i+5, ...
%                             'DesignMethod','butter','SampleRate',Fs);
%                         buttLoop = filtfilt(d,buttLoop);
%                         buttLoop = filtfilt(d50{i},buttLoop);
                        [F, E] = butter(N,[(50*i-5)/Fn (50*i+5)/Fn],'stop');
                        buttLoop = FiltFiltM(F, E, buttLoop);
                    end
                    
                    % remove 60Hz screen refresh rate (+ harmonics)
                    for i = 1:5
%                         d = designfilt('bandstopiir','FilterOrder',2, ...
%                             'HalfPowerFrequency1',60*i-5,'HalfPowerFrequency2',60*i+5, ...
%                             'DesignMethod','butter','SampleRate',Fs);
%                         buttLoop = filtfilt(d,buttLoop);
%                         buttLoop = filtfilt(d60{i},buttLoop);
                        [F, E] = butter(N,[(60*i-5)/Fn (60*i+5)/Fn],'stop');
                        buttLoop = FiltFiltM(F, E, buttLoop);
                    end
                    muafilt = buttLoop';
                    
                    % select only correct trials and downsample
                    clear  temp_chn
                    temp_chn = nan([length(temp_MAT) length(tb)]);
                    z = 1;
                    for trl=1:length(Trials)
                        if MAT(trl,end)>0
                            clear temp_int temp_data
                            temp_int = Trials(trl)+dio_onset(trl)-Fs*pre_trial:Trials(trl)+dio_onset(trl)+Fs*post_trial-1;
                            temp_chn(z,:) = downsample(muafilt(temp_int),downsample_factor);
                            z=z+1;
                        end
                    end
                    temp_block_mua(chn,:,:) = temp_chn;
                end
                if plot_avg
                    plot(squeeze(nanmean(temp_block_mua,[1,2])));
                end
                % store temps
                
                clear mapping_inst
                mapping_inst = mapping(instance_chns(instance_idx):instance_chns(instance_idx)+128-1);
                mapping_inst = mapping_inst-((instance_idx-1)*128);
                temp_block_mua = temp_block_mua(mapping_inst,:,:);

                save([datadir_gen,monkey,temp_folder,'TEMP_MUA_',dates{day},'_',partition_name,'.mat'],'temp_block_mua','-v7.3');
                pause(5)
                disp(['% (Partiton: ',num2str(partition_name),' done!)']);
            else
                restarted = 1;
            end
            disp(['%%% ... done in: ',num2str(round(toc)),'s'])
        end
    end
end