diff --git a/matlab_code/filtered_movement_regressors.m b/matlab_code/filtered_movement_regressors.m index ff58083..d11cef1 100755 --- a/matlab_code/filtered_movement_regressors.m +++ b/matlab_code/filtered_movement_regressors.m @@ -7,8 +7,6 @@ function filtered_movement_regressors(path_mov_reg, TR, filt_option, order, LP_f % filt_option filter option, options 1:6, see below % order order of the filter to be applied % LP_freq_min Low pass frequency (in minutes) to be filtered -%% defaults -head_ratio_cm = 5; %% when compiling comment out paths and uncomment this section TR = str2num(TR); @@ -27,43 +25,21 @@ function filtered_movement_regressors(path_mov_reg, TR, filt_option, order, LP_f filt_method{5} = 'FiltFilt_all'; filt_method{6} = 'FiltFilt_only_Y'; -switch filt_type +if filt_option>4 + order = order/2; % because filtfilt does forward and backward filtering +end +fs = 1 / TR; +fNy = fs / 2; +switch filt_type case 'lp' - hr_min = LP_freq_min; % recasted to reuse code hr = hr_min/60; - - fs = 1/TR; - fNy = fs/2; - - fa = abs(hr - floor((hr + fNy) / fs) * fs); - - % cutting frequency normalized between 0 and nyquist - Wn = min(fa)/fNy; - if ~isempty(order) - b_filt = fir1(order, Wn, 'low'); - a_filt = 1; - end - num_f_apply = 0; - + [b_filt,a_filt] = butter(order,rr/fNy,'stop');[b_filt,a_filt] = butter(order/2,hr/fNy,'low'); case 'notch' - fc_RR_bw = [fc_RR_min, fc_RR_max]; rr = fc_RR_bw / 60; - - fs = 1 / TR; - fNy = fs / 2; - - fa = abs(rr - floor((rr+fNy) / fs) * fs); - - W_notch = fa / fNy; - Wn = mean(W_notch); - Wd = diff(W_notch); - bw = abs(Wd); % take the absolute value in order to ensure that the difference between the min and max is not negative - [b_filt, a_filt] = iirnotch(Wn, bw); - num_f_apply = floor(order / 2); % if order<4 apply filter 1x, if order=4 2x, if order=6 3x - + [b_filt,a_filt] = butter(order,rr/fNy,'stop'); end @@ -86,8 +62,7 @@ function filtered_movement_regressors(path_mov_reg, TR, filt_option, order, LP_f MR_ld=make_friston_regressors(MR);%% Using this function to only get the linear displacements MR_ld=MR_ld(:,1:6); - switch filt_option - + switch filt_option case 1 %'None'; MR_filt = MR_ld; @@ -97,16 +72,10 @@ function filtered_movement_regressors(path_mov_reg, TR, filt_option, order, LP_f case 3 %'Filt_all'; MR_filt = filter(b_filt,a_filt,MR_ld); - for i=1:num_f_apply-1 - MR_filt = filter(b_filt,a_filt,MR_filt); - end case 4 %'Filt_only_Y'; MR_filt = MR_ld; MR_filt(:,2) = filter(b_filt,a_filt,MR_ld(:,2)); - for i=1:num_f_apply-1 - MR_filt(:,2) = filter(b_filt,a_filt,MR_filt(:,2)); - end case 5 %'FiltFilt_all'; MR_filt = filtfilt(b_filt,a_filt,MR_ld); @@ -117,15 +86,9 @@ function filtered_movement_regressors(path_mov_reg, TR, filt_option, order, LP_f case 6 %'FiltFilt_only_Y'; MR_filt = MR_ld; MR_filt(:,2) = filtfilt(b_filt,a_filt,MR_ld(:,2)); - for i=1:num_f_apply-1 - MR_filt(:,2) = filtfilt(b_filt,a_filt,MR_filt(:,2)); - end end - - hd_mm = 10 * head_ratio_cm; % cm to mm conversion MR_backed = MR_filt; - MR_backed(:, 4:end) = 180*MR_backed(:,4:end)/(pi*hd_mm); %% The last 6 movement regressors are derivatives of the first 6. Needed for task fMRI, but not used in the Fair Lab second_derivs=cat(1, [0 0 0 0 0 0], diff(MR_backed(:,1:6)));