%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Implémentation symmétrique de l'algorithme FastICA         
% Version du 6/10/2015
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Définition des constantes -> Pour l'instant, il est nécessaire de
% modifier nData et nIC selon la taille des données et le nombre de 
% composantes indépendantes à extraire respectivement. De plus, les 
% signaux à traiter, si obtenus de l'extérieur, doivent être placés
% dans une matrice de taille nIC x nData
clc;
close all;
clear all;
nData = 1000; % À adapter selon la taille des données
freq1 = 7;
freq2 = 11;
omega1 = 2*pi*freq1;
omega2 = 2*pi*freq2;
nIC = 3;
epsilon_decorr = 0.99999;
epsilon_poids = 0.00001;
nIter = 4;

% Génération de signaux de test

% Onde triangulaire et en dent de scie 
T = 10*(1/50);
Fs = 5*nData;
dt = 1/Fs;
t = 0:dt:T-dt;
sig2 = sawtooth(2*pi*50*t, 0.5);
sig3 = sawtooth(2*pi*20*t);

% Onde sinusoïdale
n = 1;
for n=1:nData
    
    sig1(1,n) = cos(n*(omega1/nData));
    %sig2(1,n) = sin(n*(omega2/nData));
    
end
n = 1;


% Restructuration des signaux de test
sigs = [sig1 ; sig2 ; sig3];

%Ajout de bruit
sigs = awgn(sigs,100);

% Matrice de mélange
Aorig=rand(size(sigs,1));
%Aorig = [0.1547 0.2569 0.3569 ; 0.8945  0.9214  0.9852 ; 0.1365  0.3214 0.1965];

% Pour des signaux réels, les placer dans mixedsig et enlever
% "(Aorig*sigs)" et utiliser load('filepath')
mixedsig=(Aorig*sigs);

%Soustraction de la moyenne 
for n=1:nIC
mixedsig(n,:) = mixedsig(n,:) - mean(mixedsig(n,:));
%mixedsig(2,:) = mixedsig(2,:) - mean(mixedsig(2,:));
%mixedsig(3,:) = mixedsig(3,:) - mean(mixedsig(3,:));
end

% Calcul de la matrice de covariance
covar = cov(mixedsig');

% "Whitening" pour avoir une variance unitaire
[V,D] = eig(covar);

for n = 1:nIC
D(n,n) = 1/sqrt(D(n,n));
%D(2,2) = 1/sqrt(D(2,2));
%D(3,3) = 1/sqrt(D(3,3));
end


%signals = [mixedsig(1,:) ; mixedsig(2,:) ; mixedsig(3,:)];

new_sigs = V*D*V'*mixedsig;
signals = new_sigs;

%new_sig1 = new_sigs(1,:);
%new_sig2 = new_sigs(2,:);

% Matrice de vecteurs de poids de l'itération précédente
prev_w = zeros(n,n);%[0 0 0 ; 0 0 0 ; 0 0 0];
%w_plus = [0 0 0 ; 0 0 0 ; 0 0 0];

%Initialisation de la matrice de vecteurs de poids
w = orth(rand(nIC));

%w = V*D*V'*w;


%w(1,:) = w(1,:)/norm(w(1,:));
%w(2,:) = w(2,:)/norm(w(2,:));

%prev_dot_prod = zeros(nIC,1); %[0 ; 0 ; 0];
%dot_prod = zeros(nIC,1); %[0 ; 0 ; 0];

k = 1;
%i = 1;
j = 0;
%fin = 0;

%  while j == 0
%     w = w/norm(w*w');    
%     w = 1.5*w - 0.5*((w*w')*w);
%     W = w*w';
%     if (W(1,1) >= epsilon_decorr) && (W(2,2) >= epsilon_decorr)
%         j = 1;
%     end
%  end
%     j = 0;

% Boucle principale pour déterminer la matrice de vecteurs de poids
while k <= nIter 

    w = w*real(sqrtm(inv(w'*w)));
    dot_prod = min(abs(diag(w'*prev_w)));
    
    if (1 - dot_prod <= epsilon_poids)
        %w = ((V')^-1*D^-1*V^-1)*w;
        %W = w' * V*D*V';
        W = w';
        break;
    end
    prev_w = w;
    
    U = new_sigs'*w;
    
    w = (new_sigs*tanh(U))/nData -  (ones(size(w,1),1) * sum(ones(size(U)) .* (1 - tanh(U).^2))).*w/nData;
  
    
%     Implémentation de l'approche "deflation"; non fonctionnelle actuellement
%     
%     for i = 1:nIC
%     w_plus(i,:) = mean(new_sigs*tanh(w(i,:)*new_sigs)') - (mean(1 - (tanh(w(i,:)*new_sigs)).^2))*w(i,:)';
%     %w(i,:) = w_plus(i,:)/norm(w_plus(i,:))
%     %prev_dot_prod = dot_prod;
%     
%     end
%     while j == 0
%     w = w/norm(w*w');    
%     w = 1.5*w - 0.5*((w*w')*w);
%     W = w*w';
%     if (W(1,1) >= epsilon_decorr) && (W(2,2) >= epsilon_decorr)
%         j = 1;
%     end
%     end
%     j = 0;
%     dot_prod(1,1) = dot(prev_w(1,:), w(1,:));
%     dot_prod(2,1) = dot(prev_w(2,:), w(2,:));
%     prev_w = w;
%     
%     if (1 - dot_prod(1,:) <= epsilon_poids) && (1 - dot_prod(2,:) <= epsilon_poids)
%         fin = 1;
%     end
%     
%     k = k + 1;
%     new_sigs = w*new_sigs;
    
end
%w = ((V')^-1*D^-1*V^-1)*w;
%res = w'*signals;
%result = res(1,:);

% Extraction des composantes indépendantes
X = W*signals;

%X = w*new_sigs;


% Affichage des composantes indépendantes
for n=1:nIC
figure(n)
plot(X(n,:))
end

% Affichage des signaux à partir desquels les composantes indépendantes
% ont été extraites
for n=1:nIC
figure(n+nIC)
plot(mixedsig(n,:))
end

