%% Just few definitions of variables and file names
maze_csv = [recording '_M.csv'];
EEG_file = [recording '_N.mat'];
channels=[2:2:12];
%channels=[1:6]
channels_name = {'F3';'F4';'C3';'C4'; 'P3'; 'P4'};


%% This block loads EEG and construct timetable to synchronize with MAze

load (EEG_file)

eeg = [ContinuousData{channels+1,4}];
Fs = ContinuousData{2,2};
eeg_timestamp = [ContinuousData{channels+1,3}];
base = datetime(RecordingDateTime,'InputFormat','MM/dd/uuuu HH:mm:ss');
aux_time = datetime(seconds(eeg_timestamp(:,1)) + base,'InputFormat','MM/dd/uuuu HH:mm:ss');

%EEG preprocessing goes here, by now butterworth filtering, CAR and
%normalization is running. ICA must be also here, but such processing will
%depend on each recording characteristic a no standard could be stablished
[b,a] = butter(4,2*[0.5 40]/Fs);
eeg = filtfilt(b, a, eeg);

eeg = detrend(eeg');
eeg = eeg';
%eeg = normalizasd(eeg);

% If BPts will be constructed such process must be implemented here prior
% to timetable construction.

eeg_TT = array2timetable(eeg, 'VariableNames', channels_name,'RowTimes', aux_time);

%%
% Maze table has several fields, it is read as a table an converted to
% timetable using first column as Time
behavior_TT =  readtable(maze_csv);
aux_time = datetime(behavior_TT.Time + base - timeofday(base));
behavior_TT.Time = aux_time;
behavior_TT = table2timetable(behavior_TT);

%% 
% InObject2 is a column indicating timestamps with animal exploring
% NewObject, only those timestamps are extracted and used to synchronize
% with EEG recordings

% time table are sliced to commonrage, MAZE time table must be first value and EEG ime table last one
combined_TT = synchronize(behavior_TT,eeg_TT,'commonrange');
combined_TT.Time = combined_TT.Time - combined_TT.Time(1);

% HeadDistanceToObjectxBorder determines when animal were exploring
% first values different from NaN are identified and sliced, explorations were considered as 2 seconds differences from last value 
% idx contains indexes from each epoch exploration.

idx1 = ~isnan(combined_TT.HeadDistanceToObject1Border);
aux = find(idx1);
idx1 = [ aux(1); aux(find(diff(seconds(combined_TT.Time(aux))) >2) +1 ) ];

idx2 = ~isnan(combined_TT.HeadDistanceToObject2Border);
aux = find(idx2);
idx2 = [ aux(1); aux(find(diff(seconds(combined_TT.Time(aux))) >2) +1 ) ];

%% EEG is resampled and events are detected 
eeg_resampled =  retime(combined_TT,'regular','spline','SampleRate',Fs);

% EEG is converted to a sinlge matrix to use MATLAB function for signal processing
eeg = eeg_resampled{:,{'F3','F4','C3','C4', 'P3' , 'P4'}};

idx1 = reindexing(idx1,combined_TT,eeg_resampled);
idx2 = reindexing(idx2,combined_TT,eeg_resampled);

% EEG epochs are conformed based on events detection
epochs1 = eeg_epochize(eeg,idx1,Fs,[-1000,3000],0);
epochs2 = eeg_epochize(eeg,idx2,Fs,[-1000,3000],0);
t_epoch = [-1000:3000];

% Spectral estimation of each epoch is based on a pwelch function
[espectro1, n1]=pwelch_mat3d(epochs1,1000,0.9);
[espectro2, n2]=pwelch_mat3d(epochs2,1000,0.9);

% now relative power is estimated
relative1 = relative_power_mat3d (espectro1);
relative2 = relative_power_mat3d (espectro2);

%%
print_exploration(t_epoch(n1),20*log10(relative1),20*log10(relative2),1);
print_qtpdf([recording '_alfa.pdf'])

print_exploration(t_epoch(n1),20*log10(relative1),20*log10(relative2),2);
print_qtpdf([recording '_beta.pdf'])

%%
print_exploration_mean(t_epoch(n1)-500,20*log10(relative1),20*log10(relative2),1,int8(contains("ACQ",recording)));
print_qtpdf([recording '_media_alfa.pdf'])

print_exploration_mean(t_epoch(n1)-500,20*log10(relative1),20*log10(relative2),2,int8(contains("ACQ",recording)));
print_qtpdf([recording '_media_beta.pdf'])


