pearson_ica_demo.m

PearsonICA程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is a demonstration of the Pearson-ICA algorithm in work.
%
% Copyright (c) Helsinki University of Technology,
% Signal Processing Laboratory,
% Juha Karvanen, Jan Eriksson, and Visa Koivunen.
%
% For details see the files readme.txt
% and gpl.txt provided within the same package.
%
% Last modified: 25.9.2000
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf(['\nThis is a demonstration of the Pearson-ICA algorithm in work.\n']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Source signals
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % The number of samples (signal length).
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  signal_length=5000;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % The number of signals.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  numsig=4; 

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % The Rayleigh distributed source signals are generated.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  sources=sqrt((randn(numsig,signal_length).^2)+(randn(numsig, ...
	 					  signal_length).^2));

  fprintf(['\nFirst, %d zero mean unit variance signals of length %d '...
           'will be\n'],numsig, signal_length);
  fprintf(['       generated from the Rayleigh' ...
	   ' distribution.\n']);
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % To illustrate the quality of separation a short sequence of signal means 
  % is substituted to the begining of every source signal. In the plot these 
  % sequences are seen as constant valued parts of the signals. Because the 
  % added sequences do not overlap, they are not seen in the mixed signals.
  % When the separation is successful the same constant sequences are
  % again visible in the separation result.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % The number of observations from the beginning to be plotted.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    lastplot=round(signal_length/4000*numsig)*20;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % The length of mean vector to be substituted.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    lag_length=round(lastplot/numsig);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % The substitution of means.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for signal_no=1:numsig,
      zero_beg(signal_no)=(signal_no-1)*lag_length+1;
      zero_end(signal_no)=signal_no*lag_length;
      source_mean=mean(sources(signal_no,...
	       [1:zero_beg(signal_no)-1 zero_end(signal_no)+1:end]));
      sources(signal_no,zero_beg(signal_no):zero_end(signal_no))=...
	       source_mean;
    end

    fprintf(['\n       A sequence of zeros (signal means) of length %d '...
	   'is substituted to\n'],lag_length); 
    fprintf(['       every source signal such that the indices '...
	   'do not overlap. This is\n']);
    fprintf(['       done in order to visualize the separation' ...
	   ' result. If the final\n']);
    fprintf(['       separation is succesful, the same constant '... 
	   'sequences should be again\n']);
    fprintf(['       visible.\n']); 
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Ploting the beginning of the source signals.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  H1=figure; 
  for signal_no=1:numsig,
    subplot(numsig,1,signal_no);
    hold on;
    plot(sources(signal_no,1:lastplot),'b');
    interval=zero_beg(signal_no):zero_end(signal_no)-1;
    plot(interval,sources(signal_no,interval),'r');  
    xlabel('sample no.'); 
    V=axis;
    for signal_no2=1:numsig-1,
      line(zero_end(signal_no2),V(3):0.2:V(4));
    end
    if signal_no==1, title('The beginning of the signals'); end;
  end
  fprintf(['\n       The first %d samples of the source signals are '...
	   'shown in Figure No. %d.\n'],lastplot,H1);
  fprintf(['       The substituted means are drawn with red.\n']);
  fprintf('\nPress any key to continue...\n');
  pause;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mixing matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  MIX_MATRIX=randn(numsig);
  fixed_matrix=[0.7396    0.9084    0.2994    0.3089; ...
                0.4898    0.2980    0.5771    0.4108; ...
                0.1096    0.7808    0.8361    0.4669; ...
                0.4199    0.8799    0.2706    0.7467];
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % To use completely random matrix, comment out the next two lines.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  MIX_MATRIX(1:min(numsig,4),1:min(numsig,4))=...
                                 fixed_matrix(1:min(numsig,4),1:min(numsig,4));
 
  fprintf(['\nSecond, the source signals are linearily mixed '...
	   'together using the matrix\n']);
  fprintf('       [ ');
  fprintf('%+1.4f ',MIX_MATRIX(1,:));
  fprintf(';\n');
  for mix_no=2:numsig-1,
    fprintf('         ');
    fprintf('%+1.4f ',MIX_MATRIX(mix_no,:));
    fprintf(';\n');
    end
  fprintf('         '); 
  fprintf('%+1.4f ',MIX_MATRIX(end,:));
  fprintf(']\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mixed signals (observation matrix)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  mixedsig=MIX_MATRIX*sources;
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Ploting the beginning of the mixed signals.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  H2=figure;
  for signal_no=1:numsig,
    subplot(numsig,1,signal_no);
    hold on;
    plot(mixedsig(signal_no,1:lastplot),'b');
    xlabel('sample no.');
    V=axis;
    for signal_no2=1:numsig-1,
      line(zero_end(signal_no2),V(3):0.2:V(4));
    end
    if signal_no==1, title('The beginning of the observed mixtures'); end;
  end
  fprintf(['\n       The first %d samples of the mixtures are '...
	   'shown in Figure No. %d.\n'],lastplot,H2);
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Separation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  fprintf('\nNow we will try to recover original signals from the\n');
  fprintf('       mixture only using the Pearson-ICA algorithm.\n');
  fprintf('\nPress any key to begin...\n');
  pause;
  fprintf('\n');
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
  % The maximum number of iterations is set to 100. Otherwise
  % the standard parameter values are used. 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  [icasig, A, W]= pearson_ica(mixedsig,'maxNumIterations',100);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Ploting the beginning of the separated signals.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  H3=figure;
  for signal_no=1:numsig,
    subplot(numsig,1,signal_no);
    hold on;
    plot(icasig(signal_no,1:lastplot),'b');
    xlabel('sample no.');
    V=axis;
    for signal_no2=1:numsig-1,
      line(zero_end(signal_no2),V(3):0.2:V(4));
    end
    if signal_no==1, title('The beginning of the separated signals'); end;
  end  
  fprintf(['\nThe first %d samples of the separated signals are '...
	   'shown in Figure No. %d.\n'],lastplot,H3);
  fprintf(['Compare the result to the originals shown in '...
	   'Figure No. %d. '],H1);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
  % In order to compare directly the results, the sign and the permutation
  % indeterminency must be resolved. This is done by selecting the largest
  % absolut row values (and the corresponding signs) from 
  % the "result matrix" W*MIX_MATRIX. 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  Y=W*MIX_MATRIX;
  [value,permutation]=max(abs(Y),[],2);
  source_sign=sign(diag(Y(:,permutation)));

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
  % Also the signals have to have the same scale. This is done
  % by normalizing the sources and the separation results.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  normalized_sources=sources-(ones(signal_length,1)*mean(sources,2)')';
  normalized_sources=normalized_sources./(ones(signal_length,1)...
		       *sqrt(mean(normalized_sources.^2,2))')';
  normalized_icasig=icasig-(ones(signal_length,1)*mean(icasig,2)')';
  normalized_icasig=normalized_icasig./(ones(signal_length,1)...
		       *sqrt(mean(normalized_icasig.^2,2))')';
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Resolving permutation and sign. 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  normalized_icasig=normalized_icasig.*(ones(signal_length,1)*source_sign')';
  normalized_icasig(permutation,:)=normalized_icasig;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Ploting the beginnings of the normalized signals. 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  H4=figure;
  for signal_no=1:numsig,
    subplot(numsig,1,signal_no);
    hold on;
    plot(normalized_sources(signal_no,1:lastplot),'r'); 
    plot(normalized_icasig(signal_no,1:lastplot),'b');
    xlabel('sample no.');
    V=axis;
    for signal_no2=1:numsig-1,
      line(zero_end(signal_no2),V(3):0.2:V(4));
    end
    if signal_no==1, title(['The beginnings of the normalized source '...
		    '(red) and separated signals (blue)']); 
    end;
  end 
  fprintf('The first %d\n',lastplot);
  fprintf(['samples of the sign and permutation resolved normalized '...
	  '(mean=0, variance=1)\n']);
  fprintf(['separation results are plotted (blue) on the top of '...
	   'the normalized\n']);
  fprintf(['source signals (red) in Figure No. %d.\n'],H4);
  
	   
% The end of the file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%