clear;
close all;

Fe = 24000; % Hz
Te = 1/Fe; % s
N = 5000; % nombre de bits envoyés
Rb = 6000; % débit binaire

bits = randi([0, 1], 1, N); % bits envoyés

M = 2^1; % signal binaire
Ts = log2(M)/Rb; % période symbole
Ns = floor(Ts/Te);
T = (0:N*Ns/log2(M)-1) * Te; % échelle temporelle

h = ones(1, Ns); % mise en forme: réponse impulsionnelle rectangulaire
h_c = [ 1 zeros(1, Ns-1) ]; % propagation: dirac
h_r = [ones(1, Ns/2) zeros(1, Ns/2)]; % réception: réponse échellon

n0 = 4; % déterminé en traçant g

%% tracé de g
g = conv(conv(h, h_c), h_r); % réponse impulsionnelle globale

figure;
plot(g);
title("Tracé de g");
ylabel("Amplitude");

%% chaine de transmission
X_m = 2 * bits - 1; % mapping binaire à moyenne nulle
X_k = kron(X_m, [1 zeros(1, Ns-1)]); % Suréchantillonnage
X_f = filter(h, 1, X_k); % signal émis
X_c = filter(h_c, 1, X_f); % signal transmis

%% tracé du diagramme de l'oeil du signal non bruité
X_r = filter(h_r, 1, X_c); % signal reçu, sans bruit
figure;
oeil = reshape(X_r, Ns, length(X_r)/Ns); % permet de tracer le diagramme de l'oeil 
plot(oeil);
title("Diagramme de l'oeil, signal non bruité");
xlabel("Temps (s)");
ylabel("Amplitude");

%% tracé de TEB = f(E_b/N_0)
P_x = mean(abs(X_c).^2);
EbN0_db = linspace(0, 8, 100);
EbN0 = 10.^(EbN0_db./10);
TEBs = [];
for e=EbN0
    
    sigma2_x = P_x * Ns / (2 * log2(M) * e); % calcul de sigma^2
    X_b = X_c + sqrt(sigma2_x) * randn(1, length(X_c)); % signal bruité
    X_r = filter(h_r, 1, X_b); % signal reçu
    
    X_e = X_r( n0:Ns:N*Ns ); % échantillonage du signal reçu
    recu = X_e > 0; % on récupère les bits, décision + "démapping"
    TEB = mean((recu - bits).^2); % on calcule l'erreur
    
    TEBs = [ TEBs TEB ];

end

TEB_theorique = qfunc(sqrt(EbN0));
TEB_theorique_ref = qfunc(sqrt(2*EbN0));

figure;
semilogy(EbN0_db, TEBs, '+');
hold;
plot(EbN0_db, TEB_theorique);
plot(EbN0_db, TEB_theorique_ref);
title("TEB = f(E_b/N_0)");
xlabel("E_b/N_0 (dB)");
ylabel("TEB");
legend("TEB numérique", "TEB théorique", "TEB référence");

%% tracé de plusieurs diagrammes de l'oeil
EbN0_db = linspace(0, 12, 4);
for e=EbN0_db
   
    sigma2_x = P_x * Ns / (2 * log2(M) * 10^(e/10)); % calcul de sigma^2
    X_b = X_c + sqrt(sigma2_x) * randn(1, length(X_c)); % signal bruité
    X_r = filter(h_r, 1, X_b); % signal reçu
    
    figure;
    oeil = reshape(X_r, Ns, length(X_r)/Ns); % permet de tracer le diagramme de l'oeil 
    plot(oeil);
    title("Diagramme de l'oeil, E_b/N_0 = " + num2str(e) + " dB");
    xlabel("Temps (s)");
    ylabel("Amplitude");

end