Little changes and debug
This commit is contained in:
@ -5,11 +5,9 @@
|
|||||||
#include <complex.h>
|
#include <complex.h>
|
||||||
#include "../wav/wav.h"
|
#include "../wav/wav.h"
|
||||||
#include "../files/files.h"
|
#include "../files/files.h"
|
||||||
#include "../plot/plot_constellation.h"
|
|
||||||
#include <string.h>
|
#include <string.h>
|
||||||
#include <SDL2/SDL.h>
|
|
||||||
|
|
||||||
#define A 10000
|
#define A 1
|
||||||
|
|
||||||
struct qam_system_s {
|
struct qam_system_s {
|
||||||
int M; // Nombre de symboles M-QAM
|
int M; // Nombre de symboles M-QAM
|
||||||
@ -38,7 +36,7 @@ void init_constellation (qam_system* qam) {
|
|||||||
double complex ip = -(sm - 1) + 2 * i;
|
double complex ip = -(sm - 1) + 2 * i;
|
||||||
for (int j = 0; j < sm; j++) {
|
for (int j = 0; j < sm; j++) {
|
||||||
double complex qp = -(sm - 1) + 2 * j;
|
double complex qp = -(sm - 1) + 2 * j;
|
||||||
qam->constellation[i][j] = A * (ip + I * qp) / norm_factor;
|
qam->constellation[i][j] = (ip + I * qp);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
@ -94,10 +92,9 @@ void demodulate(qam_system* qam, double complex* s, int nb_symbols, uint8_t* bit
|
|||||||
}
|
}
|
||||||
r /= qam->N;
|
r /= qam->N;
|
||||||
|
|
||||||
// Distance euclidien de Ir et Qr pour avoir le point le plus proche de la constellation
|
// Distance euclidien de Ir et Qr pour avoir le point le plus proche de la constellation (lent)
|
||||||
|
|
||||||
int sm = (int)sqrt(qam->M);
|
int sm = (int)sqrt(qam->M);
|
||||||
/* TEMPS INFINI
|
|
||||||
double min_d = INFINITY;
|
double min_d = INFINITY;
|
||||||
int i_cl, j_cl = 0;
|
int i_cl, j_cl = 0;
|
||||||
for (int i = 0; i < sm; i++) {
|
for (int i = 0; i < sm; i++) {
|
||||||
@ -110,39 +107,16 @@ void demodulate(qam_system* qam, double complex* s, int nb_symbols, uint8_t* bit
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
*/
|
|
||||||
double norm_factor = sqrt((double)(qam->M - 1) / 3.0);
|
|
||||||
double Ir = creal(r) * norm_factor / A;
|
|
||||||
double Qr = cimag(r) * norm_factor / A;
|
|
||||||
|
|
||||||
int i = (int)round((Ir + (sm - 1)) / 2.0);
|
// index du symbole (id) : même mappage que dans bits_to_symbols()
|
||||||
int j = (int)round((Qr + (sm - 1)) / 2.0);
|
int id = i_cl * sm + j_cl;
|
||||||
|
|
||||||
i = (i < 0) ? 0 : ((i >= sm) ? sm - 1 : i);
|
|
||||||
j = (j < 0) ? 0 : ((j >= sm) ? sm - 1 : j);
|
|
||||||
|
|
||||||
int id = i * sm + j;
|
|
||||||
|
|
||||||
//int id = i_cl * sm + j_cl;
|
|
||||||
for (int b = 0; b < qam->k; b++) {
|
for (int b = 0; b < qam->k; b++) {
|
||||||
bits_hat[k * qam->k + (qam->k - 1 - b)] = (id >> b) & 1;
|
bits_hat[k * qam->k + b] = (id >> (qam->k - 1 - b)) & 1;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
double complex* demodulate_points(qam_system* qam, double complex* s, int nb_symbols) {
|
|
||||||
double complex* points = malloc(sizeof(double complex) * nb_symbols);
|
|
||||||
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;
|
|
||||||
points[k] = r;
|
|
||||||
}
|
|
||||||
return points;
|
|
||||||
}
|
|
||||||
|
|
||||||
// Libération de la mémoire
|
// Libération de la mémoire
|
||||||
void free_constellation(qam_system* qam) {
|
void free_constellation(qam_system* qam) {
|
||||||
int sm = (int)sqrt(qam->M);
|
int sm = (int)sqrt(qam->M);
|
||||||
@ -151,6 +125,36 @@ void free_constellation(qam_system* qam) {
|
|||||||
free(qam->constellation);
|
free(qam->constellation);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
// Compare deux tableaux de bits (0/1) et retourne le pourcentage de fiabilité.
|
||||||
|
double compare_bits(const uint8_t *in_bits, size_t nb_bits_in, const uint8_t *out_bits, size_t nb_bits_out, size_t *erreurs) {
|
||||||
|
if (!in_bits || !out_bits) {
|
||||||
|
if (erreurs) *erreurs = (nb_bits_in < nb_bits_out) ? nb_bits_out : nb_bits_in;
|
||||||
|
return 0.0;
|
||||||
|
}
|
||||||
|
|
||||||
|
size_t n_min = (nb_bits_in < nb_bits_out) ? nb_bits_in : nb_bits_out;
|
||||||
|
size_t n_max = (nb_bits_in > nb_bits_out) ? nb_bits_in : nb_bits_out;
|
||||||
|
|
||||||
|
size_t err = 0;
|
||||||
|
for (size_t i = 0; i < n_min; ++i) {
|
||||||
|
if ((in_bits[i] & 1) != (out_bits[i] & 1)) err++;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (n_max != n_min) {
|
||||||
|
err += (n_max - n_min);
|
||||||
|
}
|
||||||
|
|
||||||
|
if (erreurs) *erreurs = err;
|
||||||
|
|
||||||
|
double total_compared = (double)n_max;
|
||||||
|
if (total_compared == 0.0) return 0.0;
|
||||||
|
|
||||||
|
double ber = (double)err / total_compared;
|
||||||
|
double reliability_percent = (1.0 - ber) * 100.0;
|
||||||
|
return reliability_percent;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
int main (int argc, char *argv[]) {
|
int main (int argc, char *argv[]) {
|
||||||
if (argc < 2) {
|
if (argc < 2) {
|
||||||
fprintf(stderr, "Utilisation: %s <fichier_entree>\n", argv[0]);
|
fprintf(stderr, "Utilisation: %s <fichier_entree>\n", argv[0]);
|
||||||
@ -158,7 +162,7 @@ int main (int argc, char *argv[]) {
|
|||||||
}
|
}
|
||||||
|
|
||||||
qam_system qam;
|
qam_system qam;
|
||||||
qam.M = 64;
|
qam.M = 16;
|
||||||
qam.k = (int)log2((double)(qam.M));
|
qam.k = (int)log2((double)(qam.M));
|
||||||
qam.Fs = 44100;
|
qam.Fs = 44100;
|
||||||
qam.Ts = 0.0003;
|
qam.Ts = 0.0003;
|
||||||
@ -185,18 +189,18 @@ int main (int argc, char *argv[]) {
|
|||||||
|
|
||||||
// Ajout du bruit
|
// Ajout du bruit
|
||||||
double signal_power = (2.0/3.0)*(qam.M-1); // puissance moyenne avant échelle
|
double signal_power = (2.0/3.0)*(qam.M-1); // puissance moyenne avant échelle
|
||||||
double snr_dB = -27; // Signal to noise ratio
|
double snr_dB = 10; // Signal to noise ratio
|
||||||
double snr_lin = pow(10.0, snr_dB / 10.0);
|
double snr_lin = pow(10.0, snr_dB / 10.0);
|
||||||
double sigma = sqrt(signal_power / snr_lin);
|
double sigma = sqrt(signal_power / snr_lin);
|
||||||
printf("Ajout du bruit... \n puissance du signal : %f\n SNR db : %f\n sigma : %f\n", signal_power, snr_dB, sigma);
|
printf("Ajout du bruit... \n puissance du signal : %f\n SNR db : %f\n sigma : %f\n", signal_power, snr_dB, sigma);
|
||||||
add_noise(s, total_samples, 1875);
|
add_noise(s, total_samples, 0);
|
||||||
|
|
||||||
printf("Demodulation...\n");
|
printf("Demodulation...\n");
|
||||||
|
|
||||||
// Demodulation QAM
|
// Demodulation QAM
|
||||||
bit_array output_bits;
|
bit_array output_bits;
|
||||||
output_bits.nb_bits = input_bits.nb_bits;
|
output_bits.nb_bits = input_bits.nb_bits;
|
||||||
output_bits.bits = (uint8_t*)malloc(output_bits.nb_bits);
|
output_bits.bits = (uint8_t*)malloc(output_bits.nb_bits * sizeof(uint8_t));
|
||||||
demodulate(&qam, s, nb_symbols, output_bits.bits);
|
demodulate(&qam, s, nb_symbols, output_bits.bits);
|
||||||
|
|
||||||
printf("Ecriture...\n");
|
printf("Ecriture...\n");
|
||||||
@ -204,24 +208,25 @@ int main (int argc, char *argv[]) {
|
|||||||
char *output_filename = make_output_filename(input_filename);
|
char *output_filename = make_output_filename(input_filename);
|
||||||
bits_to_file(output_filename, &output_bits);
|
bits_to_file(output_filename, &output_bits);
|
||||||
|
|
||||||
|
|
||||||
// Affichage du signal dans un .wav
|
// Affichage du signal dans un .wav
|
||||||
|
/*
|
||||||
double* si = (double*)malloc(sizeof(double) * total_samples);
|
double* si = (double*)malloc(sizeof(double) * total_samples);
|
||||||
for (int i = 0; i < total_samples; i++) {
|
for (int i = 0; i < total_samples; i++) {
|
||||||
si[i] = cimag(s[i]);
|
si[i] = cimag(s[i]);
|
||||||
}
|
}
|
||||||
write_wav("output.wav", si, total_samples);
|
write_wav("output.wav", si, total_samples);
|
||||||
|
*/
|
||||||
|
|
||||||
// Plot
|
size_t erreurs = 0;
|
||||||
printf("Ploting...");
|
double fiabilite = compare_bits( input_bits.bits, input_bits.nb_bits, output_bits.bits, output_bits.nb_bits, &erreurs);
|
||||||
plot_t plot;
|
|
||||||
plot_init(&plot, 1400, 1400, 0.04);
|
printf("Comparaison :\n");
|
||||||
plot_draw_constellation(&plot, qam.constellation, qam.M);
|
printf(" Bits d'entrée : %zu\n", input_bits.nb_bits);
|
||||||
SDL_Color red = {255, 0, 0, 255};
|
printf(" Bits de sortie: %zu\n", output_bits.nb_bits);
|
||||||
plot_draw_points_animated(&plot, demodulate_points(&qam, s, nb_symbols), nb_symbols, red, 5);
|
printf(" Erreurs : %zu\n", erreurs);
|
||||||
|
printf(" Fiabilité : %.4f %%\n", fiabilite);
|
||||||
|
|
||||||
// Libération mémoire
|
// Libération mémoire
|
||||||
plot_close(&plot);
|
|
||||||
free_bit_array(&input_bits);
|
free_bit_array(&input_bits);
|
||||||
free_bit_array(&output_bits);
|
free_bit_array(&output_bits);
|
||||||
free(symbols);
|
free(symbols);
|
||||||
237
QAM/qamold/qambis.c
Normal file
237
QAM/qamold/qambis.c
Normal file
@ -0,0 +1,237 @@
|
|||||||
|
#include <math.h>
|
||||||
|
#include <stdint.h>
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <complex.h>
|
||||||
|
#include "../wav/wav.h"
|
||||||
|
#include "../files/files.h"
|
||||||
|
#include "../plot/plot_constellation.h"
|
||||||
|
#include <string.h>
|
||||||
|
#include <SDL2/SDL.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
|
||||||
|
|
||||||
|
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;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
/*
|
||||||
|
double norm_factor = sqrt((double)(qam->M - 1) / 3.0);
|
||||||
|
double Ir = creal(r) * norm_factor / A;
|
||||||
|
double Qr = cimag(r) * norm_factor / A;
|
||||||
|
|
||||||
|
int i = (int)round((Ir + (sm - 1)) / 2.0);
|
||||||
|
int j = (int)round((Qr + (sm - 1)) / 2.0);
|
||||||
|
|
||||||
|
i = (i < 0) ? 0 : ((i >= sm) ? sm - 1 : i);
|
||||||
|
j = (j < 0) ? 0 : ((j >= sm) ? sm - 1 : j);
|
||||||
|
|
||||||
|
int id = i * sm + j;
|
||||||
|
*/
|
||||||
|
|
||||||
|
int id = i_cl * sm + j_cl;
|
||||||
|
for (int b = 0; b < qam->k; b++) {
|
||||||
|
bits_hat[k * qam->k + (qam->k - 1 - b)] = (id >> b) & 1;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
double complex* demodulate_points(qam_system* qam, double complex* s, int nb_symbols) {
|
||||||
|
double complex* points = malloc(sizeof(double complex) * nb_symbols);
|
||||||
|
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;
|
||||||
|
double norm_factor = sqrt((double)(qam->M - 1) / 3.0);
|
||||||
|
double Ir = creal(r);
|
||||||
|
double Qr = cimag(r);
|
||||||
|
|
||||||
|
points[k] = Ir + I * Qr;
|
||||||
|
}
|
||||||
|
return points;
|
||||||
|
}
|
||||||
|
|
||||||
|
// 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 (int argc, char *argv[]) {
|
||||||
|
if (argc < 2) {
|
||||||
|
fprintf(stderr, "Utilisation: %s <fichier_entree>\n", argv[0]);
|
||||||
|
return 1;
|
||||||
|
}
|
||||||
|
|
||||||
|
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);
|
||||||
|
|
||||||
|
printf("Lecture du fichier...\n");
|
||||||
|
// Lecture du fichier et conversion en bits
|
||||||
|
const char *input_filename = argv[1];
|
||||||
|
bit_array input_bits = file_to_bits(input_filename);
|
||||||
|
size_t nb_symbols = input_bits.nb_bits / qam.k;
|
||||||
|
|
||||||
|
printf("Mise en forme des symboles...\n");
|
||||||
|
// Mise en forme des symboles
|
||||||
|
double complex *symbols = malloc(sizeof(double complex) * nb_symbols);
|
||||||
|
bits_to_symbols(&qam, input_bits.bits, input_bits.nb_bits, symbols);
|
||||||
|
|
||||||
|
printf("Modulation...\n");
|
||||||
|
// Modulation QAM
|
||||||
|
int total_samples = qam.N * nb_symbols;
|
||||||
|
double complex* s = (double complex*)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 avant échelle
|
||||||
|
double snr_dB = 10; // Signal to noise ratio
|
||||||
|
double snr_lin = pow(10.0, snr_dB / 10.0);
|
||||||
|
double sigma = sqrt(signal_power / snr_lin);
|
||||||
|
printf("Ajout du bruit... \n puissance du signal : %f\n SNR db : %f\n sigma : %f\n", signal_power, snr_dB, sigma);
|
||||||
|
add_noise(s, total_samples, sigma);
|
||||||
|
|
||||||
|
printf("Demodulation...\n");
|
||||||
|
|
||||||
|
// Demodulation QAM
|
||||||
|
bit_array output_bits;
|
||||||
|
output_bits.nb_bits = input_bits.nb_bits;
|
||||||
|
output_bits.bits = (uint8_t*)malloc(output_bits.nb_bits);
|
||||||
|
demodulate(&qam, s, nb_symbols, output_bits.bits);
|
||||||
|
|
||||||
|
printf("Ecriture...\n");
|
||||||
|
// Ecriture du fichier de Demodulation
|
||||||
|
char *output_filename = make_output_filename(input_filename);
|
||||||
|
bits_to_file(output_filename, &output_bits);
|
||||||
|
|
||||||
|
|
||||||
|
// Affichage du signal dans un .wav
|
||||||
|
double* si = (double*)malloc(sizeof(double) * total_samples);
|
||||||
|
for (int i = 0; i < total_samples; i++) {
|
||||||
|
si[i] = cimag(s[i]);
|
||||||
|
}
|
||||||
|
write_wav("output.wav", si, total_samples);
|
||||||
|
|
||||||
|
// Plot
|
||||||
|
printf("Ploting...");
|
||||||
|
plot_t plot;
|
||||||
|
plot_init(&plot, 1400, 1400, 0.04);
|
||||||
|
plot_draw_constellation(&plot, qam.constellation, qam.M);
|
||||||
|
SDL_Color red = {255, 0, 0, 255};
|
||||||
|
plot_draw_points_animated(&plot, demodulate_points(&qam, s, nb_symbols), nb_symbols, red, 5);
|
||||||
|
|
||||||
|
// Libération mémoire
|
||||||
|
plot_close(&plot);
|
||||||
|
free_bit_array(&input_bits);
|
||||||
|
free_bit_array(&output_bits);
|
||||||
|
free(symbols);
|
||||||
|
free(s);
|
||||||
|
free_constellation(&qam);
|
||||||
|
free(output_filename);
|
||||||
|
|
||||||
|
return 0;
|
||||||
|
}
|
||||||
324
QAM/qamold/qambisbis.c
Normal file
324
QAM/qamold/qambisbis.c
Normal file
@ -0,0 +1,324 @@
|
|||||||
|
#include <math.h>
|
||||||
|
#include <stdint.h>
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <complex.h>
|
||||||
|
#include "../wav/wav.h"
|
||||||
|
#include "../files/files.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));
|
||||||
|
//int i = k * qam->N + n;
|
||||||
|
//r += s[i] * cexp(-2 * I * M_PI * qam->Fc * ((double)i / qam->Fs));
|
||||||
|
}
|
||||||
|
r /= qam->N;
|
||||||
|
|
||||||
|
// Distance euclidien de Ir et Qr pour avoir le point le plus proche de la constellation
|
||||||
|
|
||||||
|
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;
|
||||||
|
}
|
||||||
|
|
||||||
|
/*
|
||||||
|
double norm_factor = sqrt((double)(qam->M - 1) / 3.0);
|
||||||
|
double Ir = creal(r) * norm_factor / A;
|
||||||
|
double Qr = cimag(r) * norm_factor / A;
|
||||||
|
|
||||||
|
int i = (int)round((Ir + (sm - 1)) / 2.0);
|
||||||
|
int j = (int)round((Qr + (sm - 1)) / 2.0);
|
||||||
|
|
||||||
|
i = (i < 0) ? 0 : ((i >= sm) ? sm - 1 : i);
|
||||||
|
j = (j < 0) ? 0 : ((j >= sm) ? sm - 1 : j);
|
||||||
|
|
||||||
|
int id = i * sm + j;
|
||||||
|
|
||||||
|
int id = i_cl * sm + j_cl;
|
||||||
|
for (int b = 0; b < qam->k; b++) {
|
||||||
|
bits_hat[k * qam->k + (qam->k - 1 - b)] = (id >> b) & 1;
|
||||||
|
}
|
||||||
|
*/
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
double complex* demodulate_points(qam_system* qam, double complex* s, int nb_symbols) {
|
||||||
|
double complex* points = malloc(sizeof(double complex) * nb_symbols);
|
||||||
|
for (int k = 0; k < nb_symbols; k++) {
|
||||||
|
double complex r = 0;
|
||||||
|
for (int n = 0; n < qam->N; n++) {
|
||||||
|
int i = k * qam->N + n;
|
||||||
|
r += s[i] * cexp(-2 * I * M_PI * qam->Fc * ((double)i / qam->Fs));
|
||||||
|
}
|
||||||
|
r /= qam->N;
|
||||||
|
double norm_factor = sqrt((double)(qam->M - 1) / 3.0);
|
||||||
|
double Ir = creal(r);
|
||||||
|
double Qr = cimag(r);
|
||||||
|
|
||||||
|
|
||||||
|
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;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
double complex p = qam->constellation[i_cl][j_cl];
|
||||||
|
points[k] = creal(p) + I * cimag(p);
|
||||||
|
|
||||||
|
//points[k] = (int)Ir + I * (int)Qr;
|
||||||
|
}
|
||||||
|
return points;
|
||||||
|
}
|
||||||
|
|
||||||
|
// 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);
|
||||||
|
}
|
||||||
|
|
||||||
|
#include <stddef.h> // pour size_t
|
||||||
|
|
||||||
|
// Compare deux tableaux de bits (0/1) et retourne le pourcentage de fiabilité.
|
||||||
|
double compare_bits_and_get_reliability(const uint8_t *in_bits, size_t nb_bits_in, const uint8_t *out_bits, size_t nb_bits_out, size_t *erreurs) {
|
||||||
|
if (!in_bits || !out_bits) {
|
||||||
|
if (erreurs) *erreurs = (nb_bits_in < nb_bits_out) ? nb_bits_out : nb_bits_in;
|
||||||
|
return 0.0;
|
||||||
|
}
|
||||||
|
|
||||||
|
size_t n_min = (nb_bits_in < nb_bits_out) ? nb_bits_in : nb_bits_out;
|
||||||
|
size_t n_max = (nb_bits_in > nb_bits_out) ? nb_bits_in : nb_bits_out;
|
||||||
|
|
||||||
|
size_t err = 0;
|
||||||
|
for (size_t i = 0; i < n_min; ++i) {
|
||||||
|
if ((in_bits[i] & 1) != (out_bits[i] & 1)) err++;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (n_max != n_min) {
|
||||||
|
err += (n_max - n_min);
|
||||||
|
}
|
||||||
|
|
||||||
|
if (erreurs) *erreurs = err;
|
||||||
|
|
||||||
|
double total_compared = (double)n_max;
|
||||||
|
if (total_compared == 0.0) return 0.0;
|
||||||
|
|
||||||
|
double ber = (double)err / total_compared;
|
||||||
|
double reliability_percent = (1.0 - ber) * 100.0;
|
||||||
|
return reliability_percent;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
void affiche_constellation(qam_system* qam) {
|
||||||
|
int sm = (int)sqrt(qam->M);
|
||||||
|
for (int i = 0; i < sm; i++) {
|
||||||
|
for (int j = 0; j < sm; j++) {
|
||||||
|
double complex p = qam->constellation[i][j];
|
||||||
|
printf("(%d,%d) ", (int)creal(p), (int)cimag(p));
|
||||||
|
}
|
||||||
|
printf("\n");
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void affiche_points(double complex* r, int len, int len_samples) {
|
||||||
|
for (int i = 0; i < len; i += len_samples) {
|
||||||
|
for (int j = 0; j < len_samples; j++) {
|
||||||
|
double complex p = r[i + j];
|
||||||
|
printf("(%f,%f) ", (float)creal(p), (float)cimag(p));
|
||||||
|
}
|
||||||
|
printf("\n");
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
int main (int argc, char *argv[]) {
|
||||||
|
if (argc < 2) {
|
||||||
|
fprintf(stderr, "Utilisation: %s <fichier_entree>\n", argv[0]);
|
||||||
|
return 1;
|
||||||
|
}
|
||||||
|
|
||||||
|
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);
|
||||||
|
|
||||||
|
printf("Lecture du fichier...\n");
|
||||||
|
// Lecture du fichier et conversion en bits
|
||||||
|
const char *input_filename = argv[1];
|
||||||
|
bit_array input_bits = file_to_bits(input_filename);
|
||||||
|
size_t nb_symbols = input_bits.nb_bits / qam.k;
|
||||||
|
|
||||||
|
printf("Mise en forme des symboles...\n");
|
||||||
|
// Mise en forme des symboles
|
||||||
|
double complex *symbols = malloc(sizeof(double complex) * nb_symbols);
|
||||||
|
bits_to_symbols(&qam, input_bits.bits, input_bits.nb_bits, symbols);
|
||||||
|
|
||||||
|
printf("Modulation...\n");
|
||||||
|
// Modulation QAM
|
||||||
|
int total_samples = qam.N * nb_symbols;
|
||||||
|
double complex* s = (double complex*)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 avant échelle
|
||||||
|
double snr_dB = 10; // Signal to noise ratio
|
||||||
|
double snr_lin = pow(10.0, snr_dB / 10.0);
|
||||||
|
double sigma = sqrt(signal_power / snr_lin);
|
||||||
|
printf("Ajout du bruit... \n puissance du signal : %f\n SNR db : %f\n sigma : %f\n", signal_power, snr_dB, sigma);
|
||||||
|
add_noise(s, total_samples, 5);
|
||||||
|
|
||||||
|
printf("Demodulation...\n");
|
||||||
|
|
||||||
|
// Demodulation QAM
|
||||||
|
bit_array output_bits;
|
||||||
|
output_bits.nb_bits = input_bits.nb_bits;
|
||||||
|
output_bits.bits = (uint8_t*)malloc(output_bits.nb_bits);
|
||||||
|
demodulate(&qam, s, nb_symbols, output_bits.bits);
|
||||||
|
|
||||||
|
printf("Ecriture...\n");
|
||||||
|
// Ecriture du fichier de Demodulation
|
||||||
|
char *output_filename = make_output_filename(input_filename);
|
||||||
|
bits_to_file(output_filename, &output_bits);
|
||||||
|
|
||||||
|
//printf("Constelattion :\n");
|
||||||
|
//affiche_constellation(&qam);
|
||||||
|
//printf("Points de demodulation :\n");
|
||||||
|
//affiche_points(demodulate_points(&qam, s, nb_symbols), nb_symbols, qam.N);
|
||||||
|
|
||||||
|
// Affichage du signal dans un .wav
|
||||||
|
/*
|
||||||
|
double* si = (double*)malloc(sizeof(double) * total_samples);
|
||||||
|
for (int i = 0; i < total_samples; i++) {
|
||||||
|
si[i] = cimag(s[i]);
|
||||||
|
}
|
||||||
|
write_wav("output.wav", si, total_samples);
|
||||||
|
*/
|
||||||
|
|
||||||
|
|
||||||
|
size_t erreurs = 0;
|
||||||
|
double fiabilite = compare_bits_and_get_reliability( input_bits.bits, input_bits.nb_bits, output_bits.bits, output_bits.nb_bits, &erreurs);
|
||||||
|
|
||||||
|
printf("Résultat de la comparaison :\n");
|
||||||
|
printf(" Bits d'entrée : %zu\n", input_bits.nb_bits);
|
||||||
|
printf(" Bits de sortie: %zu\n", output_bits.nb_bits);
|
||||||
|
printf(" Erreurs : %zu\n", erreurs);
|
||||||
|
printf(" Fiabilité : %.4f %%\n", fiabilite);
|
||||||
|
|
||||||
|
// Libération mémoire
|
||||||
|
free_bit_array(&input_bits);
|
||||||
|
free_bit_array(&output_bits);
|
||||||
|
free(symbols);
|
||||||
|
free(s);
|
||||||
|
free_constellation(&qam);
|
||||||
|
free(output_filename);
|
||||||
|
|
||||||
|
return 0;
|
||||||
|
}
|
||||||
@ -3,7 +3,8 @@
|
|||||||
#include <stdlib.h>
|
#include <stdlib.h>
|
||||||
#include <string.h>
|
#include <string.h>
|
||||||
|
|
||||||
// Lire un fichier binaire et en faire un tableau de bits (0/1)
|
// files.c corrections : MSB-first within each byte
|
||||||
|
|
||||||
bit_array file_to_bits(const char *filename) {
|
bit_array file_to_bits(const char *filename) {
|
||||||
bit_array arr = {0};
|
bit_array arr = {0};
|
||||||
|
|
||||||
@ -13,7 +14,6 @@ bit_array file_to_bits(const char *filename) {
|
|||||||
return arr;
|
return arr;
|
||||||
}
|
}
|
||||||
|
|
||||||
// Taille fichier
|
|
||||||
fseek(f, 0, SEEK_END);
|
fseek(f, 0, SEEK_END);
|
||||||
long file_size = ftell(f);
|
long file_size = ftell(f);
|
||||||
rewind(f);
|
rewind(f);
|
||||||
@ -23,27 +23,19 @@ bit_array file_to_bits(const char *filename) {
|
|||||||
return arr;
|
return arr;
|
||||||
}
|
}
|
||||||
|
|
||||||
// Lire tous les octets
|
|
||||||
uint8_t *raw = (uint8_t*)malloc(file_size);
|
uint8_t *raw = (uint8_t*)malloc(file_size);
|
||||||
if (!raw) {
|
if (!raw) { fclose(f); return arr; }
|
||||||
fclose(f);
|
|
||||||
return arr;
|
|
||||||
}
|
|
||||||
fread(raw, 1, file_size, f);
|
fread(raw, 1, file_size, f);
|
||||||
fclose(f);
|
fclose(f);
|
||||||
|
|
||||||
// Convertir en bits (0/1 dans uint8_t)
|
|
||||||
arr.nb_bits = (size_t)file_size * 8;
|
arr.nb_bits = (size_t)file_size * 8;
|
||||||
arr.bits = (uint8_t*)malloc(arr.nb_bits);
|
arr.bits = (uint8_t*)malloc(arr.nb_bits * sizeof(uint8_t));
|
||||||
if (!arr.bits) {
|
if (!arr.bits) { free(raw); arr.nb_bits = 0; return arr; }
|
||||||
free(raw);
|
|
||||||
arr.nb_bits = 0;
|
|
||||||
return arr;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
// MSB-first in each byte: bit 7 -> index 0, bit 0 -> index 7
|
||||||
for (size_t i = 0; i < (size_t)file_size; i++) {
|
for (size_t i = 0; i < (size_t)file_size; i++) {
|
||||||
for (int b = 0; b < 8; b++) {
|
for (int b = 0; b < 8; b++) {
|
||||||
arr.bits[i * 8 + b] = (raw[i] >> b) & 1u;
|
arr.bits[i * 8 + b] = (raw[i] >> (7 - b)) & 1u;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -51,7 +43,6 @@ bit_array file_to_bits(const char *filename) {
|
|||||||
return arr;
|
return arr;
|
||||||
}
|
}
|
||||||
|
|
||||||
// Transformer un tableau de bits (0/1) en fichier binaire
|
|
||||||
int bits_to_file(const char *filename, const bit_array *arr) {
|
int bits_to_file(const char *filename, const bit_array *arr) {
|
||||||
if (!arr || !arr->bits) return -1;
|
if (!arr || !arr->bits) return -1;
|
||||||
|
|
||||||
@ -59,19 +50,17 @@ int bits_to_file(const char *filename, const bit_array *arr) {
|
|||||||
uint8_t *raw = (uint8_t*)calloc(nb_bytes, 1);
|
uint8_t *raw = (uint8_t*)calloc(nb_bytes, 1);
|
||||||
if (!raw) return -1;
|
if (!raw) return -1;
|
||||||
|
|
||||||
|
// MSB-first packing: index bit 0 -> place at position 7 of first byte
|
||||||
for (size_t i = 0; i < arr->nb_bits; i++) {
|
for (size_t i = 0; i < arr->nb_bits; i++) {
|
||||||
if (arr->bits[i]) {
|
if (arr->bits[i]) {
|
||||||
raw[i / 8] |= (1u << (i % 8));
|
size_t byte_idx = i / 8;
|
||||||
|
int bit_in_byte = i % 8; // 0..7
|
||||||
|
raw[byte_idx] |= (1u << (7 - bit_in_byte));
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
FILE *f = fopen(filename, "wb");
|
FILE *f = fopen(filename, "wb");
|
||||||
if (!f) {
|
if (!f) { perror("fopen"); free(raw); return -1; }
|
||||||
perror("fopen");
|
|
||||||
free(raw);
|
|
||||||
return -1;
|
|
||||||
}
|
|
||||||
|
|
||||||
fwrite(raw, 1, nb_bytes, f);
|
fwrite(raw, 1, nb_bytes, f);
|
||||||
fclose(f);
|
fclose(f);
|
||||||
free(raw);
|
free(raw);
|
||||||
|
|||||||
@ -19,3 +19,4 @@ char* make_output_filename(const char *input_filename);
|
|||||||
|
|
||||||
// Libérer la mémoire du bit_array
|
// Libérer la mémoire du bit_array
|
||||||
void free_bit_array(bit_array *arr);
|
void free_bit_array(bit_array *arr);
|
||||||
|
|
||||||
|
|||||||
Reference in New Issue
Block a user