zeefaad/DSP
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

M&M error plot

Browse Source
This commit is contained in:
zeefaad
2025-10-26 11:51:49 +01:00
parent 63f706eee5
commit ae04d5e2bd
5 changed files with 833 additions and 54 deletions
Download Patch File Download Diff File Expand all files Collapse all files

140
QAM/debug.py
View File

@ -5,91 +5,137 @@ import numpy as np
import sys, os
# -------------------- Fichiers de données --------------------
REF_FILE = "constellation_ref.dat"
RX_FILE = "constellation.dat"
ERROR_FILE = "pll_error.dat"
REF_FILE = "constellation_ref.dat"
RX_FILE = "constellation.dat"
PLL_ERROR_FILE = "pll_error.dat"
MM_ERROR_FILE = "mm_timing_error.dat" # Nouveau fichier pour l'erreur M&M
# -------------------- Fonctions de chargement --------------------
def load_points(filename):
if os.path.exists(filename):
return np.loadtxt(filename)
# Utiliser usecols=(0, 1) pour être sûr d'avoir 2 colonnes
return np.loadtxt(filename, usecols=(0, 1))
else:
return np.zeros((0,2))
def load_error(filename):
if os.path.exists(filename):
return np.loadtxt(filename)
# Colonne 0: Indice du symbole, Colonne 1: Valeur d'erreur
return np.loadtxt(filename, usecols=(0, 1))
else:
return np.zeros((0,2))
# -------------------- Application --------------------
app = QtWidgets.QApplication(sys.argv)
win = pg.GraphicsLayoutWidget(show=True, title="Constellation")
win.resize(900, 900)
win = pg.GraphicsLayoutWidget(show=True, title="Synchronisation QAM (M&M)")
win.resize(1200, 800)
win.setBackground('#0a0a0a') # fond noir profond
plot = win.addPlot()
plot.setAspectLocked(True) # axes égaux
# -------------------- Plot 1: Constellation --------------------
# Le premier plot prend la première ligne complète
plot_constel = win.addPlot(title="Constellation Corrigée (Après M&M)")
plot_constel.setAspectLocked(True) # axes égaux
plot_constel.setLabel('left', 'Quadrature (Q)', color='#EEE')
plot_constel.setLabel('bottom', 'In-phase (I)', color='#EEE')
plot_constel.getAxis('left').setPen(pg.mkPen('#AAA', width=2))
plot_constel.getAxis('bottom').setPen(pg.mkPen('#AAA', width=2))
# -------------------- Axes stylés --------------------
plot.getAxis('left').setPen(pg.mkPen('#AAA', width=2))
plot.getAxis('bottom').setPen(pg.mkPen('#AAA', width=2))
plot.getAxis('left').setTextPen(pg.mkPen('#EEE'))
plot.getAxis('bottom').setTextPen(pg.mkPen('#EEE'))
plot.setLabel('left', 'Quadrature (Q)', color='#EEE', size='12pt')
plot.setLabel('bottom', 'In-phase (I)', color='#EEE', size='12pt')
# -------------------- Grille subtile --------------------
grid = pg.GridItem()
grid.setPen(pg.mkPen('#444', width=1, style=QtCore.Qt.DotLine))
plot.addItem(grid)
grid_constel = pg.GridItem()
grid_constel.setPen(pg.mkPen('#444', width=1, style=QtCore.Qt.DotLine))
plot_constel.addItem(grid_constel)
# -------------------- Chargement initial --------------------
ref_data = load_points(REF_FILE)
rx_data = load_points(RX_FILE)
error_data = load_error(ERROR_FILE)
# Points référence
ref_plot = plot.plot(ref_data[:,0], ref_data[:,1],
pen=None,
symbol='o',
symbolSize=10,
symbolBrush=pg.mkBrush('#1E90FF'), # bleu vif
symbolPen=None,
name='Référence')
ref_plot = plot_constel.plot(ref_data[:,0], ref_data[:,1],
pen=None,
symbol='o',
symbolSize=10,
symbolBrush=pg.mkBrush('#1E90FF'), # bleu vif
symbolPen=None,
name='Référence')
# Points reçus
rx_plot = plot.plot(rx_data[:,0], rx_data[:,1],
pen=None,
symbol='x',
symbolSize=6,
symbolBrush=pg.mkBrush('#FF4500'), # orange vif
symbolPen=None,
name='Reçu')
# PLL error en vert
error_plot = plot.plot(error_data[:,0], error_data[:,1],
pen=pg.mkPen('#00FF00', width=2),
symbol=None,
name='PLL error')
rx_plot = plot_constel.plot(rx_data[:,0], rx_data[:,1],
pen=None,
symbol='x',
symbolSize=6,
symbolBrush=pg.mkBrush('#FF4500'), # orange vif
symbolPen=None,
name='Reçu')
# -------------------- Légende stylée --------------------
legend = plot.addLegend()
for item in legend.items:
legend_constel = plot_constel.addLegend()
for item in legend_constel.items:
item[1].setPen(pg.mkPen('#FFF', width=2))
# -------------------- Plot 2: Erreur de Phase PLL --------------------
win.nextRow() # On passe à la deuxième ligne
plot_pll = win.addPlot(title="Erreur de Phase PLL (Boucle de Costas)")
plot_pll.setLabel('left', 'Erreur de Phase Instant. (%)', color='#EEE')
plot_pll.setLabel('bottom', 'Symbole k', color='#EEE')
plot_pll.getAxis('left').setPen(pg.mkPen('#AAA', width=2))
plot_pll.getAxis('bottom').setPen(pg.mkPen('#AAA', width=2))
pll_error_data = load_error(PLL_ERROR_FILE)
pll_error_plot = plot_pll.plot(pll_error_data[:,0], pll_error_data[:,1],
pen=pg.mkPen('#00FF00', width=2), # Vert
symbol=None,
name='Erreur PLL')
# -------------------- Plot 3: Correction Temporelle M&M --------------------
# Ce plot est ajouté SANS win.nextRow() pour le placer à droite du Plot 2 (horizontalement)
plot_mm = win.addPlot(title="Correction Temporelle M&M (Sortie NCO)")
plot_mm.setLabel('left', 'Timing Offset (échantillons)', color='#EEE')
plot_mm.setLabel('bottom', 'Symbole k', color='#EEE')
plot_mm.getAxis('left').setPen(pg.mkPen('#AAA', width=2))
plot_mm.getAxis('bottom').setPen(pg.mkPen('#AAA', width=2))
mm_error_data = load_error(MM_ERROR_FILE)
mm_error_plot = plot_mm.plot(mm_error_data[:,0], mm_error_data[:,1],
pen=pg.mkPen('#FFFF00', width=2), # Jaune
symbol=None,
name='Correction M&M')
# Ligne zéro pour la référence de convergence
plot_mm.addLine(y=0, pen=pg.mkPen('#555', width=1, style=QtCore.Qt.DashLine))
# Ligne de l'offset initial (N/4) pour référence
N_val = 0
if 'qam' in dir():
N_val = qam.N
plot_mm.addLine(y=N_val/4, pen=pg.mkPen('#555', width=1, style=QtCore.Qt.DashLine))
# -------------------- Timer pour mise à jour dynamique --------------------
def update():
global N_val
ref_data = load_points(REF_FILE)
rx_data = load_points(RX_FILE)
error_data = load_error(ERROR_FILE)
pll_error_data = load_error(PLL_ERROR_FILE)
mm_error_data = load_error(MM_ERROR_FILE)
# Constellation
if ref_data.size > 0:
ref_plot.setData(ref_data[:,0], ref_data[:,1])
if rx_data.size > 0:
rx_plot.setData(rx_data[:,0], rx_data[:,1])
if error_data.size > 0:
error_plot.setData(error_data[:,0], error_data[:,1])
# PLL Error
if pll_error_data.size > 0:
pll_error_plot.setData(pll_error_data[:,0], pll_error_data[:,1])
# M&M Timing Correction
if mm_error_data.size > 0:
mm_error_plot.setData(mm_error_data[:,0], mm_error_data[:,1])
# Mise à jour de la ligne de référence N/4 si N est connu
if N_val == 0 and len(mm_error_data) > 0:
# N = Fs * Ts = 44100 * 0.01 = 441. Utiliser cette valeur si non récupérable
N_val = 441 # Valeur par défaut basée sur main.c
plot_mm.clear()
plot_mm.addItem(mm_error_plot)
plot_mm.addLine(y=0, pen=pg.mkPen('#555', width=1, style=QtCore.Qt.DashLine))
plot_mm.addLine(y=N_val/4.0, pen=pg.mkPen('#555', width=1, style=QtCore.Qt.DotLine, dash=[2, 2])) # Ligne de l'offset initial
timer = QtCore.QTimer()
timer.timeout.connect(update)

BIN
QAM/out
View File

Binary file not shown.

21
QAM/qam.c
View File

@ -364,7 +364,7 @@ void pll(qam_system* qam, double complex* s_with_preamble, double complex* r_cor
}
// Synchro temporelle Müller et Müller
void muller_muller_sync(qam_system* qam, double complex* s_data, int nb_symbols, double Kp_mm, double Ki_mm, double complex* r_mm_out) {
void muller_muller_sync(qam_system* qam, double complex* s_data, int nb_symbols, double Kp_mm, double Ki_mm, double complex* r_mm_out, FILE* fp_error) {
double integrator = 0.0;
double timing_offset = 0.0;
double current_idx_double = 0.0;
@ -411,6 +411,11 @@ void muller_muller_sync(qam_system* qam, double complex* s_data, int nb_symbols,
current_idx_double += qam->N - timing_offset;
d_k_minus_1 = d_k;
if (fp_error) {
fprintf(fp_error, "%d % .8f\n", k, 5 * timing_offset);
fflush(fp_error);
}
}
}
@ -474,7 +479,7 @@ int main () {
init_constellation(&qam);
// Conversion du texte en bits
char* texte = "Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux";
char* texte = "Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux,Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux,Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux,Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux,Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux,FIN";
int nb_chars = strlen(texte);
int nb_bits = nb_chars * 8;
int nb_symbols = (nb_bits + qam.k - 1) / qam.k;
@ -530,12 +535,14 @@ int main () {
fclose(fp_pll_error);
double complex* s_bande_base = malloc(sizeof(double complex) * nb_symbols * qam.N);
demodulate_carrier(&qam, s_corrected + L * qam.N, s_bande_base, total_samples);
demodulate_carrier(&qam, s_corrected + L * qam.N, s_bande_base, nb_symbols * qam.N);
double complex* r_mm_out = malloc(sizeof(double complex) * nb_bits);
double Kp_mm = 0.05;
double Ki_mm = 0.001;
muller_muller_sync(&qam, s_bande_base, nb_symbols, Kp_mm, Ki_mm, r_mm_out);
double complex* r_mm_out = malloc(sizeof(double complex) * nb_symbols);
double Kp_mm = 0.5;
double Ki_mm = 0.01;
FILE *fp_mm_error = fopen("mm_timing_error.dat", "w");
muller_muller_sync(&qam, s_bande_base, nb_symbols, Kp_mm, Ki_mm, r_mm_out, fp_mm_error);
fclose(fp_mm_error);
// Démodulation
FILE *fp_constel = fopen("constellation.dat", "w");

150
saveQAM/5/debug.py Normal file
View File

@ -0,0 +1,150 @@
#!/usr/bin/env python3
import pyqtgraph as pg
from pyqtgraph.Qt import QtWidgets, QtCore
import numpy as np
import sys, os
# -------------------- Fichiers de données --------------------
REF_FILE = "constellation_ref.dat"
RX_FILE = "constellation.dat"
PLL_ERROR_FILE = "pll_error.dat"
MM_ERROR_FILE = "mm_timing_error.dat" # Nouveau fichier pour l'erreur M&M
# -------------------- Fonctions de chargement --------------------
def load_points(filename):
if os.path.exists(filename):
# Utiliser usecols=(0, 1) pour être sûr d'avoir 2 colonnes
return np.loadtxt(filename, usecols=(0, 1))
else:
return np.zeros((0,2))
def load_error(filename):
if os.path.exists(filename):
# Colonne 0: Indice du symbole, Colonne 1: Valeur d'erreur
return np.loadtxt(filename, usecols=(0, 1))
else:
return np.zeros((0,2))
# -------------------- Application --------------------
app = QtWidgets.QApplication(sys.argv)
win = pg.GraphicsLayoutWidget(show=True, title="Synchronisation QAM (M&M)")
win.resize(1200, 800)
win.setBackground('#0a0a0a') # fond noir profond
# -------------------- Plot 1: Constellation --------------------
# Le premier plot prend la première ligne complète
plot_constel = win.addPlot(title="Constellation Corrigée (Après M&M)")
plot_constel.setAspectLocked(True) # axes égaux
plot_constel.setLabel('left', 'Quadrature (Q)', color='#EEE')
plot_constel.setLabel('bottom', 'In-phase (I)', color='#EEE')
plot_constel.getAxis('left').setPen(pg.mkPen('#AAA', width=2))
plot_constel.getAxis('bottom').setPen(pg.mkPen('#AAA', width=2))
grid_constel = pg.GridItem()
grid_constel.setPen(pg.mkPen('#444', width=1, style=QtCore.Qt.DotLine))
plot_constel.addItem(grid_constel)
# -------------------- Chargement initial --------------------
ref_data = load_points(REF_FILE)
rx_data = load_points(RX_FILE)
# Points référence
ref_plot = plot_constel.plot(ref_data[:,0], ref_data[:,1],
pen=None,
symbol='o',
symbolSize=10,
symbolBrush=pg.mkBrush('#1E90FF'), # bleu vif
symbolPen=None,
name='Référence')
# Points reçus
rx_plot = plot_constel.plot(rx_data[:,0], rx_data[:,1],
pen=None,
symbol='x',
symbolSize=6,
symbolBrush=pg.mkBrush('#FF4500'), # orange vif
symbolPen=None,
name='Reçu')
# -------------------- Légende stylée --------------------
legend_constel = plot_constel.addLegend()
for item in legend_constel.items:
item[1].setPen(pg.mkPen('#FFF', width=2))
# -------------------- Plot 2: Erreur de Phase PLL --------------------
win.nextRow() # On passe à la deuxième ligne
plot_pll = win.addPlot(title="Erreur de Phase PLL (Boucle de Costas)")
plot_pll.setLabel('left', 'Erreur de Phase Instant. (%)', color='#EEE')
plot_pll.setLabel('bottom', 'Symbole k', color='#EEE')
plot_pll.getAxis('left').setPen(pg.mkPen('#AAA', width=2))
plot_pll.getAxis('bottom').setPen(pg.mkPen('#AAA', width=2))
pll_error_data = load_error(PLL_ERROR_FILE)
pll_error_plot = plot_pll.plot(pll_error_data[:,0], pll_error_data[:,1],
pen=pg.mkPen('#00FF00', width=2), # Vert
symbol=None,
name='Erreur PLL')
# -------------------- Plot 3: Correction Temporelle M&M --------------------
# Ce plot est ajouté SANS win.nextRow() pour le placer à droite du Plot 2 (horizontalement)
plot_mm = win.addPlot(title="Correction Temporelle M&M (Sortie NCO)")
plot_mm.setLabel('left', 'Timing Offset (échantillons)', color='#EEE')
plot_mm.setLabel('bottom', 'Symbole k', color='#EEE')
plot_mm.getAxis('left').setPen(pg.mkPen('#AAA', width=2))
plot_mm.getAxis('bottom').setPen(pg.mkPen('#AAA', width=2))
mm_error_data = load_error(MM_ERROR_FILE)
mm_error_plot = plot_mm.plot(mm_error_data[:,0], mm_error_data[:,1],
pen=pg.mkPen('#FFFF00', width=2), # Jaune
symbol=None,
name='Correction M&M')
# Ligne zéro pour la référence de convergence
plot_mm.addLine(y=0, pen=pg.mkPen('#555', width=1, style=QtCore.Qt.DashLine))
# Ligne de l'offset initial (N/4) pour référence
N_val = 0
if 'qam' in dir():
N_val = qam.N
plot_mm.addLine(y=N_val/4, pen=pg.mkPen('#555', width=1, style=QtCore.Qt.DashLine))
# -------------------- Timer pour mise à jour dynamique --------------------
def update():
global N_val
ref_data = load_points(REF_FILE)
rx_data = load_points(RX_FILE)
pll_error_data = load_error(PLL_ERROR_FILE)
mm_error_data = load_error(MM_ERROR_FILE)
# Constellation
if ref_data.size > 0:
ref_plot.setData(ref_data[:,0], ref_data[:,1])
if rx_data.size > 0:
rx_plot.setData(rx_data[:,0], rx_data[:,1])
# PLL Error
if pll_error_data.size > 0:
pll_error_plot.setData(pll_error_data[:,0], pll_error_data[:,1])
# M&M Timing Correction
if mm_error_data.size > 0:
mm_error_plot.setData(mm_error_data[:,0], mm_error_data[:,1])
# Mise à jour de la ligne de référence N/4 si N est connu
if N_val == 0 and len(mm_error_data) > 0:
# N = Fs * Ts = 44100 * 0.01 = 441. Utiliser cette valeur si non récupérable
N_val = 441 # Valeur par défaut basée sur main.c
plot_mm.clear()
plot_mm.addItem(mm_error_plot)
plot_mm.addLine(y=0, pen=pg.mkPen('#555', width=1, style=QtCore.Qt.DashLine))
plot_mm.addLine(y=N_val/4.0, pen=pg.mkPen('#555', width=1, style=QtCore.Qt.DotLine, dash=[2, 2])) # Ligne de l'offset initial
timer = QtCore.QTimer()
timer.timeout.connect(update)
timer.start(500) # ms
# -------------------- Anti-aliasing --------------------
pg.setConfigOptions(antialias=True)
# -------------------- Exécution --------------------
if __name__ == '__main__':
sys.exit(app.exec_())

576
saveQAM/5/qam.c Normal file
View File

@ -0,0 +1,576 @@
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <string.h>
#define A 10
struct qam_system_s {
int M; // Nombre de symboles M-QAM
int k; // Nombre de bits/symboles
int sm; // sqrt M
double Fs; // Fréquence d'échantillionage
double Ts; // Temps d'échantillionage
int N; // Nombre d'échantillions
double Fc; // Fréquence de la porteuse
double complex* constellation; // Tableau de symboles I + j Q 1D
};
typedef struct qam_system_s qam_system;
// Initialisation de la constellation (double tableau de taille sqrt(M)),
void init_constellation (qam_system* qam) {
qam->constellation = (double complex*)malloc(sizeof(double complex) * qam->M);
double norm_factor = sqrt((double)(qam->M - 1) / 3.0); // Pour puissance unitaire
for (int i = 0; i < qam->sm; i++) {
double ip = -(qam->sm - 1) + 2 * i;
for (int j = 0; j < qam->sm; j++) {
double qp = -(qam->sm - 1) + 2 * j;
int idx = i * qam->sm + j;
qam->constellation[idx] = (ip + I * qp) / norm_factor;
}
}
}
int bin_to_gray (int x) {
return x ^ (x >> 1);
}
// Iteration de XOR
int gray_to_bin(int g) {
int b = g;
while (g >>= 1)
b ^= g;
return b;
}
// Modulation QAM
void modulate (qam_system* qam, double complex* symbols, int nb_symbols, double complex* s) {
for (int k = 0; k < nb_symbols; k++) {
double complex iq = symbols[k];
for (int n = 0; n < qam->N; n++) {
int idx = k * qam->N + n;
s[idx] = A * iq * cexp(2 * I * M_PI * qam->Fc * ((double)idx / qam->Fs));
}
}
}
// Demodulation QAM
void demodulate(qam_system* qam, double complex* s, int nb_symbols, uint8_t* bits_hat, FILE *fp_constel) {
for (int k = 0; k < nb_symbols; k++) {
double complex r = 0;
for (int n = 0; n < qam->N; n++) {
r += s[k * qam->N + n] * cexp(-2 * I * M_PI * qam->Fc * ((double)(k * qam->N + n) / qam->Fs)) / A;
}
r /= qam->N;
if (fp_constel) {
fprintf(fp_constel, "% .8f % .8f\n", creal(r), cimag(r));
fflush(fp_constel);
}
// Distance euclidien de Ir et Qr pour avoir le point le plus proche de la constellation (lent)
double min_d = INFINITY;
int i_cl = 0;
int j_cl = 0;
for (int idx = 0; idx < qam->M; idx++) {
double d = cabs(r - qam->constellation[idx]);
if (d < min_d) {
min_d = d;
i_cl = idx / qam->sm;
j_cl = idx % qam->sm;
}
}
// Gray mapping
if (qam->k % 2 != 0) {
printf("demodulate : k pair (k = %d)\n", qam->k);
exit(1);
}
int bits_per_axis = qam->k / 2;
int i_bin = gray_to_bin(i_cl);
int j_bin = gray_to_bin(j_cl);
for (int b = 0; b < bits_per_axis; b++) {
bits_hat[k * qam->k + b] = (i_bin >> (bits_per_axis - 1 - b)) & 1;
bits_hat[k * qam->k + bits_per_axis + b] = (j_bin >> (bits_per_axis - 1 - b)) & 1;
}
}
}
// Calcul du bruit gaussien pour un sigma donné
// Formule de Box-Muller
double gaussian_noise (double sigma) {
double u1 = (rand() + 1) / ((double)RAND_MAX + 2);
double u2 = (rand() + 1) / ((double)RAND_MAX + 2);
return sigma * sqrt(-2 * log(u1)) * cos(2 * M_PI * u2);
}
// Ajout du bruit
void add_noise (double complex* s, int len, double sigma) {
//double signal_power = (2.0/3.0)*(qam.M-1);
//double snr_dB = 5; // SNR en dB
//double snr_lin = pow(10.0, snr_dB / 10.0);
//double sigma = sqrt(signal_power / snr_lin);
for (int i = 0; i < len; i++) {
double nr = gaussian_noise(sigma);
double ni = gaussian_noise(sigma);
s[i] += nr + I * ni;
}
}
// Changer le tableau de bits en boolen ou alors la represenation binaire et shifter pour extraire les bits (pas bien si M plus grand)
void bits_to_symbols (qam_system* qam, uint8_t* bits, int nb_bits, double complex* symbols) {
int nb_symbols = nb_bits / qam->k;
// k pair
if (qam->k % 2 != 0) {
printf("bits_to_symbols : k doit être pair (k = %d) \n", qam->k);
exit(1);
}
int bits_per_axis = qam->k / 2;
for (int sym = 0; sym < nb_symbols; sym++) {
// Construire les indices binaires i_bin (MSB...) et j_bin (LSB...)
int i_bin = 0;
int j_bin = 0;
for (int b = 0; b < bits_per_axis; b++) {
i_bin = (i_bin << 1) | bits[sym * qam->k + b];
j_bin = (j_bin << 1) | bits[sym * qam->k + bits_per_axis + b];
}
int i_gray = bin_to_gray(i_bin);
int j_gray = bin_to_gray(j_bin);
if (i_gray < 0 || i_gray >= qam->sm || j_gray < 0 || j_gray >= qam->sm) {
printf("bits_to_symbols : IOOR i_gray = %d j_gray = %d sm = %d \n", i_gray, j_gray, qam->sm);
exit(1);
}
symbols[sym] = qam->constellation[i_gray * qam->sm + j_gray];
}
}
double compare_bits(uint8_t* bits1, uint8_t* bits2, int nb_bits) {
int errors = 0;
for (int i = 0; i < nb_bits; i++) {
if (bits1[i] != bits2[i]) errors++;
}
return (double)errors / nb_bits;
}
void add_dephasage(double complex* s, double phi_offset, int total_samples) {
for (int i = 0; i < total_samples; i++) {
s[i] *= cexp(I * phi_offset);
}
}
void add_freq(qam_system* qam, double complex* s, double freq_offset, int total_samples) {
for (int i = 0; i < total_samples; i++) {
double t = (double)i / qam->Fs;
s[i] *= cexp(I * 2 * M_PI * freq_offset * t);
}
}
void fill_constellation_data(qam_system* qam, FILE *fp_ref) {
for (int i = 0; i < qam->sm; i++) {
for (int j = 0; j < qam->sm; j++) {
int idx = i * qam->sm + j;
fprintf(fp_ref, "% .8f % .8f\n", creal(qam->constellation[idx]), cimag(qam->constellation[idx]));
}
}
}
void reconstruction_text(int nb_chars, uint8_t* output_bits, char* texte_recup) {
for(int i = 0; i < nb_chars; i++){
char c = 0;
for(int b = 0; b < 8; b++){
c |= output_bits[i*8 + b] << (7-b);
}
texte_recup[i] = c;
}
texte_recup[nb_chars] = '\0';
}
void text_to_bits(int nb_chars, int nb_bits, char* texte, uint8_t* input_bits) {
for(int i = 0; i < nb_chars; i++){
for(int b = 0; b < 8; b++){
input_bits[i*8 + b] = (texte[i] >> (7-b)) & 1;
}
}
}
// Libération de la mémoire
void free_constellation(qam_system* qam) {
free(qam->constellation);
}
// Préambule QAM
void generate_preamble(qam_system* qam, double complex* preamble, int L) {
// L 1er symboles de la constellation
for (int i = 0; i < L; i++) {
preamble[i] = qam->constellation[i % qam->M];
}
}
// Concatène le préambule avec le signal modulé
double complex* concat_preamble_signal(double complex* preamble_mod, int preamble_len, double complex* s_mod, int nb_symbols, int N) {
int total_samples = (preamble_len + nb_symbols) * N;
double complex* s_concat = malloc(sizeof(double complex) * total_samples);
// Copier le préambule modulé
for (int i = 0; i < preamble_len * N; i++) {
s_concat[i] = preamble_mod[i];
}
// Copier le signal modulé après le préambule
for (int i = 0; i < nb_symbols * N; i++) {
s_concat[preamble_len * N + i] = s_mod[i];
}
return s_concat;
}
// Coarse Frequency Offset avec préambule avec N = 1 (lag 1)
void cfo(double complex* s_with_preamble, int N, int L, double Fs, int total_samples, double Fc) {
double complex sum = 0;
int lag = 5;
for (int n = lag; n < L * N; n++) {
sum += s_with_preamble[n] * conj(s_with_preamble[n - lag]);
}
double phi = carg(sum);
double f_est = (phi / (2.0 * M_PI * lag)) * Fs;
double f_offset = f_est - Fc;
printf("CFO estimé : %f Hz \n", f_offset);
for (int n = 0; n < total_samples; n++) {
s_with_preamble[n] *= cexp(-I * 2.0 * M_PI * f_offset * n / Fs);
}
}
// Fine Frequency Offset (FFO) Correction
void ffo(qam_system* qam, double complex* s_with_preamble, double complex* preamble_ref, int L, int total_samples) {
double complex r_preamble[L];
for (int k = 0; k < L; k++) {
double complex r = 0;
for (int n = 0; n < qam->N; n++) {
int idx = k * qam->N + n;
r += s_with_preamble[idx] * cexp(-2 * I * M_PI * qam->Fc * ((double)idx / qam->Fs)) / A;
}
r /= qam->N;
r_preamble[k] = r;
}
double complex diff_phase_sum = 0;
double Ts_symbole = qam->N / qam->Fs;
for (int k = 1; k < L; k++) {
double complex error_k = r_preamble[k] * conj(preamble_ref[k]);
double complex error_k_minus_1 = r_preamble[k - 1] * conj(preamble_ref[k - 1]);
diff_phase_sum += error_k * conj(error_k_minus_1);
}
double phi_sym = carg(diff_phase_sum);
double delta_f_fine_est = phi_sym / (2.0 * M_PI * Ts_symbole);
printf("FFO estimé (Delta f fine) : %f Hz\n", delta_f_fine_est);
for (int n = 0; n < total_samples; n++) {
s_with_preamble[n] *= cexp(-I * 2.0 * M_PI * delta_f_fine_est * n / qam->Fs);
}
}
// Phase Offset (PO) Correction
void po(qam_system* qam, double complex* s_with_preamble, double complex* preamble_ref, int L, int total_samples) {
double complex r_preamble[L];
for (int k = 0; k < L; k++) {
double complex r = 0;
for (int n = 0; n < qam->N; n++) {
int idx = k * qam->N + n;
r += s_with_preamble[idx] * cexp(-2 * I * M_PI * qam->Fc * ((double)idx / qam->Fs)) / A;
}
r /= qam->N;
r_preamble[k] = r;
}
double complex phase_error_sum = 0;
for (int k = 0; k < L; k++) {
phase_error_sum += r_preamble[k] * conj(preamble_ref[k]);
}
double phi_est = carg(phase_error_sum);
printf("Phase Offset (PO) estimée : %f radians (soit %f degrés)\n", phi_est, phi_est * 180.0 / M_PI);
double complex phase_corr = cexp(-I * phi_est);
for (int n = 0; n < total_samples; n++) {
s_with_preamble[n] *= phase_corr;
}
}
// Costas loop
void pll(qam_system* qam, double complex* s_with_preamble, double complex* r_corr, int L, int total_samples, int nb_symbols, double Kp, double Ki, double alpha, FILE* fp_error) {
double phase_est = 0.0;
double integrator = 0.0;
double filtered_error = 0.0;
int start_data_idx = L * qam->N;
int N = qam->N;
for (int k = 0; k < nb_symbols; k++) {
double complex r_k = 0;
for (int n = 0; n < N; n++) {
int idx = start_data_idx + k * N + n;
double t = (double)idx / qam->Fs;
r_k += s_with_preamble[idx] * cexp(-2.0 * I * M_PI * qam->Fc * t) / A;
}
double complex r_symbol = r_k / N;
r_symbol *= cexp(-I * phase_est);
double min_dist = INFINITY;
double complex closest = 0;
for (int idx = 0; idx < qam->M; idx++) {
double dist = cabs(r_symbol - qam->constellation[idx]);
if (dist < min_dist) {
min_dist = dist;
closest = qam->constellation[idx];
}
}
double error = carg(r_symbol * conj(closest));
filtered_error = (1.0 - alpha) * filtered_error + alpha * error;
integrator += Ki * filtered_error;
phase_est += Kp * filtered_error + integrator;
if (fp_error) {
fprintf(fp_error, "%d % .8f\n", k, 100 * filtered_error);
fflush(fp_error);
}
for (int n = 0; n < N; n++) {
int idx = start_data_idx + k * N + n;
r_corr[idx] = s_with_preamble[idx] * cexp(-I * phase_est);
}
}
}
// Synchro temporelle Müller et Müller
void muller_muller_sync(qam_system* qam, double complex* s_data, int nb_symbols, double Kp_mm, double Ki_mm, double complex* r_mm_out, FILE* fp_error) {
double integrator = 0.0;
double timing_offset = 0.0;
double current_idx_double = 0.0;
double complex d_k_minus_1 = 0;
for (int k = 0; k < nb_symbols; k++) {
int opt_idx = (int)round(current_idx_double);
int tilde_idx = (int)round(current_idx_double + qam->N / 2.0);
if (opt_idx >= nb_symbols * qam->N || tilde_idx >= nb_symbols * qam->N || tilde_idx < 0)
break;
double complex r_k = s_data[opt_idx];
double complex r_tilde_k = s_data[tilde_idx];
double min_d = INFINITY;
double complex d_k = 0;
int decision_idx = 0;
for (int idx = 0; idx < qam->M; idx++) {
double d = cabs(r_k - qam->constellation[idx]);
if (d < min_d) {
min_d = d;
d_k = qam->constellation[idx];
decision_idx = idx;
}
}
r_mm_out[k] = r_k;
// Erreur M&M
double error_k = 0.0;
if (k > 0) {
// Formule M&M: e_k = Re{ r_tilde_k * (d_{k-1}^* - d_k^*) }
double complex error_term = r_tilde_k * conj(d_k_minus_1 - d_k);
error_k = creal(error_term);
}
// Costas Loop PI
integrator += Ki_mm * error_k;
timing_offset = Kp_mm * error_k + integrator;
current_idx_double += qam->N - timing_offset;
d_k_minus_1 = d_k;
if (fp_error) {
fprintf(fp_error, "%d % .8f\n", k, 5 * timing_offset);
fflush(fp_error);
}
}
}
// Demodulation M&M
void demodulate_sync_adapted(qam_system* qam, double complex* r_mm_out, int nb_symbols, uint8_t* bits_hat, FILE *fp_constel) {
for (int k = 0; k < nb_symbols; k++) {
double complex r = r_mm_out[k];
if (fp_constel) {
fprintf(fp_constel, "% .8f % .8f\n", creal(r), cimag(r));
fflush(fp_constel);
}
// Distance euclidien (quantification/décision)
double min_d = INFINITY;
int i_cl = 0;
int j_cl = 0;
for (int idx = 0; idx < qam->M; idx++) {
double d = cabs(r - qam->constellation[idx]);
if (d < min_d) {
min_d = d;
i_cl = idx / qam->sm;
j_cl = idx % qam->sm;
}
}
// Gray mapping et extraction des bits
if (qam->k % 2 != 0) {
printf("demodulate : k doit être pair (k = %d)\n", qam->k);
exit(1);
}
int bits_per_axis = qam->k / 2;
int i_bin = gray_to_bin(i_cl);
int j_bin = gray_to_bin(j_cl);
for (int b = 0; b < bits_per_axis; b++) {
bits_hat[k * qam->k + b] = (i_bin >> (bits_per_axis - 1 - b)) & 1;
bits_hat[k * qam->k + bits_per_axis + b] = (j_bin >> (bits_per_axis - 1 - b)) & 1;
}
}
}
void demodulate_carrier(qam_system* qam, double complex* s_input, double complex* s_bandebase, int total_samples) {
for (int n = 0; n < total_samples; n++) {
double t = (double)n / qam->Fs;
s_bandebase[n] = s_input[n] * cexp(-I * 2 * M_PI * qam->Fc * t) / A;
}
}
int main () {
// Initialisation du system qam
qam_system qam;
qam.M = 16;
qam.k = (int)log2((double)(qam.M));
qam.sm = (int)sqrt(qam.M);
qam.Fs = 44100;
qam.Ts = 0.01;
qam.N = (int)(qam.Fs * qam.Ts);
qam.Fc = 2000;
init_constellation(&qam);
// Conversion du texte en bits
char* texte = "Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux,Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux,Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux,Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux,Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux, Vif juge, trempez ce blond whisky aqueux,FIN";
int nb_chars = strlen(texte);
int nb_bits = nb_chars * 8;
int nb_symbols = (nb_bits + qam.k - 1) / qam.k;
uint8_t* input_bits = malloc(nb_bits * sizeof(uint8_t));
text_to_bits(nb_chars, nb_bits, texte, input_bits);
// Conversion en symboles
double complex* symbols = malloc(sizeof(double complex) * nb_symbols);
bits_to_symbols(&qam, input_bits, nb_bits, symbols);
// Initialisation du préambule
int L = 15;
int total_samples = qam.N * (nb_symbols + L);
double complex* preamble = malloc(sizeof(double complex) * L);
generate_preamble(&qam, preamble, L);
double complex* preamble_mod = malloc(sizeof(double complex) * L * qam.N);
modulate(&qam, preamble, L, preamble_mod);
double complex* s_mod = malloc(sizeof(double complex) * nb_symbols * qam.N);
modulate(&qam, symbols, nb_symbols, s_mod);
double complex* s = malloc(sizeof(double complex) * total_samples);
double complex* s_with_preamble = concat_preamble_signal(preamble_mod, L, s_mod, nb_symbols, qam.N);
// Ajout du bruit
add_noise(s_with_preamble, total_samples, 2);
FILE *fp_ref = fopen("constellation_ref.dat", "w");
fill_constellation_data(&qam, fp_ref);
fclose(fp_ref);
// Ajout de dephasage
add_dephasage(s_with_preamble, M_PI / 2.0 , total_samples);
// AJout de decalage de fréquence
add_freq(&qam, s_with_preamble, 100, total_samples);
// Estimation / correction du CFO
cfo(s_with_preamble, qam.N, L, qam.Fs, total_samples, qam.Fc);
ffo(&qam, s_with_preamble, preamble, L, total_samples);
// Correction phase
po(&qam, s_with_preamble, preamble, L, total_samples);
// PLL
double complex* s_corrected = malloc(sizeof(double complex) * total_samples);
double Kp = 0.03;
double Ki = 0.002;
double alpha = 0.1;
FILE *fp_pll_error = fopen("pll_error.dat", "w");
pll(&qam, s_with_preamble, s_corrected, L, total_samples, nb_symbols, Kp, Ki, alpha, fp_pll_error);
fclose(fp_pll_error);
double complex* s_bande_base = malloc(sizeof(double complex) * nb_symbols * qam.N);
demodulate_carrier(&qam, s_corrected + L * qam.N, s_bande_base, nb_symbols * qam.N);
double complex* r_mm_out = malloc(sizeof(double complex) * nb_symbols);
double Kp_mm = 0.5;
double Ki_mm = 0.01;
FILE *fp_mm_error = fopen("mm_timing_error.dat", "w");
muller_muller_sync(&qam, s_bande_base, nb_symbols, Kp_mm, Ki_mm, r_mm_out, fp_mm_error);
fclose(fp_mm_error);
// Démodulation
FILE *fp_constel = fopen("constellation.dat", "w");
uint8_t* output_bits = (uint8_t*)malloc(nb_bits * sizeof(uint8_t));
demodulate_sync_adapted(&qam, r_mm_out, nb_symbols, output_bits, fp_constel);
//demodulate(&qam, s_corrected + L * qam.N, nb_symbols, output_bits, fp_constel);
fclose(fp_constel);
// Reconstruction du texte
char* texte_recup = malloc(nb_chars + 1);
reconstruction_text(nb_chars, output_bits, texte_recup);
//printf("Texte original : %s\n\n", texte);
printf("Texte demodulé : %s\n", texte_recup);
// Calcul du BER
double ber = compare_bits(input_bits, output_bits, nb_bits);
printf("Taux d'erreur QAM: %.4f\n", ber * 100);
// Libération mémoire
free(input_bits);
free(output_bits);
free(symbols);
free(preamble);
free(preamble_mod);
free(s_with_preamble);
free(texte_recup);
free_constellation(&qam);
return 0;
}