% This code can be used to reproduce the fetal ECG extraction experiment published 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"
%  
%  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  )';


%% Selecting the 3 first sensors
select	= [ 1 2 3 ] ;
X	= X(select',:);
[m,T]	= size(X) ;



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


%% About the next two lines: We change the order and the signs of the source signals.
%% This is because the version of JADE (jadeR.m versin 1.5) bundled with this demo program
%% uses a different sorting convention than the one used when the ICASSP paper was
%% written.  Of course, this does not make any difference in the end so the next two lines
%% could be commented out.  However, we artificially restore the original order and signs
%% in order to produce the very same figure as the one included in the paper.  Actually,
%% this is not true strictly speaking: a careful reader will notice (unsignificant)
%% differences between the numerical values produced by this code and those reported in
%% the ICASSP paper.  This is due to the fact that the current version of JADE removes the
%% sample mean from the signals by default.  Since the data set is approximately
%% zero-mean, this makes very little difference.
B       = B([ 3 2 1 ], :);       %% changing the order of the sources
B       = diag([1 -1 1]) * B ;   %% changing the sign of the second source


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


%% who is where ?  
baby	= [ 1 ] ;
mother	= [ 2 3 ] ;


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

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


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


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,  :)


disp(' ')
disp('Note: There are minute differences between the above values and')
disp('those reported in the paper.  They are not significant (compare the ')
disp('estimated signals).  They are due to the fact that the mean value')
disp('of the signals is removed in the current implementation of JADE. ')


%%% ----------  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	] ;  

Ska	= Ska(select',:);


figure(1); clf

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

ncols	= 4 ;

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',9)
		set(gca,'FontName','Times')
		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',9)
		set(gca,'FontName','Times')
		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',9)
		set(gca,'FontName','Times')
		if cap ==1, title('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','Times')
		set(gca,'Fontsize',9)
		if cap ==1, title('Foetal MICA component'); end

end


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