sepa_l.m

盲分离算法

% SOURCE SEPARATION IN LINEAR MIXTURES

%-------------------------------------------------------------------------------
%	s  =	source signals
%	e  = 	observations
%	y  =	output signals
%-------------------------------------------------------------------------------

% Grenoble, December 2000
% This work has been partly funded by the European project BLIS (IST-1999-14190)


% initialisations
iter=1;					% incrementation
nsrc=2;					% number of sources
Id=eye(nsrc);				% identity matrix
B=Id;					% separating matrix
track_B=B(:);				% keep for further display

source1;				% source s1
source2;				% source s2
s=[s1; s2];				% sources signal (statistically independent)
matrix;					% mixing matrix (nsrc-by-nsrc)

% observations
e=A*s;					% linear system

% plots
subplot(ay(1));
plot(e(1,:),'r','MarkerSize',2);
title('observations 1','FontSize',15);
subplot(ay(2));
plot(e(2,:),'b','MarkerSize',2);
title('observations 2','FontSize',15);
subplot(ay(3));
plot(s(1,:),s(2,:),'*','MarkerSize',2);title('sources distribution','FontSize',15);
subplot(ay(4));
plot(e(1,:),e(2,:),'*','MarkerSize',2);title('observations distribution','FontSize',15);

% outputs initialisation
y=B*e;					% outputs of the linear system estimatiom

waitfor(ay(8),'Userdata',1);		% wait user press "START" on GUI

% iterations
while iter <= itmax

   set(az(8),'String',iter);		% set the number of incrementation (iter) on GUI
   
   % score function estimation
   psi=estim_psi(y,sigma);
   
   % linear system estimation
   G=(psi*y')/T;			% E(psi(yj)*yi)
   YY=y*y'/T;				% E(yi^2) (=variance(yi))
   H=G-diag(diag(G))+Id-diag(diag(YY));	% ensure a decrease of the mutual information
   					% and a normalisation (variance(yi)=1)			
   B=(Id+mu*H)*B;			% algorithm
   
   % output signal
   y=B*e;

   % keep for further display
   track_B=[track_B B(:)];
   
   % plots
   subplot(az(2));
   plot(y(1,:),'r','MarkerSize',2);title('output signal 1','FontSize',15);
   subplot(az(3));
   plot(y(2,:),'b','MarkerSize',2);title('output signal 2','FontSize',15);
   subplot(az(1));
   plot(y(1,:),y(2,:),'*','MarkerSize',2);title('outputs distribution','FontSize',15);
   drawnow;
   
   iter=iter+1;				% incrementation
   
end;
	
% plots
;
subplot(az(2));plot(s(1,:),'r','MarkerSize',2);
hold on;plot(y(1,:),'b','MarkerSize',2);title('red: source 1   blue: output 1','FontSize',15);
subplot(az(3));plot(s(2,:),'r','MarkerSize',2);
hold on;plot(y(2,:),'b','MarkerSize',2);title('red: source 2   blue: output 2','FontSize',15);

subplot(az(11));
hold on;
for k=1:nsrc*nsrc
 plot([0:iter-1],track_B(k,:),'MarkerSize',2);
end;
title('separating matrix coefficients','FontSize',15);