DSP/QAM/qam.c

#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <string.h>

#define A 1 

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); 
        }
    }
}

// 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++) {
            s[k * qam->N + n] = iq * cexp(2 * I * M_PI * qam->Fc * ((double)n / qam->Fs)); 
        }
    }
}

// Demodulation QAM 
void demodulate(qam_system* qam, double complex* s, int nb_symbols, uint8_t* bits_hat) {
    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)n / qam->Fs));
        }
        r /= qam->N;

        // 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;
}

void symbol_timing_recovery(double complex* r, int r_len, double complex* preamble, int preamble_len, int N_max, int* offset_est, int* N_est) {
    double max_corr = 0;
    int best_offset = 0;
    int best_N = 1;

    for(int N_try = 1; N_try <= N_max; N_try++) {
        for(int off = 0; off <= r_len - preamble_len * N_try; off++) {
            double complex corr = 0;
            for(int k = 0; k < preamble_len; k++) {
                for(int n = 0; n < N_try; n++) {
                    corr += r[off + k * N_try + n] * conj(preamble[k]);
                }
            }
            double mag = cabs(corr);
            if(mag > max_corr) {
                max_corr = mag;
                best_offset = off;
                best_N = N_try;
            }
        }
    }

    *offset_est = best_offset;
    *N_est = best_N;
}
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 = "Bonjour tout le monde salut !!";  
    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);

    // Démodulation
    //uint8_t* output_bits = malloc(nb_bits * sizeof(uint8_t));
    //demodulate(&qam, s, nb_symbols, output_bits);

    // Calcul du taux d'erreur
    //double ber = compare_bits(input_bits, output_bits, nb_bits);
    //printf("Taux d'erreur : %.4f\n", ber *  100);

    // Blind QAM
    int preamble_symbols = 32;
    uint8_t* preamble_bits = malloc(preamble_symbols * qam.k * sizeof(uint8_t));
    // Genere un préambule aléatoire
    for(int i = 0; i < preamble_symbols * qam.k; i++) {
        preamble_bits[i] = rand() % 2;
    }
    double complex* preamble_symbols_complex = (double complex*)malloc(sizeof(double complex) * preamble_symbols);
    bits_to_symbols(&qam, preamble_bits, preamble_symbols * qam.k, preamble_symbols_complex);
    // Ajout du préambule dans les données
    int total_symbols = preamble_symbols + nb_symbols;
    double complex* all_symbols = malloc(sizeof(double complex) * total_symbols);
    // Copier le préambule
    for(int i = 0; i < preamble_symbols; i++) {
        all_symbols[i] = preamble_symbols_complex[i];
    }
    // Copier les symboles réels
    for(int i = 0; i < nb_symbols; i++) {
        all_symbols[i + preamble_symbols] = symbols[i];
    }
    int total_samples_with_preamble = total_symbols * qam.N;
    double complex* s_with_preamble = malloc(sizeof(double complex) * total_samples_with_preamble);

    modulate(&qam, all_symbols, total_symbols, s_with_preamble);


    // Symbol Timing Recovery blind 
    int offset_est = 0;
    int N_est = 0;
    symbol_timing_recovery(s_with_preamble, total_samples_with_preamble, preamble_symbols_complex, preamble_symbols, qam.N*2, &offset_est, &N_est);

    printf("Symbol Timing Recovery : offset=%d, N_est=%d (%d)\n", offset_est, N_est, qam.N);

    /*
    // Calcul du BER
    double ber = compare_bits(input_bits, output_bits, nb_bits);
    printf("Taux d'erreur blind QAM: %.4f\n", ber * 100);

    // 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);
    */

    // Libération mémoire
    free(input_bits);
    //free(output_bits);
    free(symbols);
    free(s);
    free_constellation(&qam);

    return 0;
}