projet-telecommunications/egalisation.m
2022-09-18 16:39:14 +02:00

222 lines
6 KiB
Matlab

clear;
close all;
Fe = 24000; % fréquence d'échantillonage (Hz)
Te = 1/Fe; % période d'échantillonage (s)
N = 10000; % nombre de bits envoyés
Rb = 3000; % débit binaire
zero_pad = 2^12;
M = 2^1; % signal BPSK
Ts = log2(M)/Rb; % période symbole
Ns = floor(Ts/Te);
T = (0:N*Ns/log2(M)-1) * Te; % échelle temporelle
alpha0 = 1;
alpha1 = 0.5;
h = ones(1, Ns); % mise en forme: réponse impulsionnelle rectangulaire
h_c = [ alpha0 zeros(1, Ns-1) alpha1 ]; % propagation: sélectif
h_r = h; % réception: réponse impulsionnelle rectangulaire
n0 = Ns; % déterminé en traçant le diagramme de l'oeil ou g
%% création de Z et C
ordre = 10;
K = ordre;
bits = [ 1 zeros(1, ordre-1) ]; % dirac
X_k = kron(bits, [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
X_r = filter(h_r, 1, X_c); % signal reçu, sans bruit
X_e = X_r( n0 : Ns : length(bits)*Ns ); % échantillonage du signal reçu
Z = zeros(K, ordre);
for i=1:ordre
Z(i:end, i) = X_e(1:K-i+1);
end
Y_0 = [ 1 zeros(1, ordre-1) ]';
C = Z \ Y_0;
%% tracés réponses en fréquence
X_eg = filter(C, 1, X_e);
g = conv(conv(h, h_c), h_r); % chaine globale
figure;
plot(abs(fft(C, zero_pad)));
title("C");
figure;
plot(abs(fft(g, zero_pad)));
title("G");
figure;
plot(abs(fft(g, zero_pad).*fft(C', zero_pad)));
title("G \times C");
figure;
plot( ifft( abs(fft(g, zero_pad).*fft(C', zero_pad)) ));
title("g * c");
%% tracés réponses impulsionnelles
bits = [ 1 zeros(1, N-1) ];
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
X_r = filter(h_r, 1, X_c); % signal reçu, sans bruit
X_e = X_r( n0 : Ns : length(bits)*Ns ); % échantillonage du signal reçu
X_eg = filter(C, 1, X_e); % égaliseur ZFE
figure;
tiledlayout(3,1);
ax1 = nexttile;
stairs(bits);
legend("Dirac");
ax2 = nexttile;
plot(X_e);
legend("Réponse impulsionnel non égalisée");
ax3 = nexttile;
plot(X_eg);
legend("Réponse impulsionnel égalisée");
linkaxes([ax1 ax2 ax3], 'xy');
axis( axis + [-1 1 -1 1]/5 );
axis
%% constellations
bits = randi([0, 1], 1, N); % bits envoyés
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
X_r = filter(h_r, 1, X_c); % signal reçu, sans bruit
X_e = X_r( n0 : Ns : length(bits)*Ns ); % échantillonage du signal reçu
X_eg = filter(C, 1, X_e); % égaliseur ZFE
figure;
tiledlayout(3,1);
ax1 = nexttile;
plot(complex(X_m), ".", 'MarkerSize', 20);
symbols = unique(X_m(2:end));
for i=1:length(symbols)
text(symbols(i), 0.1, num2str(bits(find(X_m == symbols(i), 1))), 'Color', '#7E2F8E');
end
xlabel("Re");
ylabel("Im");
legend("Symboles théoriques");
grid;
ax2 = nexttile;
plot(complex(X_e(2:end)), ".", 'MarkerSize', 20);
symbols = unique(X_e(2:end));
for i=1:length(symbols)
text(symbols(i), 0.1, num2str(bits(find(X_e == symbols(i), 1))), 'Color', '#7E2F8E');
end
xlabel("Re");
ylabel("Im");
legend("Symboles reçus");
grid;
ax3 = nexttile;
plot(complex(X_eg(2:end)), ".", 'MarkerSize', 20);
symbols = unique(X_eg(2:end));
for i=1:length(symbols)
text(symbols(i), 0.1, num2str(bits(find(X_eg == symbols(i), 1))), 'Color', '#7E2F8E');
end
xlabel("Re");
ylabel("Im");
legend("Symboles égalisés");
grid;
linkaxes([ax1 ax2 ax3], 'xy');
%% tracé de constellations bruitées
P_x = mean(abs(X_c).^2); % puissance du signal en entrée de la reception
EbN0_db = linspace(0, 12, 4);
EbN0 = 10.^(EbN0_db./10);
for e=EbN0
sigma2_x = P_x * Ns / (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:length(bits)*Ns/log2(M) ); % échantillonage du signal reçu
X_eg = filter(C, 1, X_e); % égaliseur ZFE
recu = double( X_eg > 0 ); % on récupère les bits, décision + "démapping"
close all;
figure;
sp1 = subplot(2, 1, 1);
histogram(X_eg, 100);
xlim([-4 4]);
xlabel("Re");
title("Histogramme, $$E_b/N_0 =$$ " + num2str(10*log10(e)) +"dB", 'interpreter','latex');
sp2 = subplot(2, 1, 2);
plot(complex(X_eg(2:10:end)), "*");
hold;
plot(complex(X_m), ".", 'MarkerSize', 20);
xlim([-4 4]);
title("Constellations, $$E_b/N_0 =$$ " + num2str(10*log10(e)) +"dB", 'interpreter','latex');
legend("symboles bruités", "symboles théoriques");
xlabel("Re");
ylabel("Im");
sp2.Position(4) = 0.1;
sp1.Position(2) = 0.28;
sp1.Position(4) = 0.68;
pause(1);
print(gcf,"histcons_egalise_" + num2str(10*log10(e)) + "db.png", '-dpng', '-r300');
end
%% tracé de TEB = f(E_b/N_0)
N = 100000; % nombre de bits envoyés
bits = randi([0, 1], 1, N); % bits envoyés
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
P_x = mean(abs(X_c).^2); % puissance du signal en entrée de la reception
EbN0_db = linspace(0, 10, 200);
EbN0 = 10.^(EbN0_db./10);
TEBs = [];
for e=EbN0
sigma2_x = P_x * Ns / (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:length(bits)*Ns/log2(M) ); % échantillonage du signal reçu
X_eg = filter(C, 1, X_e); % égaliseur ZFE
recu = double( X_eg > 0 ); % on récupère les bits, décision + "démapping"
TEB = mean((recu - bits).^2); % on calcule l'erreur binaire
TEBs = [ TEBs TEB ];
end
TEB_theorique = qfunc(sqrt(EbN0));
close all;
figure;
semilogy(EbN0_db, TEBs, '+');
hold;
plot(EbN0_db, TEB_theorique);
title("TEB $$= f(E_b/N_0)$$", 'interpreter','latex');
xlabel("$$E_b/N_0$$ (dB)", 'interpreter','latex');
ylabel("TEB", 'interpreter','latex');
legend("TEB numérique", "TEB théorique (canal dirac)");