% test de l'ICA sur 3 sources MSK
% ecrit a partir de main.m

% 19 aout 98
% F: identifie le melange, cad A = F*Lambda*P idealement
% S: estime les sources

%%%%%%%%%% A EXECUTER POUR DES DONNEES EXPLICITES %%%%%%%%%%%%%%%%%
clear;
randn('state',0);
rand('state',0);

ii=sqrt(-1);tol=1e-9;format compact;
path(path,'/home/grellier/Matlab/Fonctions');
dur=1000;
fprintf('INTEGRATION: dur = %g \n',dur);

K1 = 2;
K2 = 5;
P = 3;
P2 = 4;
NbTests = 50;
NbR = 11;
 
% Chargement du best Mix
load BPSKK2P3.mat -ascii;
RSB = BPSKK2P3(:,1);
A1  = BPSKK2P3(:,2);
A2  = BPSKK2P3(:,3);
A3  = BPSKK2P3(:,4);
A4  = BPSKK2P3(:,5);
NbRSB = length(RSB);
Ber1 = zeros(1,NbTests);
Ber2 = zeros(1,NbTests);

% Calcul du melange
A = [A1(NbR), i*A2(NbR), 0, 0;
       0, i*A3(NbR), A4(NbR),0];

for tests = 1:NbTests,
dur = 1000;
%%% SOURCES
x(1:P,1:dur)=zeros(P,dur);
for i=1:P,
  x(i,1:dur)=msk2(dur);ect=std(x(i,1:dur));x(i,1:dur)=x(i,1:dur)/ect;
end;
bleu = (-1).^[1:1:dur];
x(4,1:dur) = bleu.*x(1,:).*x(2,:).*x(3,:);

w=-1%          %%% signe presume des kurtosis source

%%% OBSERVATIONS
brt = bruit(K1,dur,RSB(NbR));
yn(1:K1,1:dur)=A*x+brt;
n1 = 1:2:dur; n2 = 2:2:dur;
y1 = yn(1,n1);y2 = yn(2,n1);
y(1:K2,1:dur/2)=[y1;y2;y1.^3;(y1.^2).*y2;y2.^3];
dur = dur/2;

%%% OBSERVATIONS STANDARDISEES
[U,S,V]=svd(y',0);tol=1e-9;
s=diag(S);I=find(s<tol);r=K2;
if length(I)~=0, r=I(1)-1;U=U(:,1:r);S=S(1:r,1:r);V=V(:,1:r);end;
Z=U'*sqrt(dur);L=V*S'/sqrt(dur);      %%% on a y=L*Z;
F=L;                            %%%%%% INITIALISATION DE LA TRANSFORMATION

%%% IDENTIFICATION DU MELANGE PAR BALAYAGES
%%% F identifie le melange, les matrices Q le separent iterativement
%%%%%% CONTRASTE INITIAL
T=dur;contraste=0;
for i=1:r,
 gii=Z(i,:)*Z(i,:)'/T;Z2i=Z(i,:).^2;;giiii=Z2i*Z2i'/T;
 qiiii=giiii/gii/gii-3;contraste=contraste+w*qiiii;
end;
contI=contraste;
contM=0;for i=1:P,
 gii=x(i,:)*x(i,:)'/T;x2i=x(i,:).^2;;giiii=x2i*x2i'/T;
 qiiii=giiii/gii/gii-3;contM=contM+w*qiiii;
end;

%%%% TRANSFORMATIONS UNITAIRES
S=Z;
if r==2,K=1;else,K=1+round(sqrt(r));end;  %%% K= NBRE MAXI DE BALAYAGES
Rot=eye(r);K=K+2;
for k=1:K,                         %%%%%% DEBUT BALAYAGES
Q=eye(r);
  for i=1:r-1,
  for j= i+1:r,
    S1ij=[S(i,:);S(j,:)];
    [qij,theta]=tetan_nv(S1ij,w);  %%%%%% ROTATION PLANE SEPARANTE
    Sij=qij*S1ij;                  %%%%%% TRAITEMENT D'UNE PAIRE
    S(i,:)=Sij(1,:);S(j,:)=Sij(2,:);
    Qij=eye(r);Qij(i,i)=qij(1,1);Qij(i,j)=qij(1,2);
    Qij(j,i)=qij(2,1);Qij(j,j)=qij(2,2);
    Q=Qij*Q;                       %%%%%% CUMUL DE LA ROTATION SUR LE BALAY
  end;
  end;
Rot=Rot*Q';                        %%%%%% CUMUL DE TOUTES LES ROTATIONS
end;                               %%%%%% FIN BALAYAGES
F=F*Rot;                           %%%%%% CUMUL DE LA TRANSFORMATION

%%%%%% CONTRASTE FINAL
S=Rot'*Z;   %%% ON RECALCULE LES SOURCES ESTIMEES S (facultatif)
T=length(S);contraste=0;
for i=1:r,
 gii=S(i,:)*S(i,:)'/T;S2i=S(i,:).^2;;giiii=S2i*S2i'/T;
 qiiii=giiii/gii/gii-3;contraste=contraste+w*qiiii;
end;

%fprintf(' contraste initial = %g,\n contraste final = %g\n',contI,contraste)
%fprintf(' contraste maximal = %g,\n',contM)

%%% TRACES
%figure(1);clf;Tp=100;
%for i=1:P2,
%subplot(P2,3,(i-1)*3+1);plotc(x(i,1:Tp));
%title(['source ' num2str(i) ' vraie']);
%subplot(P2,3,(i-1)*3+2);plotc(y(i,1:Tp));
%title(['observation' num2str(i) ', SNR = ' num2str(RSB(NbR)) 'dB']);
%subplot(P2,3,(i-1)*3+3);plotc(S(i,1:Tp));axis([-2 2 -2 2]);
%title(['source ' num2str(i) ' estimee']);
%end;

for i=1:4,
  S1c(i,:) = projectcpa(S(i,:),[1 -1]);
end;
for i = 1:4,
  for j = 1:4,
    Ber(i,j) = bpskcalber(x(i,n1),S1c(j,:));
  end;
end;
Ber,
[H,I] = min(Ber);
Ber1(1,tests) = 0;
for i = 1:4,
  Ber1(1,tests) = Ber(I(i));
end;