mfbox_jade.m

toolbox for spm 5 for data, model free analysis

function [A,S]=mfbox_jade(X,nem)
% Source separation of complex signals with JADE.
% Jade performs `Source Separation' in the following sense:
%   X is an n x T data matrix assumed modelled as X = A S + N where
% 
% o A is an unknown n x m matrix with full rank.
% o S is a m x T data matrix (source signals) with the properties
%   a) for each t, the components of S(:,t) are statistically
%      independent
%   b) for each p, the S(p,:) is the realization of a zero-mean
%      `source signal'.
%   c) At most one of these processes has a vanishing 4th-order
%      cumulant.
% o  N is a n x T matrix. It is a realization of a spatially white
%    Gaussian noise, i.e. Cov(X) = sigma*eye(n) with unknown variance
%    sigma.  This is probably better than no modeling at all...
%
% Jade performs source separation via a 
% Joint Approximate Diagonalization of Eigen-matrices.  
%
% Input :
%   * X: Each column of X is a sample from the n sensors
%   * nem: nem is an optional argument for the number of sources.
%     If ommited, JADE assumes as many sources as sensors.
%
% Output :
%    * A is an n x nem estimate of the mixing matrix
%    * S is an nem x T naive (ie pinv(A)*X)  estimate of the source signals
%
% Version 1.5.  Copyright: JF Cardoso.  
%
% This file is free software, subject to the 
% GNU GENERAL PUBLIC LICENSE, see gpl.txt

%T=number of rows, n=number of columns
[T,n] = size(X);

if (nargin==1), m=n; end

% sum up all values of each column and calculate mean afterwards
Xmean = mean(X,1);
%subtract mean value
for i=1:n, X(:,i) = X(:,i)-Xmean(i); end

% a statistical threshold for stopping joint diag
threshold = 0.01/sqrt(T);

[E,D,N] = mfbox_pcamat(X',1,nem);
[W,IW] = mfbox_whitenv(E,D,N);

% randomize initial state
R = rand(size(W,1));
R = R*(R'*R)^(-0.5); % orthonormalize
W = R*W;
IW = IW*R';

X = X*W';

[e,m,M] = mfbox_cumindependence(X,true,nem);

V = mfbox_rjd(M,threshold);

A = IW*V;
S = X*V;

Xmean = V'*W*Xmean';
for i=1:size(S,2), S(:,i) = S(:,i)+Xmean(i); end

return