% This code illustrates the idea of multidimensional independent component analysis
% applied to the fetal ECG extraction problem as introduced in the paper:
%
%   JF Cardoso, Multidimensional independent component analysis.  Proc. ICASSP'98. Seattle
%   (available as by anonymous FTP at "ftp://sig.enst.fr/pub/jfc/Papers/icassp98.ps")
%

%  However, in the ICASSP paper, due to lack of room, the analysis was restricted to
%  processing the ouputs of 3 sensors out of 8 in the available data set.  In this demo,
%  we process the whole data set: a better decomposition is achieved.
%  
%  More things are available at http://sig.enst.fr/~cardoso
%
clear


%%  Getting data
disp('Reading data from disk');
load 	foetal_ecg.dat

X	= foetal_ecg(:,2:9)';
Time	= foetal_ecg(:,1  )';

[m,T]	= size(X) ;


disp('Running Jade on all  sensors');
flops(0)
B	= jadeR(X,m);  %%%% Assumes version 1.5 of JADE
Flops	= flops ;
disp(['JADE ran in ' int2str(Flops) ' flops' ]) ;


A	= pinv(B) ;   %% Estimate of the mixing matrix
S	= B*X ;       %% Estimate of the source signals


%% who is where ?  
mother	= [ 1 2 3 4] ;
baby	= [ 7 8  ] ;
noise   = [ 5 6  ] ;

AM	= A(:,mother) ;
AB	= A(:,baby) ;
AN	= A(:,noise) ;

XM 	= A(:,mother)*S(mother,:);
XB	= A(:,baby  )*S(baby  ,:);
XN	= A(:,noise  )*S(noise  ,:);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


disp('Mixing matrix estimated by JADE:') 
A

disp('Orthogonal projection matrix onto 2D mother subspace:') 
Pm	= A(:,mother)*pinv(A(:,mother))

disp('Orthogonal projection matrix onto 1D fetus subspace:') ;
Pb = A(:,baby)*pinv(A(:,baby))

disp('Oblique projection matrix onto 2D mother subspace')
disp(' orthogonally to 1D babysubspace:') 
PTm 	= A(:,mother)*B(mother,:)

disp('Oblique projection matrix onto 1D baby subspace')
disp(' orthogonally to 2D mother subspace:') 
PTb 	= A(:,baby)  *B(baby,  :)



%%% ----------  Fig.1 :  The complete decomposition at sensor level  -----------------


%% The scales used for display are manually chosen:
Ska = [	-50 	50	; 
	-50 	120 	; 
	-100	50	;
	-50	30	;
	-100	50 	;
	-1000	500	;
	-500	1000	;
	-500	1000	] ;  

figure(1); clf

T	= 500 ;  %% Displaying only a fraction of the data set.
Subtime	= 1:T ;
Time	= Time(Subtime);

ncols	= 6 ;

for cap=1:m
	subplot(m,ncols,(cap-1)*ncols+1); 
		plot(Time,X(cap,Subtime));
		axis([ Time(1) Time(T) Ska(cap,1) Ska(cap,2)] );
 		set(gca,'Ytick',[  Ska(cap,1)  0  Ska(cap,2) ])
		set(gca,'Fontsize',6)
		set(gca,'FontName','Helvetica')
		if cap ==1, title('Sensor signals'); end

	subplot(m,ncols,(cap-1)*ncols+2); 
		plot(Time,S(cap,Subtime));
	 	ca = axis; axis([ Time(1) Time(T) ca(3) ca(4) ] );
 		set(gca,'Xtick',[])
		set(gca,'Fontsize',6)
		set(gca,'FontName','Helvetica')
		if cap ==1, title('Source signals'); end

	subplot(m,ncols,(cap-1)*ncols+3); 
 		plot(Time,XM(cap,Subtime))
 		set(gca,'Xtick',[])
 		set(gca,'Ytick',[  Ska(cap,1)  0  Ska(cap,2) ])
		axis([ Time(1) Time(T) Ska(cap,1) Ska(cap,2)] );
		set(gca,'Fontsize',6)
		set(gca,'FontName','Helvetica')
		if cap ==1, title('4D Mother MICA component'); end

	subplot(m,ncols,(cap-1)*ncols+4); 
 		plot(Time,XB(cap,Subtime))
 		set(gca,'Xtick',[])
 		set(gca,'Ytick',[  Ska(cap,1)  0  Ska(cap,2) ])
		axis([ Time(1) Time(T) Ska(cap,1) Ska(cap,2)] );
		set(gca,'FontName','Helvetica')
		set(gca,'Fontsize',6)
		if cap ==1, title('2D Foetal MICA component'); end

	subplot(m,ncols,(cap-1)*ncols+5); 
 		plot(Time,XB(cap,Subtime))
 		set(gca,'Xtick',[])
% 		set(gca,'Ytick',[  Ska(cap,1)  0  Ska(cap,2) ])
%		axis([ Time(1) Time(T) Ska(cap,1) Ska(cap,2)] );
		set(gca,'FontName','Helvetica')
		set(gca,'Fontsize',6)
		if cap ==1, title('Full scale Foetal MICA component'); end

		
	subplot(m,ncols,(cap-1)*ncols+6); 
 		plot(Time,XN(cap,Subtime))
 		set(gca,'Xtick',[])
 		set(gca,'Ytick',[  Ska(cap,1)  0  Ska(cap,2) ])
		axis([ Time(1) Time(T) Ska(cap,1) Ska(cap,2)] );
		set(gca,'FontName','Helvetica')
		set(gca,'Fontsize',6)
		if cap ==1, title('2D Noise MICA component'); end

		
		
end


set(gcf,'PaperOrientation','Landscape')
set(gcf,'Paperunits','centimeters')
% set(gcf,'PaperPosition',[ 2 2 22 8])