www.pudn.com > zxb.rar > egld_ica_demo.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is a demonstration of the EGLD-ICA algorithm in work.
%
% Copyright (c) Helsinki University of Technology,
% Signal Processing Laboratory,
% Jan Eriksson, Juha Karvanen, 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 EGLD-ICA algorithm in work.\n']);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
egld_use;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Source signals
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The number of samples (signal length).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
signal_length=5000;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The moments of GLD distributions from which the sources are
% generated. The fist one corresponds to the normal distribution.
% Your can add or remove signals in this demo by simply changing
% the moment vectors below. E.g. try to remove the rightmost values.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha3=[0 0.5 0.5 0];
alpha4=[3 3 3.5 3.5];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The number of signals (obtained from alpha3 vector).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
numsig=size(alpha3,2);
fprintf(['\nFirst, %d zero mean unit variance signals of length %d '...
'will be\n'],numsig, signal_length);
fprintf([' generated from the Generalized Lambda' ...
' Distribution (GLD).\n']);
fprintf([' The theoretical 3rd and 4th moments of the '...
'distributions are\n']);
fprintf(' ');
fprintf('(%1.2f,%1.2f), ',[alpha3(1:end-1);alpha4(1:end-1)]);
fprintf('and (%1.2f,%1.2f).\n',alpha3(end),alpha4(end));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The parameters of the GLD distributions are obtained with the
% method of moments.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for signal_no=1:numsig,
lambda(signal_no,:)=gld_momentfit(alpha3(signal_no), alpha4(signal_no));
end
fprintf([' The corresponding lambda parameter '...
'values are given by\n']);
fprintf(' [%+1.4f, %+1.4f, %+1.4f, %+1.4f],\n',lambda(1:end-1,:)');
fprintf(' and [%+1.4f, %+1.4f, %+1.4f, %+1.4f].\n',lambda(end,:));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The GLD distributed source signals are generated.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for signal_no=1:numsig,
sources(signal_no,:)=gld_rnd(lambda(signal_no,:),[1,signal_length]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 EGLD-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] = egld_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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%