contruct and/or apply high and/or low pass filter
Copied from spm_filter.m for SPM99.
FORMAT [K] = spm_filter('set' ,K)
FORMAT [Y] = spm_filter('apply',K,Y)
FORMAT [Y] = spm_filter('high' ,K,Y)
FORMAT [Y] = spm_filter('low' ,K,Y)
Action - 'set' fills in filter structure K
Action - 'apply' applies K to Y = K*Y
Action - 'high' only high-pass component
Action - 'low ' only low-pass component
K - filter convolution matrix or:
K{s} - cell of structs containing session-specific specifications
K{s}.RT - repeat time in seconds
K{s}.row - row of Y constituting session s
K{s}.LChoice - Low-pass filtering {'hrf' 'Gaussian' 'none'}
K{s}.LParam - Gaussian parameter in seconds
K{s}.HChoice - High-pass filtering {'specify' 'none'}
K{s}.HParam - cut-off period in seconds
K{s}.HP - low frequencies to be removed
K{s}.LP - sparse toepltz low-pass convolution matrix
Y - data matrix
K - filter structure
Y - filtered data K.K*Y
___________________________________________________________________________
spm_filter implements band pass filtering in an efficient way by
using explicitly the projector matrix form of the High pass
component. spm_filter also configures the filter structure in
accord with the specification fields if required
___________________________________________________________________________
@(#)spm_filter.m 2.5 Karl Friston 00/10/04
0001 function [vargout] = pr_spm_filter(Action,K,Y) 0002 % contruct and/or apply high and/or low pass filter 0003 % Copied from spm_filter.m for SPM99. 0004 % FORMAT [K] = spm_filter('set' ,K) 0005 % FORMAT [Y] = spm_filter('apply',K,Y) 0006 % FORMAT [Y] = spm_filter('high' ,K,Y) 0007 % FORMAT [Y] = spm_filter('low' ,K,Y) 0008 % 0009 % Action - 'set' fills in filter structure K 0010 % Action - 'apply' applies K to Y = K*Y 0011 % Action - 'high' only high-pass component 0012 % Action - 'low ' only low-pass component 0013 % K - filter convolution matrix or: 0014 % K{s} - cell of structs containing session-specific specifications 0015 % 0016 % K{s}.RT - repeat time in seconds 0017 % K{s}.row - row of Y constituting session s 0018 % K{s}.LChoice - Low-pass filtering {'hrf' 'Gaussian' 'none'} 0019 % K{s}.LParam - Gaussian parameter in seconds 0020 % K{s}.HChoice - High-pass filtering {'specify' 'none'} 0021 % K{s}.HParam - cut-off period in seconds 0022 % 0023 % K{s}.HP - low frequencies to be removed 0024 % K{s}.LP - sparse toepltz low-pass convolution matrix 0025 % 0026 % Y - data matrix 0027 % 0028 % K - filter structure 0029 % Y - filtered data K.K*Y 0030 %___________________________________________________________________________ 0031 % 0032 % spm_filter implements band pass filtering in an efficient way by 0033 % using explicitly the projector matrix form of the High pass 0034 % component. spm_filter also configures the filter structure in 0035 % accord with the specification fields if required 0036 %___________________________________________________________________________ 0037 % @(#)spm_filter.m 2.5 Karl Friston 00/10/04 0038 0039 0040 % set or apply 0041 %--------------------------------------------------------------------------- 0042 switch Action 0043 0044 case 'set' 0045 %------------------------------------------------------------------- 0046 for s = 1:length(K) 0047 0048 % matrix order 0049 %----------------------------------------------------------- 0050 k = length(K{s}.row); 0051 0052 % make low pass filter 0053 %----------------------------------------------------------- 0054 switch K{s}.LChoice 0055 0056 case 'none' 0057 %--------------------------------------------------- 0058 h = 1; 0059 d = 0; 0060 0061 case 'hrf' 0062 %--------------------------------------------------- 0063 h = pr_spm_hrf(K{s}.RT); 0064 h = [h; zeros(size(h))]; 0065 g = abs(fft(h)); 0066 h = real(ifft(g)); 0067 h = fftshift(h)'; 0068 n = length(h); 0069 d = [1:n] - n/2 - 1; 0070 0071 case 'Gaussian' 0072 %--------------------------------------------------- 0073 sigma = K{s}.LParam/K{s}.RT; 0074 h = round(4*sigma); 0075 h = exp(-[-h:h].^2/(2*sigma^2)); 0076 n = length(h); 0077 d = [1:n] - (n + 1)/2; 0078 if n == 1, h = 1; end 0079 0080 otherwise 0081 %--------------------------------------------------- 0082 error('Low pass Filter option unknown'); 0083 return 0084 0085 end 0086 0087 % create and normalize low pass filter 0088 %----------------------------------------------------------- 0089 K{s}.KL = spdiags(ones(k,1)*h,d,k,k); 0090 K{s}.KL = spdiags(1./sum(K{s}.KL')',0,k,k)*K{s}.KL; 0091 0092 0093 % make high pass filter 0094 %----------------------------------------------------------- 0095 switch K{s}.HChoice 0096 0097 case 'none' 0098 %--------------------------------------------------- 0099 K{s}.KH = []; 0100 0101 case 'specify' 0102 %--------------------------------------------------- 0103 n = fix(2*(k*K{s}.RT)/K{s}.HParam + 1); 0104 X = spm_dctmtx(k,n); 0105 K{s}.KH = sparse(X(:,[2:n])); 0106 0107 otherwise 0108 %--------------------------------------------------- 0109 error('High pass Filter option unknown'); 0110 return 0111 0112 end 0113 0114 end 0115 0116 % return structure 0117 %------------------------------------------------------------------- 0118 vargout = K; 0119 0120 0121 case {'apply','high','low'} 0122 %------------------------------------------------------------------- 0123 if iscell(K) 0124 0125 0126 % ensure requisite feild are present 0127 %----------------------------------------------------------- 0128 if ~isfield(K{1},'KL') 0129 K = pr_spm_filter('set',K); 0130 end 0131 0132 for s = 1:length(K) 0133 0134 % select data 0135 %--------------------------------------------------- 0136 y = Y(K{s}.row,:); 0137 0138 % apply low pass filter 0139 %--------------------------------------------------- 0140 if ~strcmp(Action,'high') 0141 KL = K{s}.KL; 0142 if issparse(KL), KL = full(KL); end 0143 y = KL*y; 0144 end 0145 0146 % apply high pass filter 0147 %--------------------------------------------------- 0148 if ~isempty(K{s}.KH) & ~strcmp(Action,'low') 0149 KH = K{s}.KH; 0150 if issparse(KH), KH = full(KH); end 0151 y = y - KH*(KH'*y); 0152 end 0153 0154 % reset filtered data in Y 0155 %--------------------------------------------------- 0156 Y(K{s}.row,:) = y; 0157 0158 end 0159 0160 % K is simply a convolution matrix 0161 %------------------------------------------------------------------- 0162 else 0163 if issparse(K), K = full(K); end 0164 Y = K*Y; 0165 end 0166 0167 % return filtered data 0168 %------------------------------------------------------------------- 0169 vargout = Y; 0170 0171 case 'none' 0172 %------------------------------------------------------------------- 0173 vargout = Y; 0174 0175 otherwise 0176 %------------------------------------------------------------------- 0177 warning('Filter option unknown'); 0178 0179 0180 end