DSP/QAM/save/2/qam.c

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

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

        /*
        // Ancienne methode (non gray)
        // 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;
        }
        */ 

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

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 = (int)sqrt(qam->M);

    // 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 >= sm || j_gray < 0 || j_gray >= sm) {
            printf("bits_to_symbols : IOOR i_gray = %d j_gray = %d sm = %d \n", i_gray, j_gray, sm);
            exit(1);
        }
        symbols[sym] = qam->constellation[i_gray][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) {
    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]));
        }
    }
}

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) {
    int sm = (int)sqrt(qam->M);
    for (int i = 0; i < sm; i++)
        free(qam->constellation[i]);
    free(qam->constellation);
}

int main () {
    // Initialisation du system qam
    qam_system qam;
    qam.M = 16;
    qam.k = (int)log2((double)(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";  
    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);

    // 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
    add_noise(s, total_samples, 0);

    FILE *fp_ref = fopen("constellation_ref.dat", "w");
    fill_constellation_data(&qam, fp_ref);
    fclose(fp_ref);

    // Ajout de dephasage
    //add_dephasage(s, M_PI / 6.0, total_samples);
    
    // AJout de decalage de fréquence
    //add_freq(&qam, s, 1, total_samples); 

    // Démodulation
    FILE *fp_constel = fopen("constellation.dat", "w");
    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);
    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 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;
}