mfbox_stsobi.m

toolbox for spm 5 for data, model free analysis

function [Ss,St]=mfbox_stSOBI(X,varargin)
% Copyright by Peter Gruber, Fabian J. Theis
% Signal Processing & Information Theory group
% Institute of Biophysics, University of Regensburg, Germany
% Homepage: http://research.fabian.theis.name
%           http://www-aglang.uni-regensburg.de
%
% This file is free software, subject to the 
% GNU GENERAL PUBLIC LICENSE, see gpl.txt
%
% spatiotemporal decorrelation (stSOBI)
%
% decomposes the spatiotemporal data set X into
%       X = Ss' * St
% using double-sided joint diagonalization of time / spatially delayed auto-covariances.
% it is an extension of the SOBI algorithm
%
% usage: [Ss,St]=mfbox_stSOBI(X[,argID,value[...]])
%
% with input and output arguments ([]'s are optional):
% X              	(matrix) data matrix X of dimension (ms x mt) respectively of dimension (height x width x mt)
%                            representing sensor (observed) signals, (with ms:=height x width)
%                            so X is a matrix gathering ms spatial observation in the rows each with mt temporal observations (in the columns).
% [argID,           (string) Additional parameters are given as argID, value
%   value]          (varies) pairs. See below for list of valid IDs and values.
% Ss                (matrix) reconstructed spatial source matrix (of dimension n x ms) respectively height x width x n
% St                (matrix) reconstructed temporal source matrix (of dimension n x mt)
%
% Here are the valid argument IDs and corresponding values. 
%  'lastEig' or     (integer) Number n of components to be estimated. Default equals min(ms,mt). It is strongly
%  'numOfIC'                 recommended to set a much smaller n here, otherwise cumulant estimates are very bad. Dimension
%                            reduction is performed using SVD projection along eigenvectors correspoding to largest eigenvalues ('spatiotemporal PCA').
%  'alpha'          (double) weighting parameter in [0,1] (default 0.5). The larger alpha the more temporal separation is favored.
%  'verbose'        (string) report progress of algorithm in text format. 
%                            Either 'on' or 'off'. Default is 'on'.
%  'numOfAutocov'   (integer) Number of autocovariance matrices used for joint diagonalization
%                            Default is 12.
%  'radiusStepSize' (integer) Radius step size to use when calculating autocor. Default is 1.
%  'maskColor'      (integer) If a maskColor is given, all spatial data points with precisely this value are considered
%                            to be transparent and hence not used for autocovariance calculation. 
%                            Note: it is advisable that ALL component of X use the same mask i.e. 
%                            have maskcolor at the same positions.
%                            Default is no mask color i.e. all matrix entries are used.
%  'mask'           (matrix) Logical containing true at voxels where a data
%                            point is and false otherwise, also determines
%                            the data dimension
%  'onedSOBI'       (string) Use one-dimensional SOBI instead of multidimensional one 
%                            (i.e. use single-dim. autocovariance calculation). Either 'on' or 'off'. Default is 'off'.
%
% ------------------------------------------------------------------------------------------------
% REFERENCES:
%   F.J. Theis, P. Gruber, I. Keck, A.-M. Tome and E.W. Lang, 
%   'A spatiotemporal second-order algorithm for fMRI data analysis', 
%   submitted to CIMED 2005.
%
% spatiotemporal ICA (using entropies) has first been introduced by
%   J.V. Stone, J. Porrill, N.R. Porter, and I.W. Wilkinson, 'Spatiotemporal independent component 
%   analysis of event-related fmri data using skewed probability density functions',
%   NeuroImage, 15(2):407-421, 2002.
%
% multidimensional SOBI (SOBI for images and 3d-scans) is defined in
%   F.J. Theis, A. Meyer-Baese and E.W. Lang, 'Second-order blind source 
%   separation based on multi-dimensional autocovariances', in Proc. of ICA 2004 (Granada), 2004.
%
% the classical (one-dimensional) SOBI algorithm is described in
%   A. Belouchrani, K. Abed-Meraim, J.-F. Cardoso and E. Moulines, 'Second-order
%   blind separation of temporally correlated sources', in Proc. Int. Conf. on
%   Digital Sig. Proc., (Cyprus), pp. 346-351, 1993.
%  
% for joint diagonalization stJADE uses Cardoso's diagonalization algorithm based on iterative 
% Given's rotations (matlab-file mfbox_RJD) in the case of orthogonal searches, see
%   J.-F. Cardoso and A. Souloumiac, 'Jacobi angles for simultaneous diagonalization',
%   SIAM J. Mat. Anal. Appl., vol 17(1), pp. 161-164, 1995
% and Yeredors ACDC in the general case of non-orthogonal search, see
%   A. Yeredor, 'Non-Orthogonal Joint Diagonalization in the Least-Squares Sense with Application 
%   in Blind Source Separation, IEEE Trans. On Signal Processing, vol. 50 no. 7 pp. 1545-1553, July 2002.

error(nargchk(1,Inf,nargin))     % check no. of input args

s = size(X);
orgs = s;
mt = orgs(end);
ms = prod(orgs(1:(end-1)));
oned = false; % oned SOBI
verbose = 'on';
n = min(ms,mt);
alpha = 0.5;
mult = 1; % enlarge delta by factor mult - if p <= 12 it makes sense to use mult = 2 or even 3, 
usemask = false; % use mask for spatial data?
nummat = 12;
mask = [];
symmetric = 'on';

for i=1:2:length(varargin),
    %% Check that all argument types are strings
    if ~ischar(varargin{i}),
        error('Invalid input identifier names or input argument order.');
    end
  
    %% Lower/uppercase in identifier types doesn't matter: 
    identifier = lower(varargin{i});     % identifier (lowercase)
    value = varargin{i+1};
        
    switch identifier        
    case 'verbose'
        verbose = lower(value);
    case {'lasteig','numofic'}
        if (value > n)
            error('lastEig resp. numOfIC too large: Can only estimate maximally as many sources as observations!');
        end
        n = value;
    case 'alpha'
        alpha = value;
    case 'numofautocov'
        nummat = value;
    case 'onedsobi'
        oned = strcmp(lower(value),'on');
    case 'radiusstepsize'
        mult = value;
    case 'maskcolor'
        usemask = true;
        maskcolor = value; % this is the color in X only
    case 'mask'
        usemask = true;
        mask = value;
    case 'symmetric'
        symmetric = lower(value);
    otherwise
        %%% Unknown identifier
        error(['Invalid argument identifier ' identifier '!']);        
    end
end

if (strcmp(verbose,'on')), fprintf('stSOBI: centering data...'); end

X = reshape(X,[],mt);

if (usemask)
    if (~isempty(mask))
        sm = size(mask);
        sm = sm(sm>1);
        orgs = [sm,mt];
        s = orgs;
    else
        mask = (~any(X'==maskcolor));
        X = X(mask,:);
    end
else
    mask = true(1,N);
    X = X(mask,:);
end

if (oned), s = [prod(s(1:(end-1))),s(end)]; end

sms = s(1:(end-1));
if (length(sms)==1), sms = [sms,1]; end
mask = reshape(mask,sms);

% remove the spatiotemporal mean
means = mean(X,1)';
for i=1:size(X,2), X(:,i) = X(:,i)-means(i); end
meant = mean(X,2);
for i=1:size(X,1), X(i,:) = X(i,:)-meant(i); end

if (strcmp(verbose,'on'))
    fprintf('done\nstSOBI: performing SVD for spatiotemporal PCA and dimension reduction...');
end

% SVD of X and dimension reduction
[U,D,V] = svd(X,'econ');
%norm(X-U*D*V')
Dx = diag(D);
[Dx,p] = sort(Dx,'descend');
D = D(p,p);
U = U(:,p);
V = V(p,:);
if (n+1<=length(Dx))
    vx = mean(Dx(n+1:end));
    D = diag(Dx(1:n)-vx);
else
    D = diag(Dx(1:n));
end
U = U(:,1:n);
V = V(:,1:n);
%norm(X-U*D*V')

if (strcmp(verbose,'on'))
    fprintf('done\n');
    if (min(ms,mt)>n)
        fprintf('stSOBI: dimension reduction from min(%i, %i) to %i - ',ms,mt,n);
        fprintf('retained %4.2f%% of all eigenvalues\n',sum(Dx(1:n))/sum(Dx)*100);
    else
        fprintf('stSOBI: warning - no dimension reduction has been performed. It is strongly recommended to set a much smaller n using parameter ''lastEig'', otherwise cumulant estimates are very bad.');
    end
end

if strcmp(verbose,'on')
    fprintf('stSOBI: calculating 2 x %i sobi matrices...', nummat);
end
M = zeros(n,n,2*nummat);
tmat = 1:nummat;
smat = nummat+(1:nummat);

sD = sqrt(diag(D));
if (strcmp(symmetric,'on'))
    [c,C,M(:,:,tmat)] = mfbox_sobimats(V*diag(sD),NaN,nummat);
    [c,C,M(:,:,smat)] = mfbox_sobimats({U*diag(sD),mask},nummat);
else
    [c,C,M(:,:,tmat)] = mfbox_sobimats(V*diag(sD.^-1),NaN,nummat);
    [c,C,M(:,:,smat)] = mfbox_sobimats({U*diag(sD.^3),mask},nummat);
end
for i=smat, M(:,:,i) = (M(:,:,i)^(-1)); end

% normalization within the groups in order to allow for comparisons using alpha
sqrs = sqrt(sum(reshape(M,[],2*nummat).^2,1));
M(:,:,tmat) = 2*alpha*M(:,:,tmat)/mean(sqrs(tmat));
M(:,:,smat) = 26*(1-alpha)*M(:,:,smat)/mean(sqrs(smat));
% also possible: normalization by mean determinant
% 3d move -> 26 possibilities
% 1d move -> 2 possibilities

if (strcmp(verbose,'on'))
    fprintf('done\nstSOBI: performing approximate joint diagonalization of %i (%ix%i)-matrices ',2*nummat,n,n);
end

if strcmp(verbose,'on'), fprintf('using orthogonal JD...'); end
W = mfbox_rjd(M)'; % orthogonal approximate JD
if (strcmp(verbose,'on')), fprintf('done\n'); end

if (strcmp(symmetric,'on'))
    St = W*diag(sD)*V';
    Ss = (U*diag(sD)*W^(-1))';
else
    St = W*diag(sD.^-1)*V';
    Ss = (U*diag(sD.^3)*W^(-1))';
end

% add transformed means
Smeant = pinv(Ss)'*meant;
Smeans = pinv(St)'*means;

St = St+repmat(Smeant,1,mt);
tSs = Ss+repmat(Smeans,1,ms);
Ss = zeros([n,orgs(1:(end-1))]);
Ss(:,mask) = tSs;

return