clear all;
close all;
global N;
global te

tmax=10000;
te=.1;
N=tmax/te;
t=0:te:tmax-te;


%entrée gaussienne 
x=rand(1,N);
%filtrage pour lissage
xx=filter([0.01],[1 -0.99],x);
xxx=filter([0.01],[1 -0.99],xx);
s(1,:)=filter([0.01],[1 -0.99],xxx);   

x=rand(1,N);
%filtrage pour lissage
xx=filter([0.001],[1 -.985],x);
xxx=filter([0.001],[1 -.985],xx);
s(2,:)=filter([0.001],[1 -.985],xxx);   

clear xx
clear xxx




%elimination de la partie initial non significative
Nc=1000;
s=s(:,Nc+1:N);
N=N-Nc;
tmax=N*te;
t=0:te:tmax-te;


%whitening
s_s=[mean(s(1,:).^2) mean(s(1,:).*s(2,:)); mean(s(2,:).*s(1,:)) mean(s(2,:).^2)];


[u,lambda,v]=svd(s_s);
W=lambda^(-.5)*v';

s=W*s;


% figure(1)
% plot(t,s(1,:))
% figure(2)
% plot(t,s(2,:))
% figure(3)
% plot(s(1,:),s(2,:),'.');
% drawnow



%derivative des signaux
der_s=(s(:,2:N)-s(:,1:N-1))/te;


% figure(4)
% plot(0:te:tmax-2*te,der_s(1,:))
% figure(5)
% plot(0:te:tmax-2*te,der_s(2,:))
% figure(6)
% plot(der_s(2,:),der_s(1,:),'.')
 



%Melange linéaire des sources s
A=rand(2,2);
y=A*s;




%calcul de la dérivé des observations

der_y=(y(:,2:N)-y(:,1:N-1))/te;

% figure(18)
% plot(y(1,:),y(2,:),'.')
% figure(19)
% plot(der_y(1,:),der_y(2,:),'.')


%variance covariance of the observation signals and their derivatives
mom_y=[    
    mean(y(1,:).*y(1,:))
    mean(y(2,:).*y(2,:))
    mean(y(1,:).*y(2,:))

    mean(der_y(1,:).*der_y(1,:))
    mean(der_y(2,:).*der_y(2,:))
    mean(der_y(1,:).*der_y(2,:))
]



%analytical solution

alpha=(mom_y(1)*mom_y(6)-mom_y(4)*mom_y(3))/(mom_y(3)*mom_y(5)-mom_y(6)*mom_y(2));
beta=.5*(mom_y(1)*mom_y(5)-mom_y(2)*mom_y(4))/(mom_y(3)*mom_y(5)-mom_y(6)*mom_y(2));



eta(1)=-beta*(1+sqrt( 1- alpha/(beta^2)  ));
eta(2)=-beta*(1-sqrt( 1- alpha/(beta^2)  ));

b(1)=+ inv( sqrt( mom_y(1)+2*eta(1)*mom_y(3)+eta(1)^2*mom_y(2) ) );
b(2)= b(1)*eta(1);
b(3)=+ inv( sqrt( mom_y(1)+2*eta(2)*mom_y(3)+eta(2)^2*mom_y(2) ) );
b(4)= b(3)*eta(2);

B=[b(3) b(4); b(1) b(2) ];


s_estim=B*y;



figure(7)
subplot(2,2,1)
plot(t,s_estim(1,:),'b');
hold on
plot(t,s(1,:),'r');
subplot(2,2,2)
plot(t,s_estim(2,:),'b');
hold on
plot(t,s(2,:),'r');
subplot(2,2,3)
plot(t,s_estim(1,:),'b');
hold on
plot(t,s(2,:),'r');
subplot(2,2,4)
plot(t,s_estim(2,:),'b');
hold on
plot(t,s(1,:),'r');



