% A simple demo program for EASI :
% an adaptive blind source separation algorithm
% Ref :
%   @inproceedings{PFSEusipco,
%   title = "Adaptive source separation without prewhitening",
%   author = "Beate Laheld and Jean-Fran\c{c}ois Cardoso",
%   booktitle = "{Proc. EUSIPCO'94}",
%   pages = "183--186",
%   address = {Edinburgh},
%   month = sep,
%   year = "1994" }
% For any information, comments, etc... contact cardoso@sig.enst.fr
%--------------------------------------------------------------------

%
% Aug. 22, 97: corrected a bug where a comment was not commented out...
clear
blah=[
'Each experiment consists in :                                       ';
' - Randomly choosing a mixing matrix and an intitial separator      ';
' - Running 200 steps of EASI with QAM4 signals                      ';
' - Displaying the coefficients of the global system separator*mixing';
'    These should converge to 0 or 1                                 ';
];
disp(blah);
% -----------------------------------------------------------------------

nsou  = 3;	% number of sources
ncap  = 4;	% number of sensors (not smaller than nsou !)

T	= 200;	% number of samples
lambda  = 0.02;	% adaptation step

idsou	= eye(nsou)		;
suivi	= zeros(nsou*nsou,T)	; % a buffer 

while 1, fprintf('Another experiment\n');

A= randn(ncap,nsou)+i*randn(ncap,nsou);	% random mixing matrix
B= randn(nsou,ncap)+i*randn(nsou,ncap); % random initial separating matrix

for t=1:T

 % next line generates a vector of 'nsou' independent normalized QAM4 signals
 s	= (2*fix(2*rand(nsou,1))-1+i* (2*fix(2*rand(nsou,1))-1)) /sqrt(2);
 x	= A*s		; 	% the `observed'  mixture
 y	= B*x		;	% the output

%-----------------  here is the algorithm   ----------
 y2	= y*y'			;
 g	= diag(y2).*y		;
% Other nonlinearities to be tried............
% g	=   y .* diag(0.1*y2).*diag(y2)	;
% g	=   y .* sqrt(diag(y2))		;
% g	=   y .*  log(diag(y2))		;
% g	= - y .*     (diag(y2)<0.9)	;
% g	=   y ./  log(diag(y2))		;
% g	= - y ./ sqrt(diag(y2))		;
% g	= - y ./      diag(y2)		;
%----
 gy	= g * y'		;
 G	= (y2 - idsou)/(1+lambda*trace(y2)) + (gy - gy')/(1+lambda*abs(g'*y)) ;
 B	= B - lambda*G*B	; % updating the separating matrix
%-----------------  here ends the algorithm ----------

 tot	= B*A		;  % the global mixing-unmixing system 
 suivi(:,t)= tot(:)	;  % keep for further display

end;

%--- display the modulus of the coefficients of the global system
	clg;
	plot( [0 T], [0 0]);plot( [0 T], [ 1 1 ]);
	xlabel('Iteration number');
	ylabel('Modulus of mixture coefficients');
	hold on;
	for ip=1:nsou*nsou
	 plot(abs(suivi(ip,1:T)));
	end;
	hold off;
	axis([0 T 0 2]);
drawnow;

end;	% next experiment