clear all;
close all;
clc;

%déclaration de variables
Ts = 4; %sous-entendu 4*Te
K = 1000; %nombre de symboles;
alpha = 4;
R = 4;
N = 10;

%génération des data
%a = 2*randi([0,1],1,K)-1;
a=2*rand(1,K);
a=-(a<0.5)+(a>=0.5);
s = kron(a,[1,zeros(1,Ts-1)]);
g1 = [-0.049 0.482 -0.556 1 -0.171];
g2 = [0.443 1 0.921 0.189 -0.087];
g3 = [-0.211 -0.199 1 -0.284 0.136];
g4 = [0.417 1 0.873 0.285 -0.049];
g = [g1,g2,g3,g4];
for n = 1 : K - R - 1
    r1(n) = fliplr(a(n:R+n))*g1.';
    r2(n) = fliplr(a(n:R+n))*g2.';
    r3(n) = fliplr(a(n:R+n))*g3.';
    r4(n) = fliplr(a(n:R+n))*g4.';
end
r1 = kron(r1,[1,0,0,0]);
r2 = kron(r2,[0,1,0,0]);
r3 = kron(r3,[0,0,1,0]);
r4 = kron(r4,[0,0,0,1]);
r = r1+r2+r3+r4;

%implémentation de la méthode
rr = [];
for n = 1 : K - R - 1 - N   
     rn1 = fliplr(r(1+(n-1)*alpha:alpha:(n-1)*alpha+alpha*N)).';
     rn2 = fliplr(r(2+(n-1)*alpha:alpha:1 + (n-1)*alpha+alpha*N)).';
     rn3 = fliplr(r(3+(n-1)*alpha:alpha:2 + (n-1)*alpha+alpha*N)).';
     rn4 = fliplr(r(4+(n-1)*alpha:alpha:3 + (n-1)*alpha+alpha*N)).';
     rn = [rn1;rn2;rn3;rn4];
     rr = [rr,rn];
end 
Rrn = rr*rr'/(K-N); %estimation de la matrice de covariance du vecteur aléatoire rn
[V0,D0] = eig(Rrn);

%Construction de la matrice VR (26 est le nombre de VP nulles)
for k=1:25
    for n=1:R+1  
        VR1(n,:)=[zeros(1,n-1),V0(1:N,k).', zeros(1,R+1-n)];       %de taille R+1 x R+N c-à-d 5 x14  
        VR2(n,:)=[zeros(1,n-1),V0(N*1+1:2*N,k).', zeros(1,R+1-n)]; %de taille R+1 x R+N c-à-d 5 x14
        VR3(n,:)=[zeros(1,n-1),V0(N*2+1:3*N,k).', zeros(1,R+1-n)]; %de taille R+1 x R+N c-à-d 5 x14
        VR4(n,:)=[zeros(1,n-1),V0(N*3+1:4*N,k).', zeros(1,R+1-n)]; %de taille R+1 x R+N c-à-d 5 x14
    end   
        %Construction de la matrice DELTA 
        VRR=[VR1; VR2; VR3; VR4];    %de taille alpha(R+1) x (R+N) c-à-d 20 x14
        Delta_ini(:,:,k) =VRR*VRR';  %de taille 20 x 20
end

%Construction de la matrice DELTA après la somme (equation 17)
Delta=sum(Delta_ini,3);

%Identification du vecteur g
[V1,D1] = eig(Delta);
gest = V1(:,1).';
erreur = abs(g*gest')/(norm(g)*norm(gest))

%représentation graphique
figure(1);
subplot(2,1,1)
plot(s(1:100));
subplot(2,1,2)
plot(r(1:100));
