Removes low frequency confounds X0
FORMAT [Y] = pr_spm_filter(K,Y)
FORMAT [K] = pr_spm_filter(K)
K - filter matrix or:
K(s) - struct array containing partition-specific specifications
K(s).RT - observation interval in seconds
K(s).row - row of Y constituting block/partition s
K(s).LChoice - Low-pass filtering {'hrf' 'Gaussian' 'none'}
K(s).LParam - Gaussian parameter in seconds
K(s).HParam - cut-off period in seconds
K(s).X0 - low frequencies to be removed (DCT)
Y - data matrix
K - filter structure
Y - filtered data
___________________________________________________________________________
pr_spm_filter implements high-pass filtering in an efficient way by
using the residual forming matrix of X0 - low frequency confounds
.pr_spm_filter also configures the filter structure in accord with the
specification fields if called with one argument
___________________________________________________________________________
@(#)pr_spm_filter.m 2.10 Karl Friston 03/03/04
0001 function [vargout] = pr_spm_filter(K,Y) 0002 % Removes low frequency confounds X0 0003 % FORMAT [Y] = pr_spm_filter(K,Y) 0004 % FORMAT [K] = pr_spm_filter(K) 0005 % 0006 % K - filter matrix or: 0007 % K(s) - struct array containing partition-specific specifications 0008 % 0009 % K(s).RT - observation interval in seconds 0010 % K(s).row - row of Y constituting block/partition s 0011 % K(s).LChoice - Low-pass filtering {'hrf' 'Gaussian' 'none'} 0012 % K(s).LParam - Gaussian parameter in seconds 0013 % K(s).HParam - cut-off period in seconds 0014 % 0015 % K(s).X0 - low frequencies to be removed (DCT) 0016 % 0017 % Y - data matrix 0018 % 0019 % K - filter structure 0020 % Y - filtered data 0021 %___________________________________________________________________________ 0022 % 0023 % pr_spm_filter implements high-pass filtering in an efficient way by 0024 % using the residual forming matrix of X0 - low frequency confounds 0025 %.pr_spm_filter also configures the filter structure in accord with the 0026 % specification fields if called with one argument 0027 %___________________________________________________________________________ 0028 % @(#)pr_spm_filter.m 2.10 Karl Friston 03/03/04 0029 0030 0031 % set or apply 0032 %--------------------------------------------------------------------------- 0033 if nargin == 1 & isstruct(K) 0034 0035 % set K.X0 0036 %------------------------------------------------------------------- 0037 for s = 1:length(K) 0038 0039 % make high pass filter 0040 %----------------------------------------------------------- 0041 k = length(K(s).row); 0042 n = fix(2*(k*K(s).RT)/K(s).HParam + 1); 0043 X0 = spm_dctmtx(k,n); 0044 K(s).X0 = X0(:,2:end); 0045 0046 % make low pass filter 0047 %----------------------------------------------------------- 0048 if isfield(K(s), 'LChoice') 0049 switch K(s).LChoice 0050 0051 case 'none' 0052 %--------------------------------------------------- 0053 h = 1; 0054 d = 0; 0055 0056 case 'hrf' 0057 %--------------------------------------------------- 0058 h = spm_hrf(K(s).RT); 0059 h = [h; zeros(size(h))]; 0060 g = abs(fft(h)); 0061 h = real(ifft(g)); 0062 h = fftshift(h)'; 0063 n = length(h); 0064 d = [1:n] - n/2 - 1; 0065 0066 case 'Gaussian' 0067 %--------------------------------------------------- 0068 sigma = K(s).LParam/K(s).RT; 0069 h = round(4*sigma); 0070 h = exp(-[-h:h].^2/(2*sigma^2)); 0071 n = length(h); 0072 d = [1:n] - (n + 1)/2; 0073 if n == 1, h = 1; end 0074 0075 otherwise 0076 %--------------------------------------------------- 0077 error('Low pass Filter option unknown'); 0078 return 0079 end 0080 % create and normalize low pass filter 0081 %----------------------------------------------------------- 0082 K(s).KL = spdiags(ones(k,1)*h,d,k,k); 0083 K(s).KL = spdiags(1./sum(K(s).KL')',0,k,k)*K(s).KL; 0084 0085 end 0086 end 0087 0088 % return structure 0089 %------------------------------------------------------------------- 0090 vargout = K; 0091 0092 else 0093 % apply 0094 %------------------------------------------------------------------- 0095 if isstruct(K) 0096 0097 % ensure requisite feilds are present 0098 %----------------------------------------------------------- 0099 if ~isfield(K(1),'X0') | ... 0100 (isfield(K(1),'LChoice') & ~isfield(K(1), 'KL')) 0101 K = pr_spm_filter(K); 0102 end 0103 0104 for s = 1:length(K) 0105 0106 % select data 0107 %--------------------------------------------------- 0108 y = Y(K(s).row,:); 0109 0110 % apply low pass filter 0111 %--------------------------------------------------- 0112 if isfield(K(s), 'KL') 0113 y = K(s).KL*y; 0114 end 0115 0116 % apply high pass filter 0117 %--------------------------------------------------- 0118 y = y - K(s).X0*(K(s).X0'*y); 0119 0120 % reset filtered data in Y 0121 %--------------------------------------------------- 0122 Y(K(s).row,:) = y; 0123 0124 end 0125 0126 % K is simply a filter matrix 0127 %------------------------------------------------------------------- 0128 else 0129 Y = K*Y; 0130 end 0131 0132 % return filtered data 0133 %------------------------------------------------------------------- 0134 vargout = Y; 0135 0136 end 0137 0138