function [ROI, Cond] = roi_preprocYmad(Action,ROI,Cond,Vf,RT);
% Temporal preprocessing of raw fMRI data from Ymad.
%  
% FORMAT ROI = roi_preproc(Action,ROI,Cond,Vf,RT);
%
% Action - 'Trim'   : replace outliers > 2 st dev with trim mean
%        - 'GMScale': session specific grand mean scaling
%        - 'Detrend': linear detrending within session
%        - 'LP'     : low-pass Butterworth filter
%        - 'HP'     : high-pass Butterworth filter
%        - 'All'    : all of the above
%        The above can also be combined, e.g. {'Trim','LP','HP'}
%
% Yraw   - a matrix of [#V-files x #subj] with ROI timeseries  
%          each column is the timeseries for a subject
%
%       - a subject structure with the following subfields
% Cond - a condition substructure with design information, 
%          of size(Cond) = [1 x #subj]
%          see kc_get_design for subfields
% Vf   - a V-files substructure with subfields
%     .iS  - V-files indexed by session
%     .G   - global mean activity for the corresponding V-file
% RT   - repetition time (in seconds)
%
% Updates the ROI structure with the field:
%     ROI.Ypr      - preprocessed ROI.Ymad_avg
% and the Cond structure with the field
%     Cond.Ypr_avg - pst averaged ROI.Yprfor this condition
% -----------------------------------------------------------------------
% @(#) roi_preprocYmad.m 1.1  K.Christoff 2000-08-25, modified 2001-01-20
% modified 2001-08-08: "Cond.Ypr_pst_avg" replaced with "Cond.Ypr_avg"
  
if nargin == 0, error('Insufficient number of input args ...'), end
if isempty(find(strmatch(Action,...
	   {'Trim','GMScale','Detrend','LP','HP','All','None'}))),
  error('Action string not recognized ...');
end
  
Y = ROI.Ymad_avg;
nsess = max(Vf.iS);
	 
 % ---------------------------------------------
 % Trim outliers by replacing abs(zscore)>2 with 
 % the sessoin specific trimmed mean (10%)
 % --------------------------------------------
 if ~isempty(strmatch('Trim',Action)) | ~isempty(strmatch('All',Action))
 % ---------------------------------------------
	 trimcut = 5; % for the trimmean cut off 10% (5% on each side)
	 prev_idx = 0;
	 for s = 1:nsess
		sess_idx = find(Vf.iS==s);
		out_idx = find(abs(zscore(Y(sess_idx))) > 3);
		trM	  = trimmean(Y(sess_idx),trimcut);
		Y(prev_idx+out_idx) = trM;
		prev_idx = sess_idx(end);
	 end
 end
	 
 % -----------------------------------------------
 % Convert to per cent signal change
 % -----------------------------------------------
 Y = (Y-mean(Y))*100/mean(Y);
 
 
 % --------------------------------------------------------------
 % Apply session specific grand-mean scaling, dividing by the 
 % trimmed (%10) grand mean for the entire brain for this session
 % --------------------------------------------------------------
 if ~isempty(strmatch('GMScale',Action)) | ~isempty(strmatch('All',Action))
 % --------------------------------------------------------------
	 trimcut = 10; % for the trimmean cut off 10% (5% on each side)
	 for s = 1:nsess
		sess_idx = find(Vf.iS==s);
		trGM = trimmean(Vf.G(sess_idx),trimcut);
		Y(sess_idx) = Y(sess_idx)*100/trGM;
	 end
 end

 
 % --------------------------------------------------------------
 % Detrend by removing linear trend - separately for each session
 % --------------------------------------------------------------
 if ~isempty(strmatch('Detrend',Action)) | ~isempty(strmatch('All',Action))
 % --------------------------------------------------------------
	 for s = 1:nsess
		sess_idx = find(Vf.iS==s);
		Y(sess_idx) = detrend(Y(sess_idx));
	 end
 end
	 

 % --------------------------------------------------------------
 % Low-pass filter the data using Butterworth filter
 % --------------------------------------------------------------
 if ~isempty(strmatch('LP',Action)) | ~isempty(strmatch('All',Action))
 % --------------------------------------------------------------
	 Nyq = 1/(RT*2);	% Nyquist frequency in Hz
	 cutoff = 1/8; %was24 	% HRF Hz cycle (~20s cycle or ~10s half-cycle)
	 order = 1;			% second order filter

	 [b,a] = butter(order,cutoff/Nyq);
	 Y = filter(b,a,Y);
 end
	 

 % --------------------------------------------------------------
 % High-pass filter the data using Butterworth filter
 % --------------------------------------------------------------
 if ~isempty(strmatch('HP',Action)) | ~isempty(strmatch('All',Action))
 % --------------------------------------------------------------
	 Nyq = 1/(RT*2);	% Nyquist frequency in Hz
	 
	 % use (apprx.) spm's high-pass filter formula: 
	 % cutoff = 2*(min(min_elapse)+RT);
%	 min_elapse = min(min(diff([Cond(:).onsec])));
%	 PoI = 2*min_elapse+RT; % period of interest (sec)
    PoI = 100;
    cutoff = 1/PoI;	
	 order = 1;			% second order filter
	 
	 % filter session by session
    [b,a] = butter(order,cutoff/Nyq,'high');
 	 for s = 1:nsess
		sess_idx = find(Vf.iS==s);
		Y(sess_idx) = filter(b,a,Y(sess_idx));
	 end

 end
	 

  
 % update output structues
 ROI.Ypr = Y;
 for c = 1:length(Cond)
	for i=1:Cond(c).num_scans
	     nVf =  size(Vf.iS,1);
	  	  tmp_idx = Cond(c).onidx + (i-1);
		  tmp_idx = tmp_idx(find(tmp_idx<=nVf)); % if trials exceed exp end 
        Cond(c).Ypr_avg(i)=mean(ROI.Ypr(tmp_idx));
	end
 end 
 
 
if nargout == 1; varargout{1} = ROI; end
if nargout == 2; varargout{2} = Cond; end
	
 
 






