algorithm.zip_zip

function B = jadeR(X,m) % Blind separation of real signals with JADE. Version 1.5 Dec. 1997. % % Usage: % * If X is an nxT data matrix (n sensors, T samples) then % B=jadeR(X) is a nxn separating matrix such that S=B*X is an nxT % matrix of estimated source signals. % * If B=jadeR(X,m), then B has size mxn so that only m sources are % extracted. This is done by restricting the operation of jadeR % to the m first principal components. % * Also, the rows of B are ordered such that the columns of pinv(B) % are in order of decreasing norm; this has the effect that the % `most energetically significant' components appear first in the % rows of S=B*X. % % Quick notes (more at the end of this file) % % o this code is for REAL-valued signals. An implementation of JADE % for both real and complex signals is also available from % http://sig.enst.fr/~cardoso/stuff.html % % o This algorithm differs from the first released implementations of % JADE in that it has been optimized to deal more efficiently % 1) with real signals (as opposed to complex) % 2) with the case when the ICA model does not necessarily hold. % % o There is a practical limit to the number of independent % components that can be extracted with this implementation. Note % that the first step of JADE amounts to a PCA with dimensionality % reduction from n to m (which defaults to n). In practice m % cannot be `very large' (more than 40, 50, 60... depending on % available memory) % % o See more notes, references and revision history at the end of % this file and more stuff on the WEB % http://sig.enst.fr/~cardoso/stuff.html % % o This code is supposed to do a good job! Please report any % problem to cardoso@sig.enst.fr % Copyright : Jean-Francois Cardoso. cardoso@sig.enst.fr verbose = 0 ; % Set to 0 for quiet operation % Finding the number of sources [n,T] = size(X); if nargin==1, m=n ; end; % Number of sources defaults to # of sensors if m>n , fprintf('jade -> Do not ask more sources than sensors here!!!\n'), return,end if verbose, fprintf('jade -> Looking for %d sources\n',m); end ; % Self-commenting code %===================== if verbose, fprintf('jade -> Removing the mean value\n'); end X = X - mean(X')' * ones(1,T); %%% whitening & projection onto signal subspace % =========================================== if verbose, fprintf('jade -> Whitening the data\n'); end [U,D] = eig((X*X')/T) ; [puiss,k] = sort(diag(D)) ; rangeW = n-m+1:n ; % indices to the m most significant directions scales = sqrt(puiss(rangeW)) ; % scales W = diag(1./scales) * U(1:n,k(rangeW))' ; % whitener iW = U(1:n,k(rangeW)) * diag(scales) ; % its pseudo-inverse X = W*X; %%% Estimation of the cumulant matrices. % ==================================== if verbose, fprintf('jade -> Estimating cumulant matrices\n'); end dimsymm = (m*(m+1))/2; % Dim. of the space of real symm matrices nbcm = dimsymm ; % number of cumulant matrices CM = zeros(m,m*nbcm); % Storage for cumulant matrices R = eye(m); %% Qij = zeros(m); % Temp for a cum. matrix Xim = zeros(1,m); % Temp Xjm = zeros(1,m); % Temp scale = ones(m,1)/T ; % for convenience %% I am using a symmetry trick to save storage. I should write a %% short note one of these days explaining what is going on here. %% Range = 1:m ; % will index the columns of CM where to store the cum. mats. for im = 1:m Xim = X(im,:) ; Qij = ((scale* (Xim.*Xim)) .* X ) * X' - R - 2 * R(:,im)*R(:,im)' ; CM(:,Range) = Qij ; Range = Range + m ; for jm = 1:im-1 Xjm = X(jm,:) ; Qij = ((scale * (Xim.*Xjm) ) .*X ) * X' - R(:,im)*R(:,jm)' - R(:,jm)*R(:,im)' ; CM(:,Range) = sqrt(2)*Qij ; Range = Range + m ; end ; end; %%% joint diagonalization of the cumulant matrices %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Init if 1, %% Init by diagonalizing a *single* cumulant matrix. It seems to save %% some computation time `sometimes'. Not clear if initialization is %% a good idea since Jacobi rotations are very efficient. if verbose, fprintf('jade -> Initialization of the diagonalization\n'); end [V,D] = eig(CM(:,1:m)); % For instance, this one for u=1:m:m*nbcm, % updating accordingly the cumulant set given the init CM(:,u:u+m-1) = CM(:,u:u+m-1)*V ; end; CM = V'*CM; else, %% The dont-try-to-be-smart init V = eye(m) ; % la rotation initiale end; seuil = 1/sqrt(T)/100; % A statistically significant threshold encore = 1; sweep = 0; updates = 0; g = zeros(2,nbcm); gg = zeros(2,2); G = zeros(2,2); c = 0 ; s = 0 ; ton = 0 ; toff = 0 ; theta = 0 ; %% Joint diagonalization proper if verbose, fprintf('jade -> Contrast optimization by joint diagonalization\n'); end while encore, encore=0; if verbose, fprintf('jade -> Sweep #%d\n',sweep); end sweep=sweep+1; for p=1:m-1, for q=p+1:m, Ip = p:m:m*nbcm ; Iq = q:m:m*nbcm ; %%% computation of Givens angle g = [ CM(p,Ip)-CM(q,Iq) ; CM(p,Iq)+CM(q,Ip) ]; gg = g*g'; ton = gg(1,1)-gg(2,2); toff = gg(1,2)+gg(2,1); theta = 0.5*atan2( toff , ton+sqrt(ton*ton+toff*toff) ); %%% Givens update if abs(theta) > seuil, encore = 1 ; updates = updates + 1; c = cos(theta); s = sin(theta); G = [ c -s ; s c ] ; pair = [p;q] ; V(:,pair) = V(:,pair)*G ; CM(pair,:) = G' * CM(pair,:) ; CM(:,[Ip Iq]) = [ c*CM(:,Ip)+s*CM(:,Iq) -s*CM(:,Ip)+c*CM(:,Iq) ] ; %% fprintf('jade -> %3d %3d %12.8f\n',p,q,s); end%%of the if end%%of the loop on q end%%of the loop on p end%%of the while loop if verbose, fprintf('jade -> Total of %d Givens rotations\n',updates); end %%% A separating matrix % =================== B = V'*W ; %%% We permut its rows to get the most energetic components first. %%% Here the **signals** are normalized to unit variance. Therefore, %%% the sort is according to the norm of the columns of A = pinv(B) if verbose, fprintf('jade -> Sorting the components\n'); end A = iW*V ; [vars,keys] = sort(sum(A.*A)) ; B = B(keys,:); B = B(m:-1:1,:) ; % Is this smart ? % Signs are fixed by forcing the first column of B to have % non-negative entries. if verbose, fprintf('jade -> Fixing the signs\n'); end b = B(:,1) ; signs = sign(sign(b)+0.1) ; % just a trick to deal with sign=0 B = diag(signs)*B ; return ; % To do. % - Implement a cheaper/simpler whitening (is it worth it?) % % Revision history: % %- V1.5, Dec. 24 1997 % - The sign of each row of B is determined by letting the first % element be positive. % %- V1.4, Dec. 23 1997 % - Minor clean up. % - Added a verbose switch % - Added the sorting of the rows of B in order to fix in some % reasonable way the permutation indetermination. See note 2) % below. % %- V1.3, Nov. 2 1997 % - Some clean up. Released in the public domain. % %- V1.2, Oct. 5 1997 % - Changed random picking of the cumulant matrix used for % initialization to a deterministic choice. This is not because % of a better rationale but to make the ouput (almost surely) % deterministic. % - Rewrote the joint diag. to take more advantage of Matlab's % tricks. % - Created more dummy variables to combat Matlab's loose memory % management. % %- V1.1, Oct. 29 1997. % Made the estimation of the cumulant matrices more regular. This % also corrects a buglet... % %- V1.0, Sept. 9 1997. Created. % % Main reference: % @article{CS-iee-94, % title = "Blind beamforming for non {G}aussian signals", % author = "Jean-Fran\c{c}ois Cardoso and Antoine Souloumiac", % HTML = "ftp://sig.enst.fr/pub/jfc/Papers/iee.ps.gz", % journal = "IEE Proceedings-F", % month = dec, number = 6, pages = {362-370}, volume =