mfbox_mdsobi.m

toolbox for spm 5 for data, model free analysis

function varargout=mfbox_mdsobi(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
%
% spatially decorrelates the either 2- or 3-dimensional data set X using multiple spatial autocovariances.
% it is an extension of the classical SOBI algorithm of Belouchrani et al.
% usage:
%       S = mfbox_mdsobi (X [,argID,value [...]])
% or 
%       [A,W] = mfbox_mdsobi (X [,argID,value [...]])
% or
%       [S,A,W] = mfbox_mdsobi (X [,argID,value [...]])
%
% with input and output arguments ([]'s are optional):
% X              	(matrix) data matrix X of dimension (height x width x numberofsignals) respectively 
%                           (height x width x depth x numberof3dscans) representing sensor (observed) signals,
% [argID,           (string) Additional parameters are given as argID, value
%   value]          (varies) pairs. See below for list of valid IDs and values.
% S                 (matrix) reconstructed source matrix
% A                 (matrix) mixing matrix 
% W                 (matrix) unmixing matrix (S=WX)
%
% Here are the valid argument IDs and corresponding values. 
%  'lastEig' or     (integer) Number of components to be estimated. Default equals the dimension of data. As usual, dimension
%  'numOfIC'                 reduction is performed using PCA projection along eigenvectors correspoding to largest eigenvalues.
%  '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 matrix entries of X 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'.
%  'verbose'        (string) report progress of algorithm in text format. 
%                            Either 'on' or 'off'. Default is 'on'.
%
% ------------------------------------------------------------------------------------------------
% REFERENCES:
%   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 mfbox_mdsobi uses Cardoso's diagonalization algorithm based on iterative 
% Given's rotations, see
%   J.-F. Cardoso and A. Souloumiac, 'Jacobi angles for simultaneous diagonalization',
%   SIAM J. Mat. Anal. Appl., vol 17(1), pp. 161-164, 1995

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

s = size(X);
orgs = s;

m = s(end);
N = prod(s(1:(end-1)));

verbose = 'on';
n = m;
p = 12;
mult = 1; % enlarge delta by factor mult - if p <= 12 it makes sense to use mult=2 or even 3, 
        %for higher p small mult has better resolution
oned = 0;
usemask = 0;

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 'numofautocov'
            nummats = value;
        case 'onedsobi'
            oned = strcmpi(value,'on');
        case 'verbose'
            verbose = lower(value);
        case {'lasteig','numofic'}
            n = value;
            if (n>m), error('numOfIC too large: Can only estimate maximally as many sources as observations!'); end
        case 'radiusstepsize'
            mult = value;
        case 'maskcolor'
            usemask = 1;
            maskcolor = value; % this is the color in X only, in Xc and so on the mask will be set to 0
        case 'mask'
            usemask = true;
            mask = value;
        otherwise
            error(['Invalid argument identifier ' identifier '!']);        
    end
end

if (strcmp(verbose,'on'))
    fprintf('mdSOBI: starting calculation\n');
    if (oned), fprintf('using 1d-autocovariances\n');
    else, fprintf('using nd-autocovariances\n');
    end
    fprintf('centering and whitening data...');
end

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

if (usemask)
    if (~isempty(mask))
        sm = size(mask);
        sm = sm(sm>1);
        orgs = [sm,m];
        s = orgs;
        N = prod(s(1:(end-1)));
    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);

meanX = mean(X,1);
for i=1:size(X,2), X(:,i) = X(:,i)-meanX(i); end

[pcaE,pcaD,pcaN] = mfbox_pcamat(X',1,n);
[whiteM,dewhiteM] = mfbox_whitenv(pcaE,pcaD,pcaN);

X = real(X*whiteM');

if (strcmp(verbose, 'on'))
    fprintf('done\n');
    if (m>n)
        fprintf('dimension reduction from %i to %i - ',m,n);
        fprintf('retained %4.2f%% of all eigenvalues\n',sum(pcaD)/(sum(pcaD)+(m-n)*pcaN)*100);
    end
end

[c,C,M] = mfbox_sobimats({X,mask},nummats,mult);

if (strcmp(verbose,'on'))
    fprintf('done\nperforming approximate joint diagonalization of %i (%ix%i)-matrices...',p,n,n);
end
U = mfbox_rjd(reshape(M,n,[]));

W = U'*whiteM;
X = X*U+repmat(meanX*W',sum(mask(:)),1);
S = zeros(N,n);
S(mask,:) = X;

if ((nargout==1)||(nargout==3))
    varargout{1} = reshape(S,[orgs(1:(end-1)),n]); % recovered sources
end
if (nargout==2), varargout{1} = dewhiteM*U; varargout{2} = W; end
if (nargout==3), varargout{2} = dewhiteM*U; varargout{3} = W; end

return