clear;
F0 = 1180;
F1 = 980;
Fe = 48000;
Te = 1/Fe;
Ts = 1/300;
Ns = floor(Ts/Te);



% Question 3.1.1.1
% on charge les bits à transmettre, ceux-ci son placés dans la variable "bits"
load('DonneesBinome1.mat');
% on répète chaque bits Ns fois pour obtenir NRZ
NRZ = repelem(bits, Ns);
% on génère notre échelle temporelle
T = 0:Te:Ns*length(bits) /Fe - Te;

% Question 3.1.1.2
%plot(T, NRZ);
%xlabel("temps (s)");
%ylabel("NRZ(t)");

% Question 3.1.1.3
%N_fft = 48000;
%F = 0:1/N_fft:Fe;
%perio = abs(fft(NRZ, N_fft)).^2/N_fft;
%semilogy(perio);
%periodogram(NRZ);



% Question 3.1.2.1
% on génère deux phases aléatoires pour les cosinus de x(t)
phi = 2*pi*rand(2, 1);
% on peut finalement créer X via l'équation (1)
X = (1-NRZ).*cos(2*pi*F0*T + phi(1)) + NRZ.*cos(2*pi*F1*T + phi(2));

%plot(T, X);
%xlabel("temps (s)");
%ylabel("X(t)");

% Question 3.1.2.4
%periodogram(X);



% Question 3.2
SNR = 10;
Px = mean(abs(X).^2);
Pb = Px / 10^(SNR/10); % d'après l'équation (2)
bruit = sqrt(Pb)*randn(1, length(X));
Xb = X + bruit;

%periodogram(Xb);
%plot(Xb);
%semilogy([0:N-1]/N*Fe, abs(fft(Xb, N)));



% on choisit 4000 Hz en fréquence de coupure car au milieu de F0 et F1
fc_tild = 1080/Fe;
% on choisit un ordre élevé pour être sûr que les fréquences indésirables soient filtrées
ordre = 100;
% on utilise ici un filtre RIF
A = 1;
% on génère B grâce à la transformé inverse d'une fonction porte
B = 2*fc_tild*sinc(2*fc_tild*[-ordre/2:1:ordre/2]);
% on rajoute des zeros après Xb pour contrer le retard créé par le filtrage
Xb = [Xb zeros(1, ordre/2)];



% Question 3.3.1 : Passe-bas
% on utilise la fonction filter de matlab
Xf_bas = filter(B, A, Xb);
% on supprime les premières valeurs pour ne plus avoir de retard.
Xf_bas = Xf_bas(ordre/2+1:length(Xb));

%plot(Xf);

%plot([0:N-1]/N*Fe, abs(fft(B, N)));
%legend("fft(B)");
%title("Filtre passe bas d'ordre 100");
%xlabel("Fréquence (Hz)");
%ylabel("Gain (dB)");

%plot([-ordre/2:1:ordre/2], B);
%plot(B)
%title("Coefficients de B pour un passe bas d'ordre 100");
%ylabel("B");

%semilogy([0:N-1]/N*Fe, abs(fft(Xb, N)));
%hold;
%semilogy([0:N-1]/N*Fe, abs(fft(Xf, N)));
%legend("Xb", "Xf\_bas");
%title("fft du signal avant et après filtrage");
%xlabel("Fréquence (Hz)");
%ylabel("Amplitude (dB)");

%periodogram(Xb);
%hold;
%periodogram(Xf);
%legend("Xb", "Xf\_bas");
%h = get(gca, 'Children');
%set(h(1), 'Color', [0.8500, 0.3250, 0.0980]);
%title("périodogramme avant et après filtrage");



% Question 3.3.2 : Passe-haut
% on modifie B pour passer d'un passe-bas à un passe-haut
B(ordre/2+1) = 1 - 2*fc_tild;
% on utilise la fonction filter de matlab
Xf_haut = filter(B, A, Xb);
% on supprime les premières valeurs pour ne plus avoir de retard.
Xf_haut = Xf_haut(ordre/2+1:length(Xb));

%plot(Xf);

%plot([0:N-1]/N*Fe, abs(fft(B, N)));
%legend("fft(B)");
%title("Filtre passe haut d'ordre 100");
%xlabel("Fréquence (Hz)");
%ylabel("Gain (dB)");

%plot([-ordre/2:1:ordre/2], B);
%plot(B)
%title("Coefficients de B pour un passe haut d'ordre 100");
%ylabel("B");

%semilogy([0:N-1]/N*Fe, abs(fft(Xb, N)));
%hold;
%semilogy([0:N-1]/N*Fe, abs(fft(Xf, N)));
%legend("Xb", "Xf\_haut");
%title("fft du signal avant et après filtrage");
%xlabel("Fréquence (Hz)");
%ylabel("Amplitude (dB)");

%periodogram(Xb);
%hold;
%periodogram(Xf);
%legend("Xb", "Xf\_haut");
%h = get(gca, 'Children');
%set(h(1), 'Color', [0.8500, 0.3250, 0.0980]);
%title("périodogramme avant et après filtrage");



% Question 3.3.5
K = 30; % déterminé expérimentalement
% on calcule l'énergie par groupes de taille Ns
somme = sum((reshape(Xf_bas, [Ns, length(bits)])).^2);
% on filtre les valeurs supérieures à K
result = (somme > K);
% on calcule le taux d'erreur avec les bits originaux
erreur = sum((bits - result).^2)/length(bits)



% Question 4.1
% On génère les deux branches différentes
branche_haut = sum(reshape(X.*cos(2*pi*F0*T + phi(1)), [Ns, length(bits)]));
branche_bas  = sum(reshape(X.*cos(2*pi*F1*T + phi(2)), [Ns, length(bits)]));
% on somme les branches
somme = branche_bas - branche_haut;
% on filtre les valeurs positives et négatives
result = (somme > 0);
% on calcule le taux d'erreur avec les bits originaux
erreur = sum((bits - result).^2)/length(bits)



% Question 4.2
% On génère les quatre branches différentes
branche_1 = sum(reshape(X.*cos(2*pi*F0*T), [Ns, length(bits)])).^2;
branche_2 = sum(reshape(X.*sin(2*pi*F0*T), [Ns, length(bits)])).^2;
branche_3 = sum(reshape(X.*cos(2*pi*F1*T), [Ns, length(bits)])).^2;
branche_4 = sum(reshape(X.*sin(2*pi*F1*T), [Ns, length(bits)])).^2;
% on somme les branches
somme = branche_4 + branche_3 - branche_2 - branche_1;
% on filtre les valeurs positives et négatives
result = (somme > 0);
% on calcule le taux d'erreur avec les bits originaux
erreur = sum((bits - result).^2)/length(bits)



% on affiche l'image à partir des bits récupérés 
%pcode reconstitution_image;
%reconstitution_image(result);
%which reconstitution_image;