


%%%%%

% load in Ms , Ps ...

Ms = {'a01', 'a03', 'a05', 'a06', 'a07', 'a09', 'a10','a13', 'a14', 'a15', 'a17', 'a18', 'a20'};

Ps = {'a01', 'a03', 'a05', 'a06', 'a07', 'a09', 'a10','a13', 'a14', 'a15', 'a17', 'a18', 'a20'};


% for each subject: read in HR vector (3000 x #voxels) -> store in M{subject_no)  & P{subject_no) ...


for s = 1:length(Ms)
   
	%LPM
	cd(Ms{2}); cd mr_results; pwd
	clear hrf; load('Lpmot6.mat'); M{s} = hrf;
	cd ('../..')

	%RP
	cd(Ms{2}); cd mr_results; pwd
	clear hrf; load('Rpar.mat'); P{s} = hrf;
   cd ('../..')
   
end




%%%%%%

% create figure window
figure
%%%%

no_subjects = size(Ms,2);


behavior_dir = 'Behavioral_data';
results_dir = 'results_dir'; % !!!! remember to create 'results_dir' folder !!!


rot_angle = 0; % condition of interest

% lower and upper bounds for desired RT range
RT_percentile_low = 25; 
RT_percentile_high = 75;

%% 
num_scans_per_trial = 20;
num_trials = 150;
num_scans_for_exp = 3000;


%%%%%% selectively average HRs for each Sj, and the designated condition and RT range
for s = 1:length(Ms) %%% loop over subjects .... %%%
   
   hrf_m = hrf_filter(hrf_m);    hrf_p = hrf_filter(hrf_p); %pre-processing
   
   
   %%% load in data matrix for current subject
   cwd = pwd;
   
   cd(behavior_dir)
   data_filename = [ 'Data' num2str(s) ];
   load(data_filename)
   
   Data = eval(data_filename);  % for each subject, use the same variable names
   
   %%% determine indices for trials for the desired condition
   ind_cond = find(Data(:,2) == 0);
   
   %%% extract trials within the RT range for this condition
   Data_RT = Data(ind_cond, 3) ; %the RT vector for this condition
   
   RT_LB = prctile(Data_RT , RT_percentile_low );
   RT_UB = prctile(Data_RT , RT_percentile_high );
   
   %find trials within RT range
   ind_RT = find(   ( Data(:,2) == rot_angle )  &  ( Data(:,3) > RT_LB ) & ( Data(:,3) <= RT_UB )   ) ;
   
   
   
   %%%%%%%%%%%%%%%%%% create 20x150 volume index to extract HR for trials of interest --- 20 = no_scans/trial; 150 = total # of trials in exp
   trial_inds_temp = (1:num_scans_for_exp)'; %%% 3000 = # of scans in whole exp
   trial_inds = reshape( trials_inds_temp , num_scans_per_trial , num_trials );
   volume_inds = trial_inds(:,ind_RT); %%% vol indices for desired trials within the image series for whole exp
   
   
   %%% extract HRs for each of the selected trials, for ROIs M & P -> assign to hr_xxx (no_scans_trial x no_trials)
   hr_trials_m = []; hr_trials_p = [];

	for j = 1:size( volume_inds , 2)
   
      hr_trials_m = [hr_trials_m , hrf_m(volume_inds(:,j))];
      hr_trials_p = [hr_trials_p , hrf_p(volume_inds(:,j))];

	end
   
   % selectively average over trials for both ROI
   hr_m = mean(hr_trials_m');
   hr_p = mean(hr_trials_p');
   
   %%% normalize HRs to an onset of 0 and a maximum amplitude of 1
   hr_m = hr_m - hr_m(1); hr_m = hr_m/max(hr_m);
   hr_p = hr_p - hr_p(1); hr_p = hr_p/max(hr_p);
   
   %%% get canonical HR
   hr_canon = spm_hrf(1);
   hr_canon = hr_canon(1:20)/max(hr_canon(1:20));
   
   %%%% plot this .....
   no_columns = 5; %%% can be changed
   no_rows = ceil( no_subjects / no_columns );
   
   subplot(no_rows , no_rows , s) ;
   
   hold on
   
   axis( 0 [num_scans_per_trial -1.5 1.5] )
   
   grid on
   box on
   set(gca,'Fontsize',8)   %%% font size is adjustable
   
   plot(hr_m, 'LineWidth',1);
   plot(hr_p,'g', 'LineWidth',1);
   plot(hr_c,'r--', 'LineWidth',1);
   
   text(13,0.8, Ms{s}(1:3));
   
   % store HRs for all Sjs in a cell array
   meanM{s} = hr_m;
   meanP{s} = hr_p;
   
   
end

title('Whole Nine Yards')

cd(cwd)

cd(results_dir)
save meanM meanP hrfc Ms Ps 'the_whole_kitt_and_kaboodle'







   
   
   
   
   
   
   


   
   
   
   
   
   
   
   
   
   
   
   