state.m

Basic projection pursuit algorithm demonstrated on 2 sound signals。 include: atutorialintroduction

% Set rand number seeds.
seed=1;randn('state',seed);rand('state',seed);

num_sources		= 3; 
num_mixtures	= num_sources;
num_samples 	= 5000;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GET DATA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Define max mask len for convolution.
max_mask_len= 500;
% [8] n = num half lives to be used to make mask.
n			= 8; 

% Make source signals as set of increasingly smooth signals.

% Make mask.
% h= half-life of exponential in mask which is then convolved with random signal.
h=2; t = n*h; lambda = 2^(-1/h); temp = [0:t-1]'; lambdas = ones(t,1)*lambda; mask = lambda.^temp; 
mask1 = mask/sum(abs(mask));
h=4; t = n*h; lambda = 2^(-1/h); temp = [0:t-1]'; lambdas = ones(t,1)*lambda; mask = lambda.^temp; 
mask2 = mask/sum(abs(mask));
h=8; t = n*h; lambda = 2^(-1/h); temp = [0:t-1]'; lambdas = ones(t,1)*lambda; mask = lambda.^temp; 
mask3 = mask/sum(abs(mask));

sources = randn(num_samples,num_sources);
sources(:,1)=filter(mask1,1,sources(:,1));
sources(:,2)=filter(mask2,1,sources(:,2));
sources(:,3)=filter(mask3,1,sources(:,3));

% Make mixing matrix.
A = randn(num_sources,num_sources);

% Make mixtures.
mixtures = sources*A;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% COMPUTE V AND U.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set short and long half-lives.
shf 		= 1; 
lhf 		= 900000; 	

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get masks to be used to find (x_tilde-x) and (x_bar-x)
% Set mask to have -1 as first element, and remaining elements must sum to unity.

% Short-term mask.
h=shf; t = n*h; lambda = 2^(-1/h); temp = [0:t-1]'; 
lambdas = ones(t,1)*lambda; mask = lambda.^temp;
mask(1) = 0; mask = mask/sum(abs(mask));  mask(1) = -1;
s_mask=mask; s_mask_len = length(s_mask);

% Long-term mask.
h=lhf;t = n*h; t = min(t,max_mask_len); t=max(t,1);
lambda = 2^(-1/h); temp = [0:t-1]'; 
lambdas = ones(t,1)*lambda; mask = lambda.^temp;
mask(1) = 0; mask = mask/sum(abs(mask));  mask(1) = -1;
l_mask=mask; l_mask_len = length(l_mask);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Filter each column of mixtures array.
S=filter(s_mask,1,mixtures); 	
L=filter(l_mask,1,mixtures);

% Find short-term and long-term covariance matrices.
U=cov(S,1);		
V=cov(L,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIND **SINGLE** SOURCE SIGNAL USING GRADIENT ASCENT.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cl=V;
cs=U;

% Make initial weight vector.
w=randn(1,num_sources);
w=w/norm(w);
w0=w;

% Use initial w0 to extract source
y0=mixtures*w0';

% Find correlation of y0 with sources
rs0=corrcoef([y0,sources]);
fprintf('Using Grad ascent: Correlation of single with sources extracted by initial w ...\n');
abs(rs0(1,2:4))

% Set learning rate ...
eta=1e-1;
% Set max number of iterations ...
maxiter=100;

% Make arrays to store results.
gs=zeros(maxiter,1); % gradient magnitude |g|
Fs=zeros(maxiter,1); % function value F

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do gradient ascent ...
for i=1:maxiter

	% Get value of function F
	Vi = w*cl*w';
	Ui = w*cs*w';
	F = log(Vi/Ui);

	% Get gradient
	g = 2*w*V./Vi - 2*w*U./Ui;

	% Update w
	w = w + eta*g;

	% Record results ...
	Fs(i)=F;
	gs(i)=norm(g);

end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Plot results ...
figure(1); plot(Fs); xlabel('Iteration Number'); ylabel('F=log(V/U)');
figure(2); plot(gs); xlabel('Iteration Number'); ylabel('Gradient Magnitude');

% Use w to extract source
y1=mixtures*w';

% Find correlation of y1 with sources
rs=corrcoef([y1,sources]);
fprintf('Using Grad ascent: Correlation of single with sources extracted by initial w ...\n');
abs(rs(1,2:4))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NOW USE W MATRIX FROM EIG FUNCTION TO EXTRACT **ALL** SOURCES.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find optimal solution as eigenvectors W.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[W d]=eig(V,U); W=real(W);

ys = mixtures*W;
a=[sources ys]; c=corrcoef(a);  
rs=c(1:num_sources,num_sources+1:num_sources*2);
fprintf('Using EIG: Correlations between sources and all recovered signals ...\n'); 
abs(rs)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%