PLL implementation

QAM/<

@@ -0,0 +1,291 @@
+#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
+    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
+};
+typedef struct qam_system_s qam_system;
+
+// Initialisation de la constellation (double tableau de taille sqrt(M)), 
+// ToDo : changer à un tableau à 1 dimension pour éviter de calculer sqrt(M)
+void init_constellation (qam_system* qam) {
+    int sm = (int)sqrt(qam->M);
+    qam->constellation = (double complex**)malloc(sizeof(double complex*) * sm);
+    
+    for (int i = 0; i < sm; i++) {
+       qam->constellation[i] = (double complex*)malloc(sizeof(double complex) * sm);
+    }
+
+    double norm_factor = sqrt((double)(qam->M - 1) / 3.0); // Pour puissance unitaire
+
+    for (int i = 0; i < sm; i++) {
+        double complex ip = -(sm - 1) + 2 * i;
+        for (int j = 0; j < sm; j++) {
+            double complex qp = -(sm - 1) + 2 * j;
+            qam->constellation[i][j] = (ip + I * qp) / norm_factor; 
+        }
+    }
+}
+
+// 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) {
+    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;
+    int sm = sqrt(qam->M);
+    for (int k = 0; k < nb_symbols; k++) {
+        int id = 0;
+        for (int b = 0 ; b < qam->k; b++) {
+            id = id * 2 + bits[k * qam->k + b];
+        }
+        int i = id / sm; 
+        int j = id % sm;
+        symbols[k] = qam->constellation[i][j];
+    }
+}
+
+// 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)
+        int sm = (int)sqrt(qam->M);
+        double min_d = INFINITY;
+        int i_cl = 0;
+        int j_cl = 0;
+        for (int i = 0; i < sm; i++) {
+            for (int j = 0; j < sm; j++) {
+                double d = cabs(r - qam->constellation[i][j]);
+                if (d < min_d) {
+                    min_d = d;
+                    i_cl = i;
+                    j_cl = j;
+                }
+            }
+        }
+
+        // index du symbole (id) : même mappage que dans bits_to_symbols()
+        int id = i_cl * sm + j_cl;
+
+        for (int b = 0; b < qam->k; b++) {
+            bits_hat[k * qam->k + b] = (id >> (qam->k - 1 - b)) & 1;
+        }
+    }
+}
+
+// PLL pour corriger le déphasage
+// Entrées : signal reçu s, nombre total d'échantillons, gain proportionnel Kp, gain intégral Ki
+// Sortie : signal corrigé r_corr
+void fpll_mqam(qam_system* qam, double complex* s, double complex* r_corr, int total_samples, double Kp, double Ki) {
+    double phase_est = 0.0;
+    double freq_correction = 0.0;
+    int sm = (int)sqrt(qam->M);
+
+    for (int n = 0; n < total_samples; n++) {
+        // Corrige le signal avec l'estimation de phase courante
+        r_corr[n] = s[n] * cexp(-I * phase_est);
+
+        // Extraire le symbole correspondant au "échantillon central" du symbole
+        if (n % qam->N == qam->N/2) { 
+            // Trouver le symbole quantifié le plus proche
+            double complex r = r_corr[n];
+            double min_d = INFINITY;
+            double complex closest = 0;
+            for (int i = 0; i < sm; i++) {
+                for (int j = 0; j < sm; j++) {
+                    double d = cabs(r - qam->constellation[i][j]);
+                    if (d < min_d) {
+                        min_d = d;
+                        closest = qam->constellation[i][j];
+                    }
+                }
+            }
+            // Calculer l'erreur de phase
+            double error = carg(r * conj(closest));
+
+            // Mettre à jour la PLL
+            freq_correction += Ki * error;
+            phase_est += Kp * error + freq_correction;
+        }
+    }
+}
+
+
+// Libération de la mémoire
+void free_constellation(qam_system* qam) {
+    int sm = (int)sqrt(qam->M);
+    for (int i = 0; i < sm; i++)
+        free(qam->constellation[i]);
+    free(qam->constellation);
+}
+
+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;
+}
+
+int main () {
+    qam_system qam;
+    qam.M = 4;
+    qam.k = (int)log2((double)(qam.M));
+    qam.Fs = 44100;
+    //qam.Ts = 0.0003;
+    //qam.N = (int)qam.Fs * qam.Ts;
+    qam.Ts = 0.01;
+    qam.N = (int)(qam.Fs * qam.Ts);
+    qam.Fc = 2000;
+    init_constellation(&qam);
+
+    //int nb_bits = 1000;
+    //int nb_symbols = nb_bits / qam.k;
+
+    //uint8_t* input_bits = malloc(nb_bits * sizeof(uint8_t));
+    //for (int i = 0; i < nb_bits; i++) {
+    //    input_bits[i] = rand() % 2;
+    //}
+    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, ";  
+    int nb_chars = strlen(texte);
+    int nb_bits = nb_chars * 8;
+    int nb_symbols = (nb_bits + qam.k - 1) / qam.k;
+
+    // Conversion du texte en bits
+    uint8_t* input_bits = malloc(nb_bits * sizeof(uint8_t));
+    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;
+        }
+    }
+    
+    // Conversion en symboles
+    double complex* symbols = malloc(sizeof(double complex) * nb_symbols);
+    bits_to_symbols(&qam, input_bits, nb_bits, symbols);
+
+    // Modulation
+    int total_samples = qam.N * nb_symbols;
+    double complex* s = malloc(sizeof(double complex) * total_samples);
+    modulate(&qam, symbols, nb_symbols, s);
+
+    // Ajout du bruit
+    double signal_power = (2.0/3.0)*(qam.M-1); // puissance moyenne
+    double snr_dB = 5; // SNR en dB
+    double snr_lin = pow(10.0, snr_dB / 10.0);
+    double sigma = sqrt(signal_power / snr_lin);
+    add_noise(s, total_samples, 0);
+
+
+    FILE *fp_ref = fopen("constellation_ref.dat", "w");
+    int sm = (int)sqrt(qam.M);
+    for (int i = 0; i < sm; i++) {
+        for (int j = 0; j < sm; j++) {
+            fprintf(fp_ref, "% .8f % .8f\n", creal(qam.constellation[i][j]), cimag(qam.constellation[i][j]));
+        }
+    }
+    fclose(fp_ref);
+    FILE *fp_constel = fopen("constellation.dat", "w");
+
+    // Ajout de dephasage
+    //double phase_offset = M_PI / 6.0;  // 30 degrés
+    //for (int i = 0; i < total_samples; i++) {
+    //    s[i] *= cexp(I * phase_offset);
+    //}
+    
+    // AJout de decalage de fréquence
+    double freq_offset = 0; // Hz de décalage
+    for (int i = 0; i < total_samples; i++) {
+        double t = (double)i / qam.Fs;
+        s[i] *= cexp(I * 2 * M_PI * freq_offset * t);
+    }
+    
+    // Ajout de decalage entre les symbole
+    //int offset_samples = (int)(0.3 * qam.N); // décalage de 30% d’un symbole
+    //memmove(s + offset_samples, s, (total_samples - offset_samples) * sizeof(double complex));
+
+    double complex* r_corr = malloc(sizeof(double complex) * total_samples);
+    double Kp = 0.001;
+    double Ki = 0.0001;
+    fpll_mqam(&qam, s, r_corr, total_samples, Kp, Ki);
+
+    // Démodulation
+    uint8_t* output_bits = (uint8_t*)malloc(nb_bits * sizeof(uint8_t));
+    demodulate(&qam, r_corr, nb_symbols, output_bits, fp_constel);
+
+    fclose(fp_constel);
+
+    // Reconstruction du texte
+    char* texte_recup = malloc(nb_chars + 1);
+    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';
+    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 blind QAM: %.4f\n", ber * 100);
+
+    // Libération mémoire
+    free(input_bits);
+    free(output_bits);
+    free(symbols);
+    free(r_corr);
+    free(s);
+    free_constellation(&qam);
+
+    return 0;
+}

QAM/OLD/constellation_ref.dat

@@ -0,0 +1,16 @@
+-1.34164079 -1.34164079
+-1.34164079 -0.44721360
+-1.34164079  0.44721360
+-1.34164079  1.34164079
+-0.44721360 -1.34164079
+-0.44721360 -0.44721360
+-0.44721360  0.44721360
+-0.44721360  1.34164079
+ 0.44721360 -1.34164079
+ 0.44721360 -0.44721360
+ 0.44721360  0.44721360
+ 0.44721360  1.34164079
+ 1.34164079 -1.34164079
+ 1.34164079 -0.44721360
+ 1.34164079  0.44721360
+ 1.34164079  1.34164079

QAM/OLD/debug.py

@@ -0,0 +1,32 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+def plot_constellations(ref_file, rx_file, title="Constellation comparison"):
+    # Charger et forcer 2D
+    ref_data = np.atleast_2d(np.loadtxt(ref_file))
+    rx_data = np.atleast_2d(np.loadtxt(rx_file))
+    
+    x_ref, y_ref = ref_data[:,0], ref_data[:,1]
+    x_rx, y_rx = rx_data[:,0], rx_data[:,1]
+
+    plt.figure(figsize=(6,6))
+    plt.scatter(x_ref, y_ref, color='blue', s=50, marker='o', label='Référence')
+    plt.scatter(x_rx, y_rx, color='red', s=50, marker='x', label='Reçu')
+
+    # Ajustement automatique des limites
+    all_x = np.concatenate([x_ref, x_rx])
+    all_y = np.concatenate([y_ref, y_rx])
+    margin = 0.1 * max(np.ptp(all_x), np.ptp(all_y))
+    plt.xlim(min(all_x)-margin, max(all_x)+margin)
+    plt.ylim(min(all_y)-margin, max(all_y)+margin)
+
+    plt.xlabel('In-phase (I)')
+    plt.ylabel('Quadrature (Q)')
+    plt.title(title)
+    plt.grid(False)
+    plt.gca().set_aspect('equal', adjustable='box')
+    plt.legend()
+    plt.show()
+
+plot_constellations("constellation_ref.dat", "constellation.dat", title="Constellation QAM")
+

QAM/OLD/debug2.py

@@ -0,0 +1,68 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from watchdog.observers import Observer
+from watchdog.events import FileSystemEventHandler
+import time
+import os
+
+needs_update = True  # flag global
+
+def plot_constellations(ref_file, rx_file, title="Constellation Comparison", save_path=None):
+    ref_data = np.atleast_2d(np.loadtxt(ref_file))
+    rx_data = np.atleast_2d(np.loadtxt(rx_file))
+    
+    x_ref, y_ref = ref_data[:,0], ref_data[:,1]
+    x_rx, y_rx = rx_data[:,0], rx_data[:,1]
+
+    plt.clf()  # efface la figure précédente
+    plt.scatter(x_ref, y_ref, color='dodgerblue', s=80, marker='o', edgecolors='k', label='Référence')
+    plt.scatter(x_rx, y_rx, color='tomato', s=80, marker='x', alpha=0.6, label='Reçu')
+
+    all_x = np.concatenate([x_ref, x_rx])
+    all_y = np.concatenate([y_ref, y_rx])
+    margin = 0.15 * max(np.ptp(all_x), np.ptp(all_y))
+    plt.xlim(min(all_x)-margin, max(all_x)+margin)
+    plt.ylim(min(all_y)-margin, max(all_y)+margin)
+
+    plt.grid(True, which='both', linestyle='--', linewidth=0.5, alpha=0.7)
+    plt.axhline(0, color='black', linewidth=1)
+    plt.axvline(0, color='black', linewidth=1)
+    plt.xlabel('In-phase (I)', fontsize=12)
+    plt.ylabel('Quadrature (Q)', fontsize=12)
+    plt.title(title, fontsize=14, fontweight='bold')
+    plt.gca().set_aspect('equal', adjustable='box')
+    plt.legend()
+    if save_path:
+        plt.savefig(save_path, dpi=300, bbox_inches='tight')
+    plt.pause(0.1)
+
+class FileChangeHandler(FileSystemEventHandler):
+    def on_modified(self, event):
+        global needs_update
+        if event.src_path.endswith("constellation_ref.dat") or event.src_path.endswith("constellation.dat"):
+            print(f"{event.src_path} modifié")
+            needs_update = True  # on ne fait que signaler
+
+if __name__ == "__main__":
+    ref_file = "constellation_ref.dat"
+    rx_file = "constellation.dat"
+
+    plt.ion()
+    plt.figure(figsize=(7,7))
+    plot_constellations(ref_file, rx_file)
+
+    event_handler = FileChangeHandler()
+    observer = Observer()
+    observer.schedule(event_handler, path=os.path.dirname(os.path.abspath(ref_file)) or '.', recursive=False)
+    observer.start()
+
+    try:
+        while True:
+            if needs_update:
+                plot_constellations(ref_file, rx_file)
+                needs_update = False
+            time.sleep(0.2)  # boucle principale
+    except KeyboardInterrupt:
+        observer.stop()
+    observer.join()
+

QAM/OLD/qam.c

@@ -0,0 +1,240 @@
+#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
+    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
+};
+typedef struct qam_system_s qam_system;
+
+// Initialisation de la constellation (double tableau de taille sqrt(M)), 
+// ToDo : changer à un tableau à 1 dimension pour éviter de calculer sqrt(M)
+void init_constellation (qam_system* qam) {
+    int sm = (int)sqrt(qam->M);
+    qam->constellation = (double complex**)malloc(sizeof(double complex*) * sm);
+    
+    for (int i = 0; i < sm; i++) {
+       qam->constellation[i] = (double complex*)malloc(sizeof(double complex) * sm);
+    }
+
+    double norm_factor = sqrt((double)(qam->M - 1) / 3.0); // Pour puissance unitaire
+
+    for (int i = 0; i < sm; i++) {
+        double complex ip = -(sm - 1) + 2 * i;
+        for (int j = 0; j < sm; j++) {
+            double complex qp = -(sm - 1) + 2 * j;
+            qam->constellation[i][j] = (ip + I * qp) / norm_factor; 
+        }
+    }
+}
+
+// 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) {
+    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;
+    int sm = sqrt(qam->M);
+    for (int k = 0; k < nb_symbols; k++) {
+        int id = 0;
+        for (int b = 0 ; b < qam->k; b++) {
+            id = id * 2 + bits[k * qam->k + b];
+        }
+        int i = id / sm; 
+        int j = id % sm;
+        symbols[k] = qam->constellation[i][j];
+    }
+}
+
+// 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)
+        int sm = (int)sqrt(qam->M);
+        double min_d = INFINITY;
+        int i_cl, j_cl = 0;
+        for (int i = 0; i < sm; i++) {
+            for (int j = 0; j < sm; j++) {
+                double d = cabs(r - qam->constellation[i][j]);
+                if (d < min_d) {
+                    min_d = d;
+                    i_cl = i;
+                    j_cl = j;
+                }
+            }
+        }
+
+        // index du symbole (id) : même mappage que dans bits_to_symbols()
+        int id = i_cl * sm + j_cl;
+
+        for (int b = 0; b < qam->k; b++) {
+            bits_hat[k * qam->k + b] = (id >> (qam->k - 1 - b)) & 1;
+        }
+    }
+}
+
+// Libération de la mémoire
+void free_constellation(qam_system* qam) {
+    int sm = (int)sqrt(qam->M);
+    for (int i = 0; i < sm; i++)
+        free(qam->constellation[i]);
+    free(qam->constellation);
+}
+
+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;
+}
+
+int main () {
+    qam_system qam;
+    qam.M = 16;
+    qam.k = (int)log2((double)(qam.M));
+    qam.Fs = 44100;
+    qam.Ts = 0.0003;
+    qam.N = (int)qam.Fs * qam.Ts;
+    qam.Fc = 2000;
+    init_constellation(&qam);
+
+    //int nb_bits = 1000;
+    //int nb_symbols = nb_bits / qam.k;
+
+    //uint8_t* input_bits = malloc(nb_bits * sizeof(uint8_t));
+    //for (int i = 0; i < nb_bits; i++) {
+    //    input_bits[i] = rand() % 2;
+    //}
+    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, ";  
+    int nb_chars = strlen(texte);
+    int nb_bits = nb_chars * 8;
+    int nb_symbols = (nb_bits + qam.k - 1) / qam.k;
+
+    // Conversion du texte en bits
+    uint8_t* input_bits = malloc(nb_bits * sizeof(uint8_t));
+    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;
+        }
+    }
+    
+    // Conversion en symboles
+    double complex* symbols = malloc(sizeof(double complex) * nb_symbols);
+    bits_to_symbols(&qam, input_bits, nb_bits, symbols);
+
+    // Modulation
+    int total_samples = qam.N * nb_symbols;
+    double complex* s = malloc(sizeof(double complex) * total_samples);
+    modulate(&qam, symbols, nb_symbols, s);
+
+    // Ajout du bruit
+    double signal_power = (2.0/3.0)*(qam.M-1); // puissance moyenne
+    double snr_dB = 5; // SNR en dB
+    double snr_lin = pow(10.0, snr_dB / 10.0);
+    double sigma = sqrt(signal_power / snr_lin);
+    add_noise(s, total_samples, sigma);
+
+
+    FILE *fp_ref = fopen("constellation_ref.dat", "w");
+    int sm = (int)sqrt(qam.M);
+    for (int i = 0; i < sm; i++) {
+        for (int j = 0; j < sm; j++) {
+            fprintf(fp_ref, "% .8f % .8f\n", creal(qam.constellation[i][j]), cimag(qam.constellation[i][j]));
+        }
+    }
+    fclose(fp_ref);
+    FILE *fp_constel = fopen("constellation.dat", "w");
+
+    double phase_offset = M_PI / 6.0;  // 30 degrés
+    for (int i = 0; i < total_samples; i++) {
+        s[i] *= cexp(I * phase_offset);
+    }
+    double freq_offset = 0; // Hz de décalage
+    for (int i = 0; i < total_samples; i++) {
+        double t = (double)i / qam.Fs;
+        s[i] *= cexp(I * 2 * M_PI * freq_offset * t);
+    }
+    //int offset_samples = (int)(0.3 * qam.N); // décalage de 30% d’un symbole
+    //memmove(s + offset_samples, s, (total_samples - offset_samples) * sizeof(double complex));
+
+
+    // Démodulation
+    uint8_t* output_bits = (uint8_t*)malloc(nb_bits * sizeof(uint8_t));
+    demodulate(&qam, s, nb_symbols, output_bits, fp_constel);
+
+    fclose(fp_constel);
+
+    // Reconstruction du texte
+    char* texte_recup = malloc(nb_chars + 1);
+    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';
+    printf("Texte original : %s\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 blind QAM: %.4f\n", ber * 100);
+
+    // Libération mémoire
+    free(input_bits);
+    free(output_bits);
+    free(symbols);
+    free(s);
+    free_constellation(&qam);
+
+    return 0;
+}

QAM/debug.py

@@ -1,32 +1,81 @@
+#!/usr/bin/env python3
+import pyqtgraph as pg
+from pyqtgraph.Qt import QtWidgets, QtCore
 import numpy as np
-import matplotlib.pyplot as plt
+import sys, os
 
-def plot_constellations(ref_file, rx_file, title="Constellation comparison"):
+REF_FILE = "constellation_ref.dat"
-    # Charger et forcer 2D
+RX_FILE  = "constellation.dat"
-    ref_data = np.atleast_2d(np.loadtxt(ref_file))
-    rx_data = np.atleast_2d(np.loadtxt(rx_file))
-    
-    x_ref, y_ref = ref_data[:,0], ref_data[:,1]
-    x_rx, y_rx = rx_data[:,0], rx_data[:,1]
 
-    plt.figure(figsize=(6,6))
+def load_points(filename):
-    plt.scatter(x_ref, y_ref, color='blue', s=50, marker='o', label='Référence')
+    if os.path.exists(filename):
-    plt.scatter(x_rx, y_rx, color='red', s=50, marker='x', label='Reçu')
+        return np.loadtxt(filename)
+    else:
+        return np.zeros((0,2))
 
-    # Ajustement automatique des limites
+# -------------------- Application --------------------
-    all_x = np.concatenate([x_ref, x_rx])
+app = QtWidgets.QApplication(sys.argv)
-    all_y = np.concatenate([y_ref, y_rx])
+win = pg.GraphicsLayoutWidget(show=True, title="Constellation Viewer")
-    margin = 0.1 * max(np.ptp(all_x), np.ptp(all_y))
+win.resize(900, 900)
-    plt.xlim(min(all_x)-margin, max(all_x)+margin)
+win.setBackground('#0a0a0a')  # fond noir profond
-    plt.ylim(min(all_y)-margin, max(all_y)+margin)
 
-    plt.xlabel('In-phase (I)')
+plot = win.addPlot()
-    plt.ylabel('Quadrature (Q)')
+plot.setAspectLocked(True)  # axes égaux
-    plt.title(title)
-    plt.grid(False)
-    plt.gca().set_aspect('equal', adjustable='box')
-    plt.legend()
-    plt.show()
 
-plot_constellations("constellation_ref.dat", "constellation.dat", title="Constellation QAM")
+# -------------------- 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)
+
+# -------------------- Chargement initial des points --------------------
+ref_data = load_points(REF_FILE)
+rx_data  = load_points(RX_FILE)
+
+# Points sans contour
+ref_plot = plot.plot(ref_data[:,0], ref_data[:,1],
+                     pen=None,
+                     symbol='o',
+                     symbolSize=10,             # cercles un peu plus gros
+                     symbolBrush=pg.mkBrush('#1E90FF'),  # bleu vif
+                     symbolPen=None,
+                     name='Référence')
+
+rx_plot = plot.plot(rx_data[:,0], rx_data[:,1],
+                    pen=None,
+                    symbol='x',
+                    symbolSize=6,              # croix plus petites
+                    symbolBrush=pg.mkBrush('#FF4500'),  # orange vif
+                    symbolPen=None,
+                    name='Reçu')
+
+# Légende stylée
+legend = plot.addLegend()
+for item in legend.items:
+    item[1].setPen(pg.mkPen('#FFF', width=2))
+
+# -------------------- Timer pour mise à jour dynamique --------------------
+def update():
+    ref_data = load_points(REF_FILE)
+    rx_data  = load_points(RX_FILE)
+    ref_plot.setData(ref_data[:,0], ref_data[:,1])
+    rx_plot.setData(rx_data[:,0], rx_data[:,1])
+
+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_())
 

QAM/qam.c

@@ -99,7 +99,8 @@ void demodulate(qam_system* qam, double complex* s, int nb_symbols, uint8_t* bit
         // Distance euclidien de Ir et Qr pour avoir le point le plus proche de la constellation (lent)
         int sm = (int)sqrt(qam->M);
         double min_d = INFINITY;
-        int i_cl, j_cl = 0;
+        int i_cl = 0;
+        int j_cl = 0;
         for (int i = 0; i < sm; i++) {
             for (int j = 0; j < sm; j++) {
                 double d = cabs(r - qam->constellation[i][j]);
@@ -120,6 +121,50 @@ void demodulate(qam_system* qam, double complex* s, int nb_symbols, uint8_t* bit
     }
 }
 
+// PLL pour corriger le déphasage
+void pll_qam_symbol(qam_system* qam, double complex* symbols_rx, double complex* r_corr, int nb_symbols, double Kp, double Ki) {
+
+    double phase_est = 0.0;
+    double integrator = 0.0;
+
+    int sm = (int)sqrt(qam->M);
+    int N = qam->N;
+
+    for (int k = 0; k < nb_symbols; k++) {
+        double complex r_symbol = 0;
+        for (int n = 0; n < N; n++) {
+            int idx = k * N + n;
+            r_symbol += symbols_rx[idx] * cexp(2.0 * -I * M_PI * qam->Fc * ((double)idx / qam->Fs));
+        }
+        r_symbol /= N;
+
+        r_symbol *= cexp(-I * phase_est);
+
+        double min_d = INFINITY;
+        double complex closest = 0;
+        for (int i = 0; i < sm; i++) {
+            for (int j = 0; j < sm; j++) {
+                double d = cabs(r_symbol - qam->constellation[i][j]);
+                if (d < min_d) {
+                    min_d = d;
+                    closest = qam->constellation[i][j];
+                }
+            }
+        }
+
+        double error = carg(r_symbol * conj(closest));
+
+        integrator += Ki * error;
+        phase_est += Kp * error + integrator;
+
+        for (int n = 0; n < N; n++) {
+            int idx = k * N + n;
+            r_corr[idx] = symbols_rx[idx] * cexp(-I * phase_est);
+        }
+    }
+}
+
+
 // Libération de la mémoire
 void free_constellation(qam_system* qam) {
     int sm = (int)sqrt(qam->M);
@@ -136,13 +181,53 @@ double compare_bits(uint8_t* bits1, uint8_t* bits2, int nb_bits) {
     return (double)errors / nb_bits;
 }
 
+// Minimise le BER (si la pll s'est lockée de maniere déphasée de k*pi/2) 
+void demodulate2(qam_system* qam, double complex* r_corr, int nb_symbols, uint8_t* input_bits, uint8_t* output_bits, FILE* fp_constel) {
+
+    int nb_bits = nb_symbols * qam->k;
+    double best_ber = INFINITY;
+    double best_angle = 0.0;
+    uint8_t* temp_bits = (uint8_t*)malloc(nb_bits * sizeof(uint8_t));
+
+    for (int r = 0; r < 4; r++) {
+        double angle = r * M_PI/2;
+
+        double complex* rotated = (double complex*)malloc(sizeof(double complex) * nb_symbols * qam->N);
+        for (int i = 0; i < nb_symbols * qam->N; i++) {
+            rotated[i] = r_corr[i] * cexp(I * angle);
+        }
+
+        demodulate(qam, rotated, nb_symbols, temp_bits, fp_constel);
+
+        double ber = compare_bits(input_bits, temp_bits, nb_bits);
+
+        if (ber < best_ber) {
+            best_ber = ber;
+            best_angle = angle;
+        }
+
+        free(rotated);
+    }
+
+    for (int i = 0; i < nb_symbols * qam->N; i++) {
+        r_corr[i] *= cexp(I * best_angle);
+    }
+
+    demodulate(qam, r_corr, nb_symbols, output_bits, fp_constel);
+
+    free(temp_bits);
+}
+
+
 int main () {
     qam_system qam;
     qam.M = 16;
     qam.k = (int)log2((double)(qam.M));
     qam.Fs = 44100;
-    qam.Ts = 0.0003;
+    //qam.Ts = 0.0003;
-    qam.N = (int)qam.Fs * qam.Ts;
+    //qam.N = (int)qam.Fs * qam.Ts;
+    qam.Ts = 0.01;
+    qam.N = (int)(qam.Fs * qam.Ts);
     qam.Fc = 2000;
     init_constellation(&qam);
 
@@ -153,7 +238,7 @@ int main () {
     //for (int i = 0; i < nb_bits; i++) {
     //    input_bits[i] = rand() % 2;
     //}
-    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, ";  
+    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";  
     int nb_chars = strlen(texte);
     int nb_bits = nb_chars * 8;
     int nb_symbols = (nb_bits + qam.k - 1) / qam.k;
@@ -180,7 +265,7 @@ int main () {
     double snr_dB = 5; // SNR en dB
     double snr_lin = pow(10.0, snr_dB / 10.0);
     double sigma = sqrt(signal_power / snr_lin);
-    add_noise(s, total_samples, sigma);
+    add_noise(s, total_samples, 20);
 
 
     FILE *fp_ref = fopen("constellation_ref.dat", "w");
@@ -193,22 +278,32 @@ int main () {
     fclose(fp_ref);
     FILE *fp_constel = fopen("constellation.dat", "w");
 
-    double phase_offset = M_PI / 6.0;  // 30 degrés
+    // Ajout de dephasage
-    for (int i = 0; i < total_samples; i++) {
+    //double phase_offset = M_PI / 6.0;  // 30 degrés
-        s[i] *= cexp(I * phase_offset);
+    //for (int i = 0; i < total_samples; i++) {
-    }
+    //    s[i] *= cexp(I * phase_offset);
-    double freq_offset = 0; // Hz de décalage
+    //}
+    
+    // AJout de decalage de fréquence
+    double freq_offset = 1; // Hz de décalage
     for (int i = 0; i < total_samples; i++) {
         double t = (double)i / qam.Fs;
         s[i] *= cexp(I * 2 * M_PI * freq_offset * t);
     }
+    
+    // Ajout de decalage entre les symbole
     //int offset_samples = (int)(0.3 * qam.N); // décalage de 30% d’un symbole
     //memmove(s + offset_samples, s, (total_samples - offset_samples) * sizeof(double complex));
 
+    double complex* r_corr = malloc(sizeof(double complex) * total_samples);
+    double Kp = 0.25;
+    double Ki = 0.02;
+    pll_qam_symbol(&qam, s, r_corr, nb_symbols, Kp, Ki);
 
     // Démodulation
     uint8_t* output_bits = (uint8_t*)malloc(nb_bits * sizeof(uint8_t));
-    demodulate(&qam, s, nb_symbols, output_bits, fp_constel);
+    demodulate2(&qam, r_corr, nb_symbols, input_bits, output_bits, fp_constel);
+    //demodulate(&qam, r_corr, nb_symbols, output_bits, fp_constel);
 
     fclose(fp_constel);
 
@@ -222,7 +317,7 @@ int main () {
         texte_recup[i] = c;
     }
     texte_recup[nb_chars] = '\0';
-    printf("Texte original : %s\n", texte);
+    printf("Texte original : %s\n\n", texte);
     printf("Texte demodulé : %s\n", texte_recup);
 
     // Calcul du BER
@@ -233,6 +328,7 @@ int main () {
     free(input_bits);
     free(output_bits);
     free(symbols);
+    free(r_corr);
     free(s);
     free_constellation(&qam);
 

QAM/todo.md

@@ -0,0 +1,157 @@
+
+# 🧭 Ordre d’implémentation d’une chaîne de réception QAM (réelle)
+
+## **Étape 0 — Environnement et base de test**
+
+🎯 Objectif : avoir une base stable de simulation avant toute boucle adaptative.
+
+**À faire :**
+
+* Lire un **flux IQ** (fichier ou SDR)
+* Implémenter un **filtre RRC de réception**
+* Visualiser constellations, spectre, symboles
+
+🧠 But : voir des symboles flous mais reconnaissables — aucun algorithme adaptatif encore.
+
+---
+
+## **Étape 1 — Correction de fréquence (CFO / Carrier Recovery)**
+
+🎯 Objectif : supprimer le décalage de fréquence de la porteuse avant le timing recovery.
+
+**Pourquoi en premier ?**
+
+* Si ton signal tourne dans le plan complexe, **Mueller & Müller échouera** complètement.
+* Il faut un signal “quasistationnaire” avant de chercher le bon instant d’échantillonnage.
+
+**Méthodes à implémenter :**
+
+* Estimation grossière de CFO (par corrélation / FFT)
+* Boucle de Costas ou PLL de phase
+
+📘 Résultat attendu : constellation QAM fixe mais “floue” (problème de timing restant).
+
+---
+
+## **Étape 2 — Synchronisation temporelle (Timing Recovery : M&M)**
+
+🎯 Objectif : trouver l’instant exact d’échantillonnage par symbole.
+
+**Tu connais déjà :**
+
+* Interpolateur fractionnaire
+* Détecteur d’erreur M&M
+* Boucle PI
+
+💡 Astuce : commence avec un signal *parfaitement corrigé en fréquence* avant d’activer la boucle timing.
+
+📘 Résultat attendu : points de constellation bien centrés, toujours un peu brouillés (canal non corrigé).
+
+---
+
+## **Étape 3 — Égalisation adaptative (FFE / CMA / DD-LMS)**
+
+🎯 Objectif : supprimer l’ISI et compenser la distorsion de canal.
+
+**Méthodes typiques :**
+
+* **CMA** (Constant Modulus Algorithm) pour pré-verrouillage (aveugle)
+* Puis **DD-LMS** (Decision Directed) une fois la décision fiable
+
+💡 L’égaliseur doit venir **après** le timing (sinon le signal est mal échantillonné).
+
+📘 Résultat attendu : constellation nette, points regroupés correctement autour des symboles 16-QAM.
+
+---
+
+## **Étape 4 — Décision + Mapping**
+
+🎯 Objectif : convertir les symboles QAM en bits.
+
+**À faire :**
+
+* Implémenter la décision dure (±1, ±3)
+* Mapping / demapping Gray
+* Vérifier BER par rapport à trame connue
+
+💡 Teste d’abord sans bruit pour valider le mapping bit ↔ symbole.
+
+---
+
+## **Étape 5 — Boucles de phase fines (Costas loop fine)**
+
+🎯 Objectif : corriger la phase résiduelle après timing et égalisation.
+
+* Souvent intégrée dans la boucle M&M ou séparée (PLL de phase fine)
+* Sert à verrouiller la dernière rotation du plan IQ
+
+---
+
+## **Étape 6 — Correction d’erreurs / décodage (FEC)**
+
+🎯 Objectif : terminer la chaîne par le décodage des bits.
+
+* LDPC, Viterbi, Turbo selon ton système
+* C’est la couche “bitstream”, plus logique que DSP
+
+---
+
+## **Étape 7 — Optimisation et intégration**
+
+🎯 Objectif : passer du prototype à la version embarquée.
+
+* Conversion float → fixe (Q-format)
+* Pipeline temps réel (DMA, buffers circulaires)
+* Profiling CPU / mémoire
+* Test sur matériel SDR, puis en RF réelle
+
+---
+
+# ⚙️ En résumé — Ordre d’implémentation
+
+| Étape | Bloc                                 | Type            | Pourquoi cet ordre          |
+| ----- | ------------------------------------ | --------------- | --------------------------- |
+| 0     | RRC + acquisition IQ                 | statique        | Base stable et visualisable |
+| 1     | **Correction de fréquence (CFO)**    | boucle 1        | sinon M&M échoue            |
+| 2     | **Synchronisation temporelle (M&M)** | boucle 2        | aligner les symboles        |
+| 3     | **Égalisation adaptative**           | boucle 3        | corriger canal              |
+| 4     | **Décision + mapping bits**          | logique         | extraire données            |
+| 5     | **Boucle de phase fine (Costas)**    | ajustement      | phase finale                |
+| 6     | **Décodage FEC**                     | post-traitement | fiabiliser le bitstream     |
+| 7     | **Optimisation C/FPGA**              | système         | rendre temps réel           |
+
+---
+
+# 🧠 Ordre de test conseillé
+
+1. Simule tout en **float** dans Python/MATLAB
+2. Valide chaque bloc **indépendamment**
+3. Assemble et teste avec bruit / décalage
+4. **Ensuite seulement**, porte en C (ou sur DSP/FPGA)
+
+---
+
+---
+
+| # | Bloc / Étape                            | Objectif principal                              | Priorité   | Ce qu’il faut coder / comprendre                | Test / Validation                                |
+| - | --------------------------------------- | ----------------------------------------------- | ---------- | ----------------------------------------------- | ------------------------------------------------ |
+| 0 | **Acquisition & Filtrage RRC**          | Lire le signal IQ et filtrer pour limiter l’ISI | Très haute | FIR RRC, buffer IQ                              | Visualiser constellation, spectre                |
+| 1 | **Correction de fréquence (CFO)**       | Supprimer offset de fréquence de la porteuse    | Très haute | PLL / Costas loop, corrélation, FFT             | Constellation immobile, pas de rotation          |
+| 2 | **Synchronisation temporelle (M&M)**    | Aligner échantillons sur symboles               | Très haute | Interpolateur fractionnaire, TED M&M, boucle PI | Points centrés, vérification de tau_hat          |
+| 3 | **Égalisation adaptative**              | Supprimer ISI et compenser canal                | Haute      | FFE ou DFE, algorithme CMA / DD-LMS             | Constellation nette, erreur moyenne faible       |
+| 4 | **Décision symbolique et mapping**      | Convertir symbole → bits                        | Haute      | Hard decision QAM, Gray mapping                 | Vérifier BER avec trame connue                   |
+| 5 | **Boucle de phase fine (Costas / PLL)** | Corriger la phase résiduelle                    | Moyenne    | PLL numérique, phase fine                       | Points de constellation fixes, phase verrouillée |
+| 6 | **Décodage FEC**                        | Extraire flux de bits fiable                    | Moyenne    | LDPC / Viterbi / Turbo                          | Comparer bits reçus / transmis, BER              |
+| 7 | **Optimisation & passage temps réel**   | Adapter pour C / DSP / FPGA                     | Moyenne    | Point fixe, buffers circulaires, pipeline       | Profil CPU / mémoire, latence, test en SDR réel  |
+
+---
+
+### 💡 Notes pratiques :
+
+* **Test bloc par bloc** avant d’intégrer la chaîne complète
+* Toujours **simuler en float** avant passage en C ou point fixe
+* Chaque boucle (CFO, Timing, Phase) doit être **réglée indépendamment** pour éviter l’instabilité
+* Commencer avec **trames simples** avant bruit réel, puis ajouter AWGN / jitter / offsets
+
+---
+