function varargout = spm_filter(K,Y)
% Removes low frequency confounds X0
% FORMAT [Y] = spm_filter(K,Y)
% FORMAT [K] = 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).HParam - cut-off period in seconds
%
% K(s).X0 - low frequencies to be removed (DCT)
%
% Y - data matrix
%
% K - filter structure
% Y - filtered data
%__________________________________________________________________________
%
% spm_filter implements high-pass filtering in an efficient way by
% using the residual forming matrix of X0 - low frequency confounds.
% spm_filter also configures the filter structure in accord with the
% specification fields if called with one argument.
%__________________________________________________________________________
% Copyright (C) 1999-2015 Wellcome Trust Centre for Neuroimaging
% Karl Friston
% $Id: spm_filter.m 6416 2015-04-21 15:34:10Z guillaume $
%-Configure filter
%==========================================================================
if nargin == 1 && isstruct(K)
% set K.X0
%----------------------------------------------------------------------
for s = 1:length(K)
% make high pass filter
%------------------------------------------------------------------
k = length(K(s).row);
n = fix(2*(k*K(s).RT)/K(s).HParam + 1);
X0 = spm_dctmtx(k,n);
K(s).X0 = X0(:,2:end);
end
% return filter structure
%----------------------------------------------------------------------
varargout = { K };
%-Apply filter
%==========================================================================
else
% K is a filter structure
%----------------------------------------------------------------------
if isstruct(K)
% ensure requisite fields are present
%------------------------------------------------------------------
if ~isfield(K(1),'X0')
K = spm_filter(K);
end
if numel(K) == 1 && length(K.row) == size(Y,1)
% apply high pass filter
%--------------------------------------------------------------
Y = Y - K.X0*(K.X0'*Y);
else
for s = 1:length(K)
% select data
%----------------------------------------------------------
y = Y(K(s).row,:);
% apply high pass filter
%----------------------------------------------------------
y = y - K(s).X0*(K(s).X0'*y);
% reset filtered data in Y
%----------------------------------------------------------
Y(K(s).row,:) = y;
end
end
% K is simply a filter matrix
%----------------------------------------------------------------------
else
Y = K * Y;
end
% return filtered data
%----------------------------------------------------------------------
varargout = { Y };
end