%Very basic Matlab code for the unifying model for BSS
%as proposed in my article in Signal Processing, 85(7):1419-1427, 2005
%Aapo Hyvarinen, May 2003 -- Mar 2006

randn('state',1)  %fix random number generator for published results

clear error  %this stores the errors for artificial data

%Main structural parameters
N=5000; %data size
n=6; %number of mixtures and sources 


%THE WHOLE STUFF WILL BE REPEATED 100 TIMES TO EVALUATE THE ALGORITHM

for iter=1:100

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CREATE ARTIFICIAL DATA FOR TESTING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%create source signals 
%two nongaussian and two gaussian time-correlated signals
%using only two different time-correlations
%and two nonstatvar signals
%thus only unibss criteria estimates all!


%autoregressive constants

alpha1=.33;  
alpha2=.75;  

iteraatio=iter %show this to user


S=[];
S(:,1)=randn(n,1);
for t=2:N;
  S(1,t)=alpha1*S(1,t-1)+abs(randn).^1.3.*sign(randn);
  S(2,t)=alpha2*S(2,t-1)+abs(randn).^1.3.*sign(randn);
  S(3,t)=alpha1*S(3,t-1)+randn;
  S(4,t)=alpha2*S(4,t-1)+randn;
  S(5,t)=.9*S(5,t-1)+randn;
  S(6,t)=.9*S(6,t-1)+randn;
end
%make nonstatvar signals
S(5:6,:)=S(5:6,:).*sign(randn(2,N));

%normalize

S=S-mean(S')'*ones(1,N);

for t=1:n
S(t,:)=S(t,:)/std(S(t,:));
end


%create mixtures

A=randn(n,n); %mixing matrix
X=A*S; %observed data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  THE UNIBSS ALGORITHM           %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%The input to this part is the data matrix X
%and the structural variables n and N (which are simply given by size(X))

%1. PREPROCESSING OF THE DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%

%(it is assumed that the data is already zero-mean)

%whitening

C=cov(X');
V=inv(sqrtm(C));
Z=V*X;

%compute time-delayed data

Zd=[Z(:,N),Z(:,1:(N-1))]; %we could also zero-pad
CZd=Z*Zd'/N; %time-delayed covariance


%2. INITIALIZE ALGORITHM
%%%%%%%%%%%%%%%%%%%%%


W=randn(n,n);
W=inv(sqrtm(W*W'))*W;



%3. THE MAIN ALGORITHM PART:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%this just takes an certain number of steps, 
%I didn't bother to implement a proper stopping criterion.

	    %Also in other respects, the code is extremely basic:
%- a single nonlinearity (-sign(.)) is implemented; this assumes supergaussian innovations
%- no effort is made to estimate the filter "h"
%- only a first-order autoregressive model is used

for t=1:999 


	for i=1:n %for each source do separately
          %compute AR(1) coefficient
  	  alpha=W(i,:)*CZd*W(i,:)';
          %compute estimate of innovations process
	  n_est=W(i,:)*(Z-alpha*Zd);
          %compute local standard deviations
	  sigma=sqrt(filter([1,1,1,1,1,1],1,n_est.^2));
          %compute gradient
          grad(i,:)=((Z-alpha*Zd)*(-sign(n_est)./sigma)')'/N;
	end
        
        %project gradient and do gradient step
        W=W+.1*(grad-W*grad'*W);

        %orthogonalize W
        W=inv(sqrtm(W*W'))*W;

	%evaluate results
	comb=W*V*A;
        %a very simple measure of estimation error, not unlike Amari error
        errornow=(n-sum(max(comb.^2)))^2;
	error(iter,t)=errornow;


end %of algorithm iteration "for t"

end %of the outer iteration that repeats estimation 100 times


%At this point, you can plot the rows of "error", or a median of them,
%to get the figure in the paper

%Here's another simple evaluation of results:
number_of_iterations_not_converged=sum(error(:,999)>.01)