TP-telecommunications/TP4/code/QPSK.m
2023-06-20 21:23:24 +02:00

121 lines
3.8 KiB
Matlab

clear;
close all;
Fe = 12000; % fréquence d'échantillonage (Hz)
Te = 1/Fe; % période d'échantillonage (s)
N = 50000; % nombre de bits envoyés
Rb = 6000; % débit binaire
Fp = 2000; % fréquence porteuse
load('bits');
%bits = randi([0, 1], 1, N); % bits envoyés
M = 2^2; % signal 16-QAM
Ts = log2(M)/Rb; % période symbole
Ns = floor(Ts/Te);
T = (0:N*Ns/log2(M)-1) * Te; % échelle temporelle
fc = 1.25*Fe; % fréquence de coupure/
ordre = 100; % ordre du filtre
alpha = 0.35; % roll off
span = 10; % largeur du racine de cosinus surélevé
h = rcosdesign(alpha, span, Ns); % mise en forme: racine de cosinus surélevé
h_c = [ 1 zeros(1, Ns-1) ]; % propagation: dirac
h_r = h; % réception: racine de cosinus surélevé
g = conv(conv(h, h_c), h_r); % réponse impulsionnelle globale
n0 = 1; % déterminé en traçant le diagramme de l'oeil de real(X_env)
X_regrp = reshape(bits, log2(M), N/log2(M)).'; % on regroupe les bits
X_co = qammod(bi2de(X_regrp), M, 'gray');
X_kr = kron(X_co.', [1 zeros(1, Ns-1)]); % Suréchantillonnage
X_env = [ X_kr zeros(1, span*Ns/2) ]; % on rajoute des zéros pour le retard
X_env = filter(h, 1, X_env); % enveloppe complexe
X_env = X_env(span*Ns/2 + 1 : end); % on dégage les zéros
X_tran = filter(h_c, 1, X_env); % signal transmis -> passage dans le canal
DSP = pwelch(real(X_env));
P_x = mean(abs(X_env).^2);
EbN0_db = linspace(0, 8, 200);
EbN0 = 10.^(EbN0_db./10);
TEBs = [];
TESs = [];
for e=EbN0
sigma2_x = P_x * Ns / (2 * log2(M) * e); % calcul de sigma^2
X_bruit = sqrt(sigma2_x) * randn(2, length(X_tran)); % bruit
X_bruit = X_tran + X_bruit(1,:) + j * X_bruit(2,:); % signal bruité
X_recu = [ X_bruit zeros(1, span*Ns/2) ]; % on rajoute des zéros pour le retard
X_recu = filter(h_r, 1, X_recu); % signal reçu
X_recu = X_recu(span*Ns/2 + 1:end); % on dégage les zéros
X_sample = X_recu( n0 : Ns : N*Ns/log2(M) ); % signal échantillonné
X_demod = qamdemod(X_sample, M, 'gray'); % demodulation
recu = de2bi(X_demod, log2(M)); % on reconstitue les bits à partir des symbols
TES = mean((bi2de(recu) - bi2de(X_regrp)).^2); % on calcule l'erreur symbole
recu = reshape(recu', size(bits)); % on transforme les symboles en bits
TEB = mean((recu - bits).^2); % on calcule l'erreur binaire
TESs = [ TESs TES ];
TEBs = [ TEBs TEB ];
end
TEB_theorique = qfunc(sqrt(2*EbN0));
figure;
semilogy(EbN0_db, TESs, '+');
hold;
semilogy(EbN0_db, TEBs, '+');
plot(EbN0_db, TEB_theorique);
title("TEB = f(E_b/N_0)");
xlabel("E_b/N_0 (dB)");
ylabel("erreur (dB)");
legend("TES numérique", "TEB numérique", "TEB théorique");
save QPSK TEBs DSP
save EbN0_db EbN0_db
%%
EbN0_db = [2 4 6 8];
EbN0 = 10.^(EbN0_db./10);
for e=EbN0
sigma2_x = P_x * Ns / (2 * log2(M) * e); % calcul de sigma^2
X_bruit = sqrt(sigma2_x) * randn(2, length(X_tran)); % bruit
X_bruit = X_tran + X_bruit(1,:) + j * X_bruit(2,:); % signal bruité
X_recu = [ X_bruit zeros(1, span*Ns/2) ]; % on rajoute des zéros pour le retard
X_recu = filter(h_r, 1, X_recu); % signal reçu
X_recu = X_recu(span*Ns/2 + 1:end); % on dégage les zéros
X_sample = X_recu( n0 : Ns : N*Ns/log2(M) ); % signal échantillonné
X_demod = qamdemod(X_sample, M, 'gray'); % demodulation
recu = de2bi(X_demod, log2(M)); % on reconstitue les bits à partir des symbols
TES = mean((bi2de(recu) - bi2de(X_regrp)).^2); % on calcule l'erreur symbole
recu = reshape(recu', size(bits)); % on transforme les symboles en bits
TEB = mean((recu - bits).^2); % on calcule l'erreur binaire
figure;
plot(X_sample, '+');
hold;
plot(complex(X_co), '.', 'MarkerSize', 25);
axis([-2, 2, -2, 2]);
title(sprintf("Eb/N0 = %.2fdB, erreur = %.2f%s", 10*log10(e), TEB*100, "%"));
legend("X\_sample", "X\_co");
axis square;
end