sepa_pnl_par.m

盲分离算法

% SOURCE SEPARATION IN NONLINEAR POST MIXTURES (PARAMETRIC METHOD)

%-------------------------------------------------------------------------------
%	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
a=zeros(nsrc,N);			% parameter
a(:,2)=1;
source1;				% source s1
source2;				% source s2
s=[s1; s2];				% sources signal (statistically independent)
matrix;					% mixing matrix (nsrc-by-nsrc)

% observations 
tu=A*s;					% linear subsystem
e(1,:)=f1(tu(1,:));			% nonlinear subsystem (f1 : nonlinear function)
e(2,:)=f2(tu(2,:));			% nonlinear subsystem (f2 : nonlinear function)

% plots
subplot(ay(10))
plot(tu(1,:),e(1,:),'*r','MarkerSize',2);title('nonlinear function 1','FontSize',15);
subplot(ay(9))
plot(tu(2,:),e(2,:),'*','MarkerSize',2);title('nonlinear function 2','FontSize',15);
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
x=e;					% output of the nonlinear subsystem estimation
y=B*x;					% output of the linear subsystem estimation

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);
   
   % nonlinear subsystem estimation
   alpha=B'*psi;			% matrix (2-by-T) sum(psi(yi)*bik)
   eps=eps_pnl_par(alpha,e,a);		% ensure a decrease of the mutual information
   a=a+nu*eps;				% algorithm	
   for t=1:nsrc
   	x(t,:)=polyval(fliplr(a(t,:)),e(t,:));	% calcul of the output
   end;
   	
   % normalisation (scale and shift factor indeterminacies)
   x(1,:)=x(1,:)-mean(x(1,:));
   x(2,:)=x(2,:)-mean(x(2,:));
   a=inv(diag([std(x(1,:)) std(x(2,:))]))*a;
   x(1,:)=x(1,:)/std(x(1,:));
   x(2,:)=x(2,:)/std(x(2,:));
	
   % 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*x;

  

  % 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);
   subplot(az(5));
   plot(tu(1,:),x(1,:),'*r','MarkerSize',2);title('nonlinear function 1 compensation','FontSize',15);
   subplot(az(4));
   plot(tu(2,:),x(2,:),'*','MarkerSize',2);title('nonlinear function 2 compensation','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);