www.pudn.com > twosignals.rar > state.m

% 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) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%