%% NeuroDOT Workshop Pre-Processing Script
% This script will be used during the hands-on portion of the
% Pre-Processing section of the workshop.

% The dataset used in this script is ND_LUMO_Data_Sample_MOT001.mat, which
% can be downloaded from the shared folder. 

% This dataset is from a motor finger tappping experiment where a 
% participant was instructed to tap their fingers in response to an 'L' or
% % an 'R' in the center of a screen. The stimulus consisted of 12 34 second blocks, 
% where the participant tapped either their left or right hand for 10 seconds, 
% % followed by 24 seconds of rest.


%% PREPROCESSING PIPELINE
load(['ND_LUMO_Data_Sample_MOT001.mat']); % data, info, flags

flags.RsdMax = 35;
flags.AbsDiff = 0.5;
flags.FloorTh = 10e-8;
flags.DtKeep = 1;
dt=34;  % Block length


%% Crop data to match A Matrix' measurement list
% Load measurement list from A Matrix
temp = load(['A_LUMO_51_tile_on_participant_mesh_r3dlt50mm.mat'],'info'); % Sensitivity matrix

% Crop data
mA=cat(2,temp.info.pairs.Src,temp.info.pairs.Det,temp.info.pairs.WL);
mD=cat(2,info.pairs.Src,info.pairs.Det,info.pairs.WL);
keepInfo=ismember(mD,mA,'rows');
data=data(keepInfo,:);                                % Choose measurements from data
info.pairs=temp.info.pairs;                           % Grab relevant measurements from pairs with A
info.MEAS.Clipped = info.MEAS.Clipped(keepInfo);


%% General Data Quality Assessment 
Plot_RawData_Time_Traces_Overview(data,info,flags);   % Time traces
info = Plot_RawData_Metrics_II_DQC(data,info,flags);  % Spectrum, falloff, and good signal metrics
info = Plot_RawData_Cap_DQC(data,info,flags);         % Cap-relevant views


%% View data before filtering
lmdata_b4_filt = logmean(data);                                           % Logmean Light Levels
info_b4_filt  = FindGoodMeas(lmdata_b4_filt , info, 0.075, flags);        % Detect Noisy Channels
lmdata_b4_filt  = detrend_tts(lmdata_b4_filt );                           % Detrend Data
info_b4_filt.GVTD = CalcGVTD(lmdata_b4_filt(info_b4_filt.MEAS.GI & info_b4_filt.pairs.r2d<20,:));         % Calculate GVTD

% Measurements to include
keep = info_b4_filt.pairs.WL==2 & info_b4_filt.pairs.r2d < 35;

% Visualize
figure('Position',[100 100 550 780])
subplot(3,1,1); semilogy(data(keep,:)'); 
set(gca,'XLimSpec','tight'), xlabel('Time (samples)'), 
ylabel('\phi') 
m=max(max(abs(lmdata_b4_filt(keep,:))));
subplot(3,1,2); imagesc(lmdata_b4_filt(keep,:),[-1,1].*m); 
colorbar('Location','northoutside');
xlabel('Time (samples)');ylabel('Measurement #')
[ftmag0,ftdomain] = fft_tts(lmdata_b4_filt(keep,:),info_b4_filt.system.framerate); % Generate average spectrum
ftmag = rms(ftmag0,1);
subplot(3,1,3); semilogx(ftdomain,ftmag);
xlabel('Frequency (Hz)');ylabel('|X(f)|');xlim([1e-3 10])

nlrGrayPlots_180818(lmdata_b4_filt,info_b4_filt); % Gray Plot with synch points


%% Look at block-averaged data without preprocessing
badata_b4_filt = BlockAverage(lmdata_b4_filt, info_b4_filt.paradigm.synchpts(info_b4_filt.paradigm.Pulse_2), dt);

figure('Position',[100 100 550 780])
m=max(max(abs(badata_b4_filt(keep,:))));
subplot(2,1,1); plot(badata_b4_filt(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (seconds)'), ylabel('log(\phi/\phi_0)') % plot signals 
subplot(2,1,2); imagesc(badata_b4_filt(keep,:),[-1,1].*m), colorbar('Location','northoutside'), xlabel('Time (seconds)'), ylabel('Measurement #') % show signals as image

keepd1=info_b4_filt.MEAS.GI & info_b4_filt.pairs.r2d<20 & info_b4_filt.pairs.WL==2;
keepd2=info_b4_filt.MEAS.GI & info_b4_filt.pairs.r2d>=20 & info_b4_filt.pairs.r2d<30 & info_b4_filt.pairs.WL==2;
keepd3=info_b4_filt.MEAS.GI & info_b4_filt.pairs.r2d>=30 & info_b4_filt.pairs.r2d<35 & info_b4_filt.pairs.WL==2;

figure('Position',[100 100 300 700])
subplot(3,1,1), plot(badata_b4_filt(keepd1,:)'), ylabel('R_{sd} < 20 mm'); ylim([-0.1,0.1]); xlim([1,30])
subplot(3,1,2), plot(badata_b4_filt(keepd2,:)'), ylabel('R_{sd} \in [20 30] mm');ylim([-0.1,0.1]); xlim([1,30])
subplot(3,1,3), plot(badata_b4_filt(keepd3,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (seconds)');ylim([-0.1,0.1]); xlim([1,30])


%% Logmean Light Levels
[lmdata, info.MEAS.Phi_0] = logmean(data);


%% Detect Noisy Channels
info = FindGoodMeas(lmdata, info, 0.075, flags);

% Example visualization
keep = info.pairs.WL==2 & info.pairs.r2d < 35 & info.MEAS.GI; % measurements to include

figure('Position',[100 100 550 780])
subplot(3,1,1); plot(lmdata(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (samples)'), ylabel('log(\phi/\phi_0)') % plot signals 
m=max(max(abs(lmdata(keep,:))));
subplot(3,1,2); imagesc(lmdata(keep,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (samples)'), ylabel('Measurement #') % show signals as image
[ftmag0,ftdomain] = fft_tts(lmdata(keep,:),info.system.framerate); % Generate average spectrum
ftmag = rms(ftmag0,1);
subplot(3,1,3); semilogx(ftdomain,ftmag), xlabel('Frequency (Hz)'), ylabel('|X(f)|') % plot vs. log frequency
xlim([1e-3 10])


%% Show nn1, nn2, nn3 (plots)
keepd1=info.MEAS.GI & info.pairs.r2d<20 & info.pairs.WL==2;
keepd2=info.MEAS.GI & info.pairs.r2d>=20 & info.pairs.r2d<30 & info.pairs.WL==2;
keepd3=info.MEAS.GI & info.pairs.r2d>=30 & info.pairs.r2d<35 & info.pairs.WL==2;

figure('Position',[100 100 1000 700])
subplot(3,1,1), plot(lmdata(keepd1,:)'), ylabel('R_{sd} < 20 mm')
subplot(3,1,2), plot(lmdata(keepd2,:)'), ylabel('R_{sd} \in [20 30] mm')
subplot(3,1,3), plot(lmdata(keepd3,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (samples)')


%% Detrend and High-pass Filter the Data
ddata = detrend_tts(lmdata);

% High Pass Filter
hpdata = highpass(ddata, 0.02, info.system.framerate);
% hpdata = highpass(ddata, 0.05, info.system.framerate); % problematic cutoff frequency example

figure('Position',[100 100 550 780])
subplot(3,1,1); plot(hpdata(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (samples)'), ylabel('log(\phi/\phi_0)') % plot signals 
m=max(max(abs(hpdata(keep,:))));
subplot(3,1,2); imagesc(hpdata(keep,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (samples)'), ylabel('Measurement #') % show signals as image
[ftmag0,ftdomain] = fft_tts(hpdata(keep,:),info.system.framerate); % Generate average spectrum
ftmag = rms(ftmag0,1);
subplot(3,1,3); semilogx(ftdomain,ftmag), xlabel('Frequency (Hz)'), ylabel('|X(f)|') % plot vs. log frequency
xlim([1e-3 10])


%%
figure('Position',[100 100 1000 700])
subplot(3,1,1), plot(hpdata(keepd1,:)'), ylabel('R_{sd} < 20 mm')
subplot(3,1,2), plot(hpdata(keepd2,:)'), ylabel('R_{sd} \in [20 30] mm')
subplot(3,1,3), plot(hpdata(keepd3,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (samples)')


%% Low Pass Filter 1
lp1data = lowpass(hpdata, 1, info.system.framerate);

figure('Position',[100 100 550 780])
subplot(3,1,1); plot(lp1data(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (samples)'), ylabel('log(\phi/\phi_0)') % plot signals 
m=max(max(abs(lp1data(keep,:))));
subplot(3,1,2); imagesc(lp1data(keep,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (samples)'), ylabel('Measurement #') % show signals as image
[ftmag0,ftdomain] = fft_tts(lp1data(keep,:),info.system.framerate); % Generate average spectrum
ftmag = rms(ftmag0,1);
subplot(3,1,3); semilogx(ftdomain,ftmag), xlabel('Frequency (Hz)'), ylabel('|X(f)|') % plot vs. log frequency
xlim([1e-3 10])


%%
figure('Position',[100 100 1000 700])
subplot(3,1,1), plot(lp1data(keepd1,:)'), ylabel('R_{sd} < 20 mm')
subplot(3,1,2), plot(lp1data(keepd2,:)'), ylabel('R_{sd} \in [20 30] mm')
subplot(3,1,3), plot(lp1data(keepd3,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (samples)')


%% Try Block Averaging Before Superficial Signal Regression (Left finger tapping)
badata_L_before_SSR = BlockAverage(lp1data, info.paradigm.init_synchpts(info.paradigm.Pulse_2), dt*round(info.system.framerate));
badata_L_before_SSR=bsxfun(@minus,badata_L_before_SSR,mean(badata_L_before_SSR(:,1:4),2));

figure('Position',[100 100 550 780])
m=max(max(abs(badata_L_before_SSR(keep,:))));
subplot(2,1,1); plot(badata_L_before_SSR(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (seconds)'), ylabel('log(\phi/\phi_0)') % plot signals 
subplot(2,1,2); imagesc(badata_L_before_SSR(keep,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (seconds)'), ylabel('Measurement #') % show signals as image

%
figure('Position',[100 100 300 700])
subplot(3,1,1), plot(badata_L_before_SSR(keepd1,:)'), ylabel('R_{sd} < 20 mm'); ylim([-0.13,0.13]); 
subplot(3,1,2), plot(badata_L_before_SSR(keepd2,:)'), ylabel('R_{sd} \in [20 30] mm');ylim([-0.13,0.13]); 
subplot(3,1,3), plot(badata_L_before_SSR(keepd3,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (seconds)');ylim([-0.13,0.13]); 


%% Try Block Averaging Before Superficial Signal Regression (Right finger tapping)
badata_R_before_SSR = BlockAverage(lp1data, info.paradigm.init_synchpts(info.paradigm.Pulse_2), dt*round(info.system.framerate));
badata_R_before_SSR=bsxfun(@minus,badata_R_before_SSR,mean(badata_R_before_SSR(:,1:4),2));

figure('Position',[100 100 550 780])
m=max(max(abs(badata_R_before_SSR(keep,:))));
subplot(2,1,1); plot(badata_R_before_SSR(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (seconds)'), ylabel('log(\phi/\phi_0)') % plot signals 
subplot(2,1,2); imagesc(badata_R_before_SSR(keep,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (seconds)'), ylabel('Measurement #') % show signals as image

%
figure('Position',[100 100 300 700])
subplot(3,1,1), plot(badata_R_before_SSR(keepd1,:)'), ylabel('R_{sd} < 20 mm'); ylim([-0.13,0.13]); 
subplot(3,1,2), plot(badata_R_before_SSR(keepd2,:)'), ylabel('R_{sd} \in [20 30] mm');ylim([-0.13,0.13]); 
subplot(3,1,3), plot(badata_R_before_SSR(keepd3,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (seconds)');ylim([-0.13,0.13]); 


%% Superficial Signal Regression
hem = gethem(lp1data, info, 'r2d', [1,15]);                                % Superficial Signal Regression
[SSRdata, ~] = regcorr(lp1data, info, hem);
% SSRdata = lp1data; % example to ignore SSR

figure('Position',[100 100 550 780])
subplot(3,1,1); plot(SSRdata(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (samples)'), ylabel('log(\phi/\phi_0)') % plot signals 
m=max(max(abs(SSRdata(keep,:))));
subplot(3,1,2); imagesc(SSRdata(keep,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (samples)'), ylabel('Measurement #') % show signals as image
[ftmag0,ftdomain] = fft_tts(SSRdata(keep,:),info.system.framerate); % Generate average spectrum
ftmag = rms(ftmag0,1);
subplot(3,1,3); semilogx(ftdomain,ftmag), xlabel('Frequency (Hz)'), ylabel('|X(f)|') % plot vs. log frequency
xlim([1e-3 10])


%% Try Block Averaging Right After Superficial Signal Regression (Left finger tapping)
badata_L_after_SSR = BlockAverage(SSRdata, info.paradigm.init_synchpts(info.paradigm.Pulse_2), dt*round(info.system.framerate));
badata_L_after_SSR=bsxfun(@minus,badata_L_after_SSR,mean(badata_L_after_SSR(:,1:4),2));

figure('Position',[100 100 550 780])
m=max(max(abs(badata_L_after_SSR(keep,:))));
subplot(2,1,1); plot(badata_L_after_SSR(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (seconds)'), ylabel('log(\phi/\phi_0)') % plot signals 
subplot(2,1,2); imagesc(badata_L_after_SSR(keep,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (seconds)'), ylabel('Measurement #') % show signals as image
%
figure('Position',[100 100 300 700])
subplot(3,1,1), plot(badata_L_after_SSR(keepd1,:)'), ylabel('R_{sd} < 20 mm'); ylim([-0.13,0.13]); 
subplot(3,1,2), plot(badata_L_after_SSR(keepd2,:)'), ylabel('R_{sd} \in [20 30] mm');ylim([-0.13,0.13]); 
subplot(3,1,3), plot(badata_L_after_SSR(keepd3,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (seconds)');ylim([-0.13,0.13]); 


%% Plot Superficial Signal
figure('Position',[100 100 800 430])
subplot(2,1,1)
plot(hem(2,:))
title('Estimated common superficial signal')
xlabel('Time (samples)')
subplot(2,1,2)
[ftmag,ftdomain] = fft_tts(hem(2,:),info.system.framerate); 
semilogx(ftdomain,ftmag), xlabel('Frequency (Hz)'), ylabel('|X(f)|') % plot vs. log frequency
xlim([1e-3 10])


%%
figure('Position',[100 100 1000 700])
subplot(3,1,1), plot(SSRdata(keepd1,:)'), ylabel('R_{sd} < 20 mm')
subplot(3,1,2), plot(SSRdata(keepd2,:)'), ylabel('R_{sd} \in [20 30] mm')
subplot(3,1,3), plot(SSRdata(keepd3,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (samples)')


%% Calculate GVTD
[info.GVTD, info.DQ_metrics.med_GVTD] = CalcGVTD(SSRdata(info.MEAS.GI & info.pairs.r2d<20,:));         % Calculate GVTD


%% Low Pass Filter 2
lp2data = lowpass(SSRdata, 0.05, info.system.framerate);

figure('Position',[100 100 550 780])
m=max(max(abs(lp2data(keep,:))));
subplot(3,1,1); plot(lp2data(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (samples)'), ylabel('log(\phi/\phi_0)') % plot signals 
subplot(3,1,2); imagesc(lp2data(keep,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (samples)'), ylabel('Measurement #') % show signals as image
[ftmag0,ftdomain] = fft_tts(lp2data(keep,:),info.system.framerate); % Generate average spectrum
ftmag = rms(ftmag0,1);
subplot(3,1,3); semilogx(ftdomain,ftmag), xlabel('Frequency (Hz)'), ylabel('|X(f)|') % plot vs. log frequency
xlim([1e-3 10])


%% Global variance of the temporal derivative (GVTD)
info.GVTD_filt_rs=resample_tts(info.GVTD',info,1,1e-5,...
    info.system.init_framerate)';


%% 1 Hz Resampling
[rdata, info] = resample_tts(lp2data, info, 1, 1e-5);
nlrGrayPlots_180818(rdata,info) % Plot GVTD 

figure('Position',[100 100 550 780])
m=max(max(abs(rdata(keep,:))));
subplot(3,1,1); plot(rdata(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (samples)'), ylabel('log(\phi/\phi_0)') % plot signals 
subplot(3,1,2); imagesc(rdata(keep,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (samples)'), ylabel('Measurement #') % show signals as image
[ftmag0,ftdomain] = fft_tts(rdata(keep,:),info.system.framerate); % Generate average spectrum
ftmag = rms(ftmag0,1);
subplot(3,1,3); semilogx(ftdomain,ftmag), xlabel('Frequency (Hz)'), ylabel('|X(f)|') % plot vs. log frequency
xlim([1e-3 1])


%%
figure('Position',[100 100 450 700])
subplot(3,1,1), plot(rdata(keepd1,:)'), ylabel('R_{sd} < 20 mm'),xlim([0 219])
subplot(3,1,2), plot(rdata(keepd2,:)'), ylabel('R_{sd} \in [20 30] mm'),xlim([0 219])
subplot(3,1,3), plot(rdata(keepd3,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (samples)'),xlim([0 219])


%% Block Averaging (Left finger tapping)
badata_L = BlockAverage(rdata, info.paradigm.synchpts(info.paradigm.Pulse_2), dt);
badata_L=bsxfun(@minus,badata_L,mean(badata_L(:,1:4),2));

figure('Position',[100 100 550 780])
m=max(max(abs(badata_L(keep,:))));
subplot(2,1,1); plot(badata_L(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (seconds)'), ylabel('log(\phi/\phi_0)') % plot signals 
subplot(2,1,2); imagesc(badata_L(keep,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (seconds)'), ylabel('Measurement #') % show signals as image


%%
figure('Position',[100 100 300 700])
subplot(3,1,1), plot(badata_L(keepd1,:)'), ylabel('R_{sd} < 20 mm'); ylim([-0.13,0.13]); xlim([1,34])
subplot(3,1,2), plot(badata_L(keepd2,:)'), ylabel('R_{sd} \in [20 30] mm');ylim([-0.13,0.13]); xlim([1,34])
subplot(3,1,3), plot(badata_L(keepd3,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (seconds)');ylim([-0.13,0.13]); xlim([1,34])

tp = 15;
params = struct; params.line_thickness = 1.5; params.Th.P=0; params.Th.N=0;
params.markBadOptodes=0; params.Scale=0.05;
infod1 = info;infod2 = info; infod3 = info;
infod1.MEAS.GI = infod1.MEAS.GI & keepd1; 
infod2.MEAS.GI = infod2.MEAS.GI & keepd2;
infod3.MEAS.GI = infod3.MEAS.GI & keepd3;  
PlotCapMeasSpaceData_Pairs(badata_L(:,tp), infod1,params)
PlotCapMeasSpaceData_Pairs(badata_L(:,tp), infod2,params)
PlotCapMeasSpaceData_Pairs(badata_L(:,tp), infod3,params)


%% Block Averaging (Right finger tapping)
badata_R = BlockAverage(rdata, info.paradigm.synchpts(info.paradigm.Pulse_3), dt);
badata_R=bsxfun(@minus,badata_R,mean(badata_R(:,1:4),2));

figure('Position',[100 100 550 780])
m=max(max(abs(badata_R(keep,:))));
subplot(2,1,1); plot(badata_R(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (seconds)'), ylabel('log(\phi/\phi_0)') % plot signals 
subplot(2,1,2); imagesc(badata_R(keep,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (seconds)'), ylabel('Measurement #') % show signals as image


%%
figure('Position',[100 100 300 700])
subplot(3,1,1), plot(badata_R(keepd1,:)'), ylabel('R_{sd} < 20 mm'); ylim([-0.1,0.1]); xlim([1,34])
subplot(3,1,2), plot(badata_R(keepd2,:)'), ylabel('R_{sd} \in [20 30] mm');ylim([-0.1,0.1]); xlim([1,34])
subplot(3,1,3), plot(badata_R(keepd3,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (seconds)');ylim([-0.1,0.1]); xlim([1,34])

tp = 15;
params = struct; params.line_thickness = 1.5; params.Th.P=0; params.Th.N=0;
params.markBadOptodes=0; params.Scale=0.05;
infod1 = info;infod2 = info; infod3 = info;
infod1.MEAS.GI = infod1.MEAS.GI & keepd1; 
infod2.MEAS.GI = infod2.MEAS.GI & keepd2;
infod3.MEAS.GI = infod3.MEAS.GI & keepd3;  
PlotCapMeasSpaceData_Pairs(badata_R(:,tp), infod1,params)
PlotCapMeasSpaceData_Pairs(badata_R(:,tp), infod2,params)
PlotCapMeasSpaceData_Pairs(badata_R(:,tp), infod3,params)


%% Try Block Averaging with GVTD Motion Thresholding
tkeep = info.GVTD < 1e-3; % Default threshold is 1e-3

% Block Averaging (Left finger tapping)
badata_L = BlockAverage(rdata, info.paradigm.synchpts(info.paradigm.Pulse_2), dt,tkeep);
badata_L=bsxfun(@minus,badata_L,mean(badata_L(:,1:4),2));

figure('Position',[100 100 550 780])
m=max(max(abs(badata_L(keep,:))));
subplot(2,1,1); plot(badata_L(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (seconds)'), ylabel('log(\phi/\phi_0)') % plot signals 
subplot(2,1,2); imagesc(badata_L(keep,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (seconds)'), ylabel('Measurement #') % show signals as image

%
figure('Position',[100 100 300 700])
subplot(3,1,1), plot(badata_L(keepd1,:)'), ylabel('R_{sd} < 20 mm'); ylim([-0.13,0.13]); xlim([1,34])
subplot(3,1,2), plot(badata_L(keepd2,:)'), ylabel('R_{sd} \in [20 30] mm');ylim([-0.13,0.13]); xlim([1,34])
subplot(3,1,3), plot(badata_L(keepd3,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (seconds)');ylim([-0.13,0.13]); xlim([1,34])

% Block Averaging (Right finger tapping)
badata_R = BlockAverage(rdata, info.paradigm.synchpts(info.paradigm.Pulse_3), dt,tkeep);
badata_R=bsxfun(@minus,badata_R,mean(badata_R(:,1:4),2));

figure('Position',[100 100 550 780])
m=max(max(abs(badata_R(keep,:))));
subplot(2,1,1); plot(badata_R(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (seconds)'), ylabel('log(\phi/\phi_0)') % plot signals 
subplot(2,1,2); imagesc(badata_R(keep,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (seconds)'), ylabel('Measurement #') % show signals as image

%
figure('Position',[100 100 300 700])
subplot(3,1,1), plot(badata_R(keepd1,:)'), ylabel('R_{sd} < 20 mm'); ylim([-0.1,0.1]); xlim([1,34])
subplot(3,1,2), plot(badata_R(keepd2,:)'), ylabel('R_{sd} \in [20 30] mm');ylim([-0.1,0.1]); xlim([1,34])
subplot(3,1,3), plot(badata_R(keepd3,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (seconds)');ylim([-0.1,0.1]); xlim([1,34])


%% Use the Modified Beer-Lambert Law to Block-Average HbO and HbR
ppf = [1, 1]; % select partial path length factors for each wavelength.
load('E.mat');
[HbO,HbR] = MBLL(rdata, info, ppf, E);

% Block Averaging HbO and HbR - Left finger tapping
badata_L_HbX = BlockAverage([HbO;HbR], info.paradigm.synchpts(info.paradigm.Pulse_2), dt);
badata_L_HbX=bsxfun(@minus,badata_L_HbX,mean(badata_L_HbX(:,1:4),2));

keep =  info.pairs.r2d < 35 & info.MEAS.GI; % measurements to include

% Plot HbO - Left Finger Tapping
nMeas = size(HbO,1);
keep_HbO = keep(1:nMeas);
keep_HbR = keep(nMeas+1:end);

badata_L_HbO = badata_L_HbX(1:nMeas,:);
badata_L_HbR = badata_L_HbX (nMeas+1:end,:);
figure('Position',[100 100 550 780])
m=max(max(abs(badata_L_HbO(keep_HbO,:))));
subplot(2,1,1); plot(badata_L_HbO(keep_HbO,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (seconds)'), ylabel('\Delta HbO') % plot signals 
subplot(2,1,2); imagesc(badata_L_HbO(keep_HbO,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (seconds)'), ylabel('Measurement #') % show signals as image

%
keepd1=info.MEAS.GI & info.pairs.r2d<20;
keepd2=info.MEAS.GI & info.pairs.r2d>=20 & info.pairs.r2d<30;
keepd3=info.MEAS.GI & info.pairs.r2d>=30 & info.pairs.r2d<35;

keepd1_HbO = keepd1(1:nMeas);
keepd2_HbO = keepd2(1:nMeas);
keepd3_HbO = keepd3(1:nMeas);
figure('Position',[100 100 300 700])
subplot(3,1,1), plot(badata_L_HbO(keepd1_HbO,:)'), ylabel('R_{sd} < 20 mm'); ylim([-0.13,0.13]); xlim([1,34])
subplot(3,1,2), plot(badata_L_HbO(keepd2_HbO,:)'), ylabel('R_{sd} \in [20 30] mm');ylim([-0.13,0.13]); xlim([1,34])
subplot(3,1,3), plot(badata_L_HbO(keepd3_HbO,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (seconds)');ylim([-0.13,0.13]); xlim([1,34])

%to use on hb data need to concatenate hbo and hbr and run pcmss. maybe
%make half not using empty
HbX_L = [badata_L_HbO;badata_L_HbR];
HbO_L = HbX_L;
tp = 15;
params = struct; params.line_thickness = 1.5; params.Th.P=0; params.Th.N=0;
params.markBadOptodes=0; params.Scale=0.05; params.WLs = 1;
infod1 = info; infod2 = info; infod3 = info;
infod1.MEAS.GI = infod1.MEAS.GI & keepd1; 
infod2.MEAS.GI = infod2.MEAS.GI& keepd2;
infod3.MEAS.GI = infod3.MEAS.GI& keepd3;  
PlotCapMeasSpaceData_Pairs(HbO_L(:,tp), infod1, params)
PlotCapMeasSpaceData_Pairs(HbO_L(:,tp), infod2, params)
PlotCapMeasSpaceData_Pairs(HbO_L(:,tp), infod3, params)

% Plot HbR - Left Finger Tapping
figure('Position',[100 100 550 780])
m=max(max(abs(badata_L_HbR(keep_HbR,:))));
subplot(2,1,1); plot(badata_L_HbR(keep_HbR,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (seconds)'), ylabel('\Delta HbR') % plot signals 
subplot(2,1,2); imagesc(badata_L_HbR(keep_HbR,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (seconds)'), ylabel('Measurement #') % show signals as image

%
keepd1=info.MEAS.GI & info.pairs.r2d<20;
keepd2=info.MEAS.GI & info.pairs.r2d>=20 & info.pairs.r2d<30;
keepd3=info.MEAS.GI & info.pairs.r2d>=30 & info.pairs.r2d<35;
keepd1_HbO = keepd1(1:nMeas);
keepd2_HbO = keepd2(1:nMeas);
keepd3_HbO = keepd3(1:nMeas);
figure('Position',[100 100 300 700])
subplot(3,1,1), plot(badata_L_HbR(keepd1_HbO,:)'), ylabel('R_{sd} < 20 mm'); ylim([-0.05,0.05]); xlim([1,34])
subplot(3,1,2), plot(badata_L_HbR(keepd2_HbO,:)'), ylabel('R_{sd} \in [20 30] mm');ylim([-0.05,0.05]); xlim([1,34])
subplot(3,1,3), plot(badata_L_HbR(keepd3_HbO,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (seconds)');ylim([-0.05,0.05]); xlim([1,34])

HbX_L = [badata_L_HbO;badata_L_HbR];
HbO_L = HbX_L;
tp = 15;
params = struct; params.line_thickness = 1.5; params.Th.P=0; params.Th.N=0;
params.markBadOptodes=0; params.Scale=0.05; params.WLs = 2;
infod1 = info; infod2 = info; infod3 = info;
infod1.MEAS.GI = infod1.MEAS.GI & keepd1; 
infod2.MEAS.GI = infod2.MEAS.GI& keepd2;
infod3.MEAS.GI = infod3.MEAS.GI& keepd3;  
PlotCapMeasSpaceData_Pairs(HbO_L(:,tp), infod1, params)
PlotCapMeasSpaceData_Pairs(HbO_L(:,tp), infod2, params)
PlotCapMeasSpaceData_Pairs(HbO_L(:,tp), infod3, params)

% Block Averaging HbO and HbR - Right finger tapping
badata_R_HbX = BlockAverage([HbO;HbR], info.paradigm.synchpts(info.paradigm.Pulse_3), dt);
badata_R_HbX=bsxfun(@minus,badata_R_HbX,mean(badata_R_HbX(:,1:4),2));

keep =  info.pairs.r2d < 35 & info.MEAS.GI; % measurements to include

% Plot HbO - Right Finger Tapping
nMeas = size(HbO,1);
keep_HbO = keep(1:nMeas);
keep_HbR = keep(nMeas+1:end);
badata_R_HbO = badata_R_HbX(1:nMeas,:);

badata_R_HbR = badata_R_HbX(nMeas+1:end,:);
figure('Position',[100 100 550 780])
m=max(max(abs(badata_R_HbO(keep_HbO,:))));
subplot(2,1,1); plot(badata_R_HbO(keep_HbO,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (seconds)'), ylabel('\Delta HbO') % plot signals 
subplot(2,1,2); imagesc(badata_R_HbO(keep_HbO,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (seconds)'), ylabel('Measurement #') % show signals as image

%
keepd1=info.MEAS.GI & info.pairs.r2d<20;
keepd2=info.MEAS.GI & info.pairs.r2d>=20 & info.pairs.r2d<30;
keepd3=info.MEAS.GI & info.pairs.r2d>=30 & info.pairs.r2d<35;

keepd1_HbO = keepd1(1:nMeas);
keepd2_HbO = keepd2(1:nMeas);
keepd3_HbO = keepd3(1:nMeas);
figure('Position',[100 100 300 700])
subplot(3,1,1), plot(badata_R_HbO(keepd1_HbO,:)'), ylabel('R_{sd} < 20 mm'); ylim([-0.05,0.05]); xlim([1,34])
subplot(3,1,2), plot(badata_R_HbO(keepd2_HbO,:)'), ylabel('R_{sd} \in [20 30] mm');ylim([-0.05,0.05]); xlim([1,34])
subplot(3,1,3), plot(badata_R_HbO(keepd3_HbO,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (seconds)');ylim([-0.05,0.05]); xlim([1,34])

HbX_R = [badata_R_HbO;badata_R_HbR];
HbO_R = HbX_R;
tp = 15;
params = struct; params.line_thickness = 1.5; params.Th.P=0; params.Th.N=0;
params.markBadOptodes=0; params.Scale=0.05; params.WLs = 1;
infod1 = info; infod2 = info; infod3 = info;
infod1.MEAS.GI = infod1.MEAS.GI & keepd1; 
infod2.MEAS.GI = infod2.MEAS.GI & keepd2;
infod3.MEAS.GI = infod3.MEAS.GI & keepd3;  
PlotCapMeasSpaceData_Pairs(HbO_R(:,tp), infod1, params)
PlotCapMeasSpaceData_Pairs(HbO_R(:,tp), infod2, params)
PlotCapMeasSpaceData_Pairs(HbO_R(:,tp), infod3, params)

% Plot HbR - Right Finger Tapping
figure('Position',[100 100 550 780])
m=max(max(abs(badata_R_HbR(keep_HbR,:))));
subplot(2,1,1); plot(badata_R_HbR(keep_HbR,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (seconds)'), ylabel('\Delta HbR') % plot signals 
subplot(2,1,2); imagesc(badata_R_HbR(keep_HbR,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (seconds)'), ylabel('Measurement #') % show signals as image

%
keepd1=info.MEAS.GI & info.pairs.r2d<20;
keepd2=info.MEAS.GI & info.pairs.r2d>=20 & info.pairs.r2d<30;
keepd3=info.MEAS.GI & info.pairs.r2d>=30 & info.pairs.r2d<35;
keepd1_HbR = keepd1(nMeas+1:end);
keepd2_HbR = keepd2(nMeas+1:end);
keepd3_HbR = keepd3(nMeas+1:end);
figure('Position',[100 100 300 700])
subplot(3,1,1), plot(badata_R_HbR(keepd1_HbR,:)'), ylabel('R_{sd} < 20 mm'); ylim([-0.05,0.05]); xlim([1,34])
subplot(3,1,2), plot(badata_R_HbR(keepd2_HbR,:)'), ylabel('R_{sd} \in [20 30] mm');ylim([-0.05,0.05]); xlim([1,34])
subplot(3,1,3), plot(badata_R_HbR(keepd3_HbR,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (seconds)');ylim([-0.05,0.05]); xlim([1,34])

HbX_R = [badata_R_HbO;badata_R_HbR];
HbR_R = HbX_R;
tp = 15;
params = struct; params.line_thickness = 1.5; params.Th.P=0; params.Th.N=0;
params.markBadOptodes=0; params.Scale=0.05; params.WLs = 2;
infod1 = info; infod2 = info; infod3 = info;
infod1.MEAS.GI = infod1.MEAS.GI & keepd1; 
infod2.MEAS.GI = infod2.MEAS.GI & keepd2;
infod3.MEAS.GI = infod3.MEAS.GI & keepd3;  
PlotCapMeasSpaceData_Pairs(HbR_R(:,tp), infod1, params)
PlotCapMeasSpaceData_Pairs(HbR_R(:,tp), infod2, params)
PlotCapMeasSpaceData_Pairs(HbR_R(:,tp), infod3, params)


%% Preprocessing Pipeline Without Superficial Signal Regression
clear all;close all;clc;
load(['ND_LUMO_Data_Sample_MOT001.mat']); % data, info, flags
A_fn='A_LUMO_51tile_KKI_on_HD_Mesh0_r3dlt50mm.mat';   % Sensitivity Matrix
flags.RsdMax = 35;
flags.AbsDiff = 0.5;
flags.FloorTh = 10e-8;
flags.DtKeep = 1;
dt=34;  % Block length
tp=10;  % Example (block averaged) time point

% Crop data to match A Matrix' measurement list
% Load measurement list from A Matrix
temp = load(['A_LUMO_51tile_KKI_on_HD_Mesh0_r3dlt50mm.mat'],'info'); 

% Crop data
mA=cat(2,temp.info.pairs.Src,temp.info.pairs.Det,temp.info.pairs.WL);
mD=cat(2,info.pairs.Src,info.pairs.Det,info.pairs.WL);
keepInfo=ismember(mD,mA,'rows');
data=data(keepInfo,:);                                                    % Choose measurementps from data
info.pairs=temp.info.pairs;                                               % Grab relevant measurements from pairs with A
info.MEAS.Clipped = info.MEAS.Clipped(keepInfo);

[lmdata, info.MEAS.Phi_0] = logmean(data);                                % Logmean Light Levels
info = FindGoodMeas(lmdata, info, 0.075, flags);  
keep = info.pairs.WL==2 & info.pairs.r2d < 35 & info.MEAS.GI;             % measurements to include
lmdata = detrend_tts(lmdata);                                             % Detrend Data
lmdata = highpass(lmdata, .02, info.system.framerate);                    % High Pass Filter (0.02 Hz)
lmdata = lowpass(lmdata, 1, info.system.framerate);                       % Low Pass Filter 1 (1.0 Hz)
lmdata = lowpass(lmdata, 0.2, info.system.framerate);                     % Low Pass Filter 2 (0.2 Hz)
info.GVTD = CalcGVTD(lmdata(info.MEAS.GI & info.pairs.r2d<20,:));         % Calculate GVTD
[lmdata, info] = resample_tts(lmdata, info, 1, 1e-5);                     % 1 Hz Resampling (1 Hz)
info.GVTD_filt_rs = resample_tts(info.GVTD', info, 1, 1e-5,...
info.system.init_framerate)';                                             % Resample GVTD to 1 Hz

% Block Averaging (Left finger tapping)
badata_L = BlockAverage(lmdata, info.paradigm.synchpts(info.paradigm.Pulse_2), dt);
badata_L=bsxfun(@minus,badata_L,mean(badata_L(:,1:4),2));
figure('Position',[100 100 550 780])
m=max(max(abs(badata_L(keep,:))));
subplot(2,1,1); plot(badata_L(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (seconds)'), ylabel('log(\phi/\phi_0)') % plot signals 
subplot(2,1,2); imagesc(badata_L(keep,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (seconds)'), ylabel('Measurement #') % show signals as image

keepd1=info.MEAS.GI & info.pairs.r2d<20 & info.pairs.WL==2;
keepd2=info.MEAS.GI & info.pairs.r2d>=20 & info.pairs.r2d<30 & info.pairs.WL==2;
keepd3=info.MEAS.GI & info.pairs.r2d>=30 & info.pairs.r2d<35 & info.pairs.WL==2;

figure('Position',[100 100 300 700])
subplot(3,1,1), plot(badata_L(keepd1,:)'), ylabel('R_{sd} < 20 mm'); ylim([-0.06,0.06]); xlim([1,34])
subplot(3,1,2), plot(badata_L(keepd2,:)'), ylabel('R_{sd} \in [20 30] mm');ylim([-0.06,0.06]); xlim([1,34])
subplot(3,1,3), plot(badata_L(keepd3,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (seconds)');ylim([-0.06,0.06]); xlim([1,34])

% Block Averaging (Right finger tapping)
badata_R = BlockAverage(lmdata, info.paradigm.synchpts(info.paradigm.Pulse_3), dt);
badata_R=bsxfun(@minus,badata_R,mean(badata_R(:,1:4),2));

figure('Position',[100 100 550 780])
m=max(max(abs(badata_R(keep,:))));
subplot(2,1,1); plot(badata_R(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (seconds)'), ylabel('log(\phi/\phi_0)') % plot signals 
subplot(2,1,2); imagesc(badata_R(keep,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (seconds)'), ylabel('Measurement #') % show signals as image

figure('Position',[100 100 300 700])
subplot(3,1,1), plot(badata_R(keepd1,:)'), ylabel('R_{sd} < 20 mm'); ylim([-0.1,0.1]); xlim([1,34])
subplot(3,1,2), plot(badata_R(keepd2,:)'), ylabel('R_{sd} \in [20 30] mm');ylim([-0.1,0.1]); xlim([1,34])
subplot(3,1,3), plot(badata_R(keepd3,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (seconds)');ylim([-0.1,0.1]); xlim([1,34])


%% Preprocessing Pipeline With Higher Low-Pass Filter (0.5 Hz)
clear all;close all;clc;
load(['ND_LUMO_Data_Sample_MOT001.mat']); % data, info, flags
A_fn='A_LUMO_51tile_KKI_on_HD_Mesh0_r3dlt50mm.mat';   % Sensitivity Matrix
flags.RsdMax = 35;
flags.AbsDiff = 0.5;
flags.FloorTh = 10e-8;
flags.DtKeep = 1;
dt=34;  % Block length
tp=10;  % Example (block averaged) time point

% Crop data to match A Matrix' measurement list
% Load measurement list from A Matrix
temp = load(['A_LUMO_51tile_KKI_on_HD_Mesh0_r3dlt50mm.mat'],'info'); 

% Crop data
mA=cat(2,temp.info.pairs.Src,temp.info.pairs.Det,temp.info.pairs.WL);
mD=cat(2,info.pairs.Src,info.pairs.Det,info.pairs.WL);
keepInfo=ismember(mD,mA,'rows');
data=data(keepInfo,:);                                                    % Choose measurementps from data
info.pairs=temp.info.pairs;                                               % Grab relevant measurements from pairs with A
info.MEAS.Clipped = info.MEAS.Clipped(keepInfo);

[lmdata, info.MEAS.Phi_0] = logmean(data);                                % Logmean Light Levels
info = FindGoodMeas(lmdata, info, 0.075, flags);  
keep = info.pairs.WL==2 & info.pairs.r2d < 35 & info.MEAS.GI;             % measurements to include
lmdata = detrend_tts(lmdata);                                             % Detrend Data
lmdata = highpass(lmdata, .02, info.system.framerate);                    % High Pass Filter (0.02 Hz)
lmdata = lowpass(lmdata, 1, info.system.framerate);                       % Low Pass Filter 1 (1.0 Hz)
hem = gethem(lmdata, info, 'r2d', [1,15]);                                % Superficial Signal Regression
[lmdata, ~] = regcorr(lmdata, info, hem);
lmdata = lowpass(lmdata, 0.5, info.system.framerate);                     % Low Pass Filter 2 (0.2 Hz)
info.GVTD = CalcGVTD(lmdata(info.MEAS.GI & info.pairs.r2d<20,:));         % Calculate GVTD
[lmdata, info] = resample_tts(lmdata, info, 1, 1e-5);                     % 1 Hz Resampling (1 Hz)
info.GVTD_filt_rs = resample_tts(info.GVTD', info, 1, 1e-5,...
info.system.init_framerate)';                                         % Resample GVTD to 1 Hz

% Block Averaging (Left finger tapping)
badata_L = BlockAverage(lmdata, info.paradigm.synchpts(info.paradigm.Pulse_2), dt);
badata_L=bsxfun(@minus,badata_L,mean(badata_L(:,1:4),2));
figure('Position',[100 100 550 780])
m=max(max(abs(badata_L(keep,:))));
subplot(2,1,1); plot(badata_L(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (seconds)'), ylabel('log(\phi/\phi_0)') % plot signals 
subplot(2,1,2); imagesc(badata_L(keep,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (seconds)'), ylabel('Measurement #') % show signals as image

keepd1=info.MEAS.GI & info.pairs.r2d<20 & info.pairs.WL==2;
keepd2=info.MEAS.GI & info.pairs.r2d>=20 & info.pairs.r2d<30 & info.pairs.WL==2;
keepd3=info.MEAS.GI & info.pairs.r2d>=30 & info.pairs.r2d<35 & info.pairs.WL==2;

figure('Position',[100 100 300 700])
subplot(3,1,1), plot(badata_L(keepd1,:)'), ylabel('R_{sd} < 20 mm'); ylim([-0.06,0.06]); xlim([1,34])
subplot(3,1,2), plot(badata_L(keepd2,:)'), ylabel('R_{sd} \in [20 30] mm');ylim([-0.06,0.06]); xlim([1,34])
subplot(3,1,3), plot(badata_L(keepd3,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (seconds)');ylim([-0.06,0.06]); xlim([1,34])

% Block Averaging (Right finger tapping)
badata_R = BlockAverage(lmdata, info.paradigm.synchpts(info.paradigm.Pulse_3), dt);
badata_R=bsxfun(@minus,badata_R,mean(badata_R(:,1:4),2));

figure('Position',[100 100 550 780])
m=max(max(abs(badata_R(keep,:))));
subplot(2,1,1); plot(badata_R(keep,:)'), set(gca,'XLimSpec','tight'), xlabel('Time (seconds)'), ylabel('log(\phi/\phi_0)') % plot signals 
subplot(2,1,2); imagesc(badata_R(keep,:),[-0.2,0.2].*m), colorbar('Location','northoutside'), xlabel('Time (seconds)'), ylabel('Measurement #') % show signals as image

figure('Position',[100 100 300 700])
subplot(3,1,1), plot(badata_R(keepd1,:)'), ylabel('R_{sd} < 20 mm'); ylim([-0.1,0.1]); xlim([1,34])
subplot(3,1,2), plot(badata_R(keepd2,:)'), ylabel('R_{sd} \in [20 30] mm');ylim([-0.1,0.1]); xlim([1,34])
subplot(3,1,3), plot(badata_R(keepd3,:)'), ylabel('R_{sd} \in [30 35] mm'), xlabel('Time (seconds)');ylim([-0.1,0.1]); xlim([1,34])



