mfbox_stjade.m

toolbox for spm 5 for data, model free analysis

function [Ss,St]=mfbox_stJADE(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 JADE (stJADE)
%
% decomposes the spatiotemporal data set X into
%       X = Ss' * St
% using double-sided joint diagonalization of fourth-order cumulants.
% it is an extension of the Cardoso's JADE algorithm
%
% usage: [Ss,St]=mfbox_stJADE(X[,argID,value[...]])
%
% with input and output arguments ([]'s are optional):
% X              	(matrix) data matrix X of dimension (ms x mt) representing sensor (observed) signals,
%                            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)
% 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'.
%
% ------------------------------------------------------------------------------------------------
% REFERENCES:
%   F.J. Theis, P. Gruber, I. Keck, A. Meyer-Baese and E.W. Lang, 
%   'Spatiotemporal blind source separation using double-sided approximate joint diagonalization', 
%   submitted to EUSIPCO 2005 (Antalya), 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.
%
% the classical ('temporal') JADE algorithm is described in
%   J.-F. Cardoso and A. Souloumiac, 'Blind beamforming for non
%   gaussian signals', IEE Proceedings - F, 140(6):362-370, 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

[ms,mt] = size(X);

verbose = 'on';
n = min(ms,mt);
alpha = 0.5;
nummat = Inf;
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 'nummat'
        nummat = value;
    case 'symmetric'
        symmetric = lower(value);
    otherwise
        %%% Unknown identifier
        error(['Invalid argument identifier ''' identifier '''!']);        
    end
end

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

% remove the spatiotemporal mean
means = mean(X,1)';
meant = mean(X,2);
X = X-repmat(means',ms,1);
X = X-repmat(mean(X,2),1,mt);

if (strcmp(verbose,'on'))
    fprintf('done\nstJADE: 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('stJADE: 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('stJADE: 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

nummat = min(nummat,n*(n+1)/2);
if strcmp(verbose,'on')
    fprintf('stJADE: calculating 2 x %i contracted cumulant 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_cumindependence(V*diag(sD),false,nummat);
    [c,C,M(:,:,smat)] = mfbox_cumindependence(U*diag(sD),false,nummat);
else
    [c,C,M(:,:,tmat)] = mfbox_cumindependence(V*diag(sD.^-1),false,nummat);
    [c,C,M(:,:,smat)] = mfbox_cumindependence(U*diag(sD.^3),false,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) = alpha*M(:,:,tmat)/mean(sqrs(tmat));
M(:,:,smat) = (1-alpha)*M(:,:,smat)/mean(sqrs(smat));
% also possible: normalization by mean determinant

if (strcmp(verbose,'on'))
    fprintf('done\nstJADE: 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);
Ss = Ss+repmat(Smeans,1,ms);

return