diff --git a/FFT/ffttp.c b/FFT/ffttp.c new file mode 100644 index 0000000..30552f9 --- /dev/null +++ b/FFT/ffttp.c @@ -0,0 +1,79 @@ +#include +#include +#include +#include +#include + +complex float* fft(complex float* p, uint32_t N) { + complex float w = cexpf(-2 * I * M_PI / N); + + if (N == 1) + return p; + + complex float* pe = (complex float*)malloc(sizeof(complex float) * (N / 2)); + complex float* po = (complex float*)malloc(sizeof(complex float) * (N / 2)); + for (int i = 0; i < N - 1; i += 2) { + pe[i / 2] = p[i]; + po[i / 2] = p[i + 1]; + } + + complex float* xe = fft(pe, N / 2); + complex float* xo = fft(po, N / 2); + + complex float* ps = (complex float*)malloc(sizeof(complex float) * N); + + for (int i = 0; i < N / 2; i++) { + ps[i] = xe[i] + cpowf(w, (complex float)i) * xo[i]; + ps[i + N / 2] = xe[i] - cpowf(w, (complex float)i) * xo[i]; + } + + return ps; +} + +int main() { + FILE *f = fopen("signal.csv","r"); + if (!f) { perror("signal.csv"); return 1; } + + double t, a; + int capacity = 1024, n = 0; + complex float *samples = malloc(sizeof(complex float) * capacity); + + while (fscanf(f, "%lf;%lf", &t, &a) == 2) { + if (n >= capacity) { + capacity *= 2; + samples = realloc(samples, sizeof(complex float)*capacity); + } + samples[n++] = a; + } + fclose(f); + + int fsize = 1 << (int)floorf(log2f((float)n)); + + complex float* cwav_data = (complex float*)malloc(sizeof(complex float) * fsize); + for (int i = 0; i < fsize; i++) + cwav_data[i] = samples[i]; + + complex float* ft = fft(cwav_data, fsize); + + double Fe = 400.0; + + FILE *gnuplot = popen("gnuplot -persistent", "w"); + fprintf(gnuplot, "set title 'fft'\n"); + fprintf(gnuplot, "set logscale x\n"); + fprintf(gnuplot, "set logscale y\n"); + fprintf(gnuplot, "plot '-' with linespoints title 'fft'\n"); + + for (int i = 0; i < fsize/2; i++) { + double freq = (double)i * Fe / fsize; + fprintf(gnuplot, "%f %f\n", freq ,cabsf(ft[i])); + } + fprintf(gnuplot, "e\n"); + pclose(gnuplot); + + free(samples); + free(cwav_data); + free(ft); + + return 0; +} + diff --git a/FFT/fftwav.c b/FFT/fftwav.c index 9fe8db6..2cbd6d0 100644 --- a/FFT/fftwav.c +++ b/FFT/fftwav.c @@ -3,7 +3,7 @@ #include #include #include -#include "../WAV/wav.h" +#include "../wav/wav.h" complex float* fft(complex float* p, uint32_t N) { complex float w = cexpf(-2 * I * M_PI / N); diff --git a/FFT/signal.csv b/FFT/signal.csv new file mode 100644 index 0000000..265db52 --- /dev/null +++ b/FFT/signal.csv @@ -0,0 +1,476 @@ +0;0.164159775 +0.000315;0.164159775 +0.00063;0.164159775 +0.000945;0.164159775 +0.00126;0.164159775 +0.001575;0.164159775 +0.00189;1.60616973 +0.002205;1.60616973 +0.00252;1.60616973 +0.002835;1.60616973 +0.00315;1.601180077 +0.003465;1.60616973 +0.00378;1.60616973 +0.004095;2.47935915 +0.00441;2.439441919 +0.004725;2.439441919 +0.00504;2.439441919 +0.005355;2.434452266 +0.00567;2.434452266 +0.005985;2.434452266 +0.0063;2.439441919 +0.006615;2.34962815 +0.00693;2.34962815 +0.007245;2.34962815 +0.00756;2.34962815 +0.007875;2.34962815 +0.00819;2.34962815 +0.008505;2.34962815 +0.00882;2.34962815 +0.009135;1.371656001 +0.00945;1.371656001 +0.009765;1.376645654 +0.01008;1.371656001 +0.010395;1.371656001 +0.01071;1.376645654 +0.011025;1.371656001 +0.01134;1.371656001 +0.011655;-0.120250493 +0.01197;-0.115260839 +0.012285;-0.110271186 +0.0126;-0.110271186 +0.012915;-0.115260839 +0.01323;-0.110271186 +0.013545;-0.115260839 +0.01386;-0.110271186 +0.014175;-1.552281141 +0.01449;-1.547291487 +0.014805;-1.552281141 +0.01512;-1.547291487 +0.015435;-1.547291487 +0.01575;-1.547291487 +0.016065;-1.547291487 +0.01638;-1.547291487 +0.016695;-2.395532638 +0.01701;-2.390542984 +0.017325;-2.390542984 +0.01764;-2.390542984 +0.017955;-2.390542984 +0.01827;-2.390542984 +0.018585;-2.390542984 +0.0189;-2.390542984 +0.019215;-2.305718869 +0.01953;-2.300729215 +0.019845;-2.300729215 +0.02016;-2.300729215 +0.020475;-2.305718869 +0.02079;-2.305718869 +0.021105;-2.300729215 +0.02142;-2.300729215 +0.021735;-1.322757065 +0.02205;-1.322757065 +0.022365;-1.322757065 +0.02268;-1.322757065 +0.022995;-1.322757065 +0.02331;-1.322757065 +0.023625;-1.322757065 +0.02394;-1.322757065 +0.024255;0.169149429 +0.02457;0.164159775 +0.024885;0.164159775 +0.0252;0.169149429 +0.025515;0.164159775 +0.02583;0.164159775 +0.026145;0.164159775 +0.02646;0.164159775 +0.026775;1.60616973 +0.02709;1.60616973 +0.027405;1.60616973 +0.02772;1.60616973 +0.028035;1.60616973 +0.02835;1.601180077 +0.028665;1.601180077 +0.02898;1.601180077 +0.029295;2.439441919 +0.02961;2.439441919 +0.029925;2.439441919 +0.03024;2.439441919 +0.030555;2.439441919 +0.03087;2.439441919 +0.031185;2.439441919 +0.0315;2.444431573 +0.031815;2.34962815 +0.03213;2.34962815 +0.032445;2.34962815 +0.03276;2.34962815 +0.033075;2.34962815 +0.03339;2.34962815 +0.033705;2.34962815 +0.03402;2.34962815 +0.034335;1.371656001 +0.03465;1.371656001 +0.034965;1.371656001 +0.03528;1.371656001 +0.035595;1.371656001 +0.03591;1.371656001 +0.036225;1.371656001 +0.03654;1.371656001 +0.036855;-0.115260839 +0.03717;-0.115260839 +0.037485;-0.115260839 +0.0378;-0.110271186 +0.038115;-0.110271186 +0.03843;-0.110271186 +0.038745;-0.110271186 +0.03906;-0.110271186 +0.039375;-1.557270795 +0.03969;-1.557270795 +0.040005;-1.557270795 +0.04032;-1.552281141 +0.040635;-1.557270795 +0.04095;-1.552281141 +0.041265;-1.552281141 +0.04158;-1.552281141 +0.041895;-2.390542984 +0.04221;-2.390542984 +0.042525;-2.38555333 +0.04284;-2.390542984 +0.043155;-2.390542984 +0.04347;-2.390542984 +0.043785;-2.390542984 +0.0441;-2.295739561 +0.044415;-2.300729215 +0.04473;-2.300729215 +0.045045;-2.300729215 +0.04536;-2.300729215 +0.045675;-2.300729215 +0.04599;-2.305718869 +0.046305;-2.300729215 +0.04662;-1.317767411 +0.046935;-1.322757065 +0.04725;-1.322757065 +0.047565;-1.322757065 +0.04788;-1.322757065 +0.048195;-1.317767411 +0.04851;-1.322757065 +0.048825;-1.322757065 +0.04914;0.174139082 +0.049455;0.169149429 +0.04977;0.164159775 +0.050085;0.164159775 +0.0504;0.169149429 +0.050715;0.164159775 +0.05103;0.164159775 +0.051345;0.164159775 +0.05166;1.60616973 +0.051975;1.601180077 +0.05229;1.601180077 +0.052605;1.601180077 +0.05292;1.601180077 +0.053235;1.601180077 +0.05355;1.596190423 +0.053865;1.601180077 +0.05418;2.439441919 +0.054495;2.439441919 +0.05481;2.439441919 +0.055125;2.434452266 +0.05544;2.434452266 +0.055755;2.439441919 +0.05607;2.434452266 +0.056385;2.434452266 +0.0567;2.34962815 +0.057015;2.34962815 +0.05733;2.34962815 +0.057645;2.34962815 +0.05796;2.34962815 +0.058275;2.34962815 +0.05859;2.34962815 +0.058905;2.34962815 +0.05922;1.366666347 +0.059535;1.371656001 +0.05985;1.371656001 +0.060165;1.371656001 +0.06048;1.371656001 +0.060795;1.376645654 +0.06111;1.371656001 +0.061425;1.371656001 +0.06174;-0.115260839 +0.062055;-0.115260839 +0.06237;-0.115260839 +0.062685;-0.115260839 +0.063;-0.110271186 +0.063315;-0.115260839 +0.06363;-0.110271186 +0.063945;-0.110271186 +0.06426;-1.552281141 +0.064575;-1.552281141 +0.06489;-1.552281141 +0.065205;-1.547291487 +0.06552;-1.547291487 +0.065835;-1.547291487 +0.06615;-1.547291487 +0.066465;-1.547291487 +0.06678;-2.390542984 +0.067095;-2.390542984 +0.06741;-2.390542984 +0.067725;-2.380563676 +0.06804;-2.390542984 +0.068355;-2.38555333 +0.06867;-2.38555333 +0.068985;-2.390542984 +0.0693;-2.305718869 +0.069615;-2.300729215 +0.06993;-2.305718869 +0.070245;-2.300729215 +0.07056;-2.305718869 +0.070875;-2.305718869 +0.07119;-2.305718869 +0.071505;-2.300729215 +0.07182;-1.322757065 +0.072135;-1.322757065 +0.07245;-1.317767411 +0.072765;-1.322757065 +0.07308;-1.322757065 +0.073395;-1.322757065 +0.07371;-1.332736373 +0.074025;-1.322757065 +0.07434;0.169149429 +0.074655;0.169149429 +0.07497;0.169149429 +0.075285;0.164159775 +0.0756;0.169149429 +0.075915;0.164159775 +0.07623;0.164159775 +0.076545;0.164159775 +0.07686;1.60616973 +0.077175;1.601180077 +0.07749;1.601180077 +0.077805;1.601180077 +0.07812;1.601180077 +0.078435;1.601180077 +0.07875;1.601180077 +0.079065;1.601180077 +0.07938;2.439441919 +0.079695;2.439441919 +0.08001;2.439441919 +0.080325;2.434452266 +0.08064;2.439441919 +0.080955;2.434452266 +0.08127;2.434452266 +0.081585;2.314700574 +0.0819;2.34962815 +0.082215;2.34962815 +0.08253;2.34962815 +0.082845;2.34962815 +0.08316;2.34962815 +0.083475;2.34962815 +0.08379;2.34962815 +0.084105;1.361676693 +0.08442;1.371656001 +0.084735;1.371656001 +0.08505;1.371656001 +0.085365;1.371656001 +0.08568;1.376645654 +0.085995;1.371656001 +0.08631;1.371656001 +0.086625;-0.120250493 +0.08694;-0.115260839 +0.087255;-0.115260839 +0.08757;-0.115260839 +0.087885;-0.110271186 +0.0882;-0.110271186 +0.088515;-0.110271186 +0.08883;-0.115260839 +0.089145;-1.557270795 +0.08946;-1.552281141 +0.089775;-1.547291487 +0.09009;-1.547291487 +0.090405;-1.547291487 +0.09072;-1.547291487 +0.091035;-1.547291487 +0.09135;-1.552281141 +0.091665;-2.395532638 +0.09198;-2.390542984 +0.092295;-2.390542984 +0.09261;-2.390542984 +0.092925;-2.390542984 +0.09324;-2.390542984 +0.093555;-2.38555333 +0.09387;-2.390542984 +0.094185;-2.300729215 +0.0945;-2.300729215 +0.094815;-2.300729215 +0.09513;-2.300729215 +0.095445;-2.300729215 +0.09576;-2.300729215 +0.096075;-2.300729215 +0.09639;-2.300729215 +0.096705;-1.317767411 +0.09702;-1.317767411 +0.097335;-1.322757065 +0.09765;-1.322757065 +0.097965;-1.322757065 +0.09828;-1.322757065 +0.098595;-1.327746719 +0.09891;-1.322757065 +0.099225;0.169149429 +0.09954;0.169149429 +0.099855;0.164159775 +0.10017;0.164159775 +0.100485;0.164159775 +0.1008;0.164159775 +0.101115;0.159170121 +0.10143;0.169149429 +0.101745;1.601180077 +0.10206;1.596190423 +0.102375;1.601180077 +0.10269;1.601180077 +0.103005;1.601180077 +0.10332;1.601180077 +0.103635;1.596190423 +0.10395;1.596190423 +0.104265;2.434452266 +0.10458;2.439441919 +0.104895;2.439441919 +0.10521;2.434452266 +0.105525;2.434452266 +0.10584;2.434452266 +0.106155;2.434452266 +0.10647;2.429462612 +0.106785;2.354617804 +0.1071;2.34962815 +0.107415;2.34962815 +0.10773;2.34962815 +0.108045;2.34962815 +0.10836;2.354617804 +0.108675;2.34962815 +0.10899;2.34962815 +0.109305;1.371656001 +0.10962;1.376645654 +0.109935;1.371656001 +0.11025;1.371656001 +0.110565;1.371656001 +0.11088;1.371656001 +0.111195;1.376645654 +0.11151;1.371656001 +0.111825;-0.115260839 +0.11214;-0.115260839 +0.112455;-0.110271186 +0.11277;-0.110271186 +0.113085;-0.110271186 +0.1134;-0.115260839 +0.113715;-0.110271186 +0.11403;-0.115260839 +0.114345;-1.552281141 +0.11466;-1.557270795 +0.114975;-1.552281141 +0.11529;-1.552281141 +0.115605;-1.552281141 +0.11592;-1.552281141 +0.116235;-1.552281141 +0.11655;-1.552281141 +0.116865;-2.390542984 +0.11718;-2.390542984 +0.117495;-2.390542984 +0.11781;-2.390542984 +0.118125;-2.390542984 +0.11844;-2.390542984 +0.118755;-2.390542984 +0.11907;-2.390542984 +0.119385;-2.305718869 +0.1197;-2.300729215 +0.120015;-2.300729215 +0.12033;-2.305718869 +0.120645;-2.300729215 +0.12096;-2.300729215 +0.121275;-2.305718869 +0.12159;-1.277850181 +0.121905;-1.322757065 +0.12222;-1.322757065 +0.122535;-1.322757065 +0.12285;-1.322757065 +0.123165;-1.322757065 +0.12348;-1.322757065 +0.123795;-1.327746719 +0.12411;0.18411839 +0.124425;0.164159775 +0.12474;0.164159775 +0.125055;0.164159775 +0.12537;0.164159775 +0.125685;0.164159775 +0.126;0.164159775 +0.126315;0.164159775 +0.12663;1.616149038 +0.126945;1.60616973 +0.12726;1.60616973 +0.127575;1.60616973 +0.12789;1.601180077 +0.128205;1.601180077 +0.12852;1.601180077 +0.128835;1.601180077 +0.12915;2.439441919 +0.129465;2.439441919 +0.12978;2.434452266 +0.130095;2.439441919 +0.13041;2.439441919 +0.130725;2.439441919 +0.13104;2.434452266 +0.131355;2.434452266 +0.13167;2.34962815 +0.131985;2.34962815 +0.1323;2.34962815 +0.132615;2.34962815 +0.13293;2.34962815 +0.133245;2.34962815 +0.13356;2.34962815 +0.133875;2.34962815 +0.13419;1.371656001 +0.134505;1.371656001 +0.13482;1.371656001 +0.135135;1.371656001 +0.13545;1.371656001 +0.135765;1.371656001 +0.13608;1.376645654 +0.136395;1.371656001 +0.13671;-0.120250493 +0.137025;-0.115260839 +0.13734;-0.115260839 +0.137655;-0.115260839 +0.13797;-0.115260839 +0.138285;-0.110271186 +0.1386;-0.110271186 +0.138915;-0.115260839 +0.13923;-1.552281141 +0.139545;-1.547291487 +0.13986;-1.552281141 +0.140175;-1.547291487 +0.14049;-1.547291487 +0.140805;-1.552281141 +0.14112;-1.552281141 +0.141435;-1.552281141 +0.14175;-2.390542984 +0.142065;-2.390542984 +0.14238;-2.390542984 +0.142695;-2.390542984 +0.14301;-2.390542984 +0.143325;-2.390542984 +0.14364;-2.390542984 +0.143955;-2.390542984 +0.14427;-2.300729215 +0.144585;-2.300729215 +0.1449;-2.305718869 +0.145215;-2.300729215 +0.14553;-2.300729215 +0.145845;-2.300729215 +0.14616;-2.300729215 +0.146475;-2.300729215 +0.14679;-1.317767411 +0.147105;-1.317767411 +0.14742;-1.322757065 +0.147735;-1.322757065 +0.14805;-1.322757065 +0.148365;-1.317767411 +0.14868;-1.322757065 +0.148995;-1.322757065 +0.14931;0.169149429 +0.149625;0.169149429 diff --git a/QAM/qam b/QAM/qam index 470d3b1..836a1e4 100755 Binary files a/QAM/qam and b/QAM/qam differ diff --git a/QAM/qam.c b/QAM/qam.c new file mode 100644 index 0000000..17351f1 --- /dev/null +++ b/QAM/qam.c @@ -0,0 +1,273 @@ +#include +#include +#include +#include +#include +#include + +#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; +} diff --git a/QAM/qamold/qamclean.c b/QAM/qamold/qamclean.c new file mode 100644 index 0000000..a63deb4 --- /dev/null +++ b/QAM/qamold/qamclean.c @@ -0,0 +1,184 @@ +#include +#include +#include +#include +#include +#include + +#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; +} + +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; + } + + // 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 = 10; // 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); + + // 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); + + // Libération mémoire + free(input_bits); + free(output_bits); + free(symbols); + free(s); + free_constellation(&qam); + + return 0; +} diff --git a/QAM/qamold/qamtest.c b/QAM/qamold/qamtest.c new file mode 100644 index 0000000..09bd9d4 --- /dev/null +++ b/QAM/qamold/qamtest.c @@ -0,0 +1,328 @@ +#include +#include +#include +#include +#include +#include "../files/files.h" +#include + +#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; + } + } +} + +// Demodulation QAM avec porteuse estimé +void demodulate2(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]; + } + 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; + } + } +} +/* +// Pour j = 1 dans la def de wiki sur l'autocorrelation +double fc_autocorrelation_1_rad(double complex* s, int N) { + double complex r = 0; + for (int n = 0; n < N - 1; n++) { + r += s[n] * conj(s[n + 1]); + } + r /= (N - 1); + double om = -carg(r); // rad / sample + return om; +} + +// Autocorrelation pour trouver la Fc (fréquence de la porteuse) +double fc_autocorrelation_multilag_rad(double complex* s, int N, int max_lag) { + double sum_phase = 0; + for (int k = 1; k <= max_lag; k++) { + double complex rk = 0; + for (int n = 0; n < N - k; n++) { + rk += s[n] * conj(s[n + k]); + } + rk /= (double)(N - k); + sum_phase += -carg(rk) / k; + } + return sum_phase / max_lag; // rad/sample +} + +// Correction du signal pour enlever l'offset de la porteuse +void signal_correction(double complex* s, int total_samples, double om_hat) { + for (int n = 0; n < total_samples; n++) { + s[n] *= cexp(-I * om_hat * n); + } +} +*/ + +void blind_carrier(double complex* s, int N, int M) { + int k = (int)sqrt(M); // puissance à élever + double complex sum = 0; + for(int n = 0; n < N; n++) + sum += cpow(s[n], k); + + double phi_hat = carg(sum) / k; + + for(int n = 0; n < N; n++) + s[n] *= cexp(-I * phi_hat); +} + +// 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); +} + +// Compare deux tableaux de bits (0/1) et retourne le pourcentage de fiabilité. +double taux_erreur_bits(const uint8_t *in_bits, size_t nb_bits_in, + const uint8_t *out_bits, size_t nb_bits_out) { + if (!in_bits || !out_bits) + return 1.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 = n_max - n_min; + + for (size_t i = 0; i < n_min; i++) + err += ((in_bits[i] ^ out_bits[i]) & 1); + + return n_max ? (double)err / n_max : 0.0; +} + + +int main (int argc, char *argv[]) { + if (argc < 2) { + fprintf(stderr, "Utilisation: %s \n", argv[0]); + return 1; + } + + qam_system qam; + qam.M = 16; + qam.k = (int)log2((double)(qam.M)); + qam.Fs = 44100; + qam.Ts = 0.3; + 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, 0); + + + // Demodulation QAM + //printf("Demodulation...\n"); + //bit_array output_bits; + //output_bits.nb_bits = input_bits.nb_bits; + //output_bits.bits = (uint8_t*)malloc(output_bits.nb_bits * sizeof(uint8_t)); + //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); + + //double erreurs = taux_erreur_bits(input_bits.bits, input_bits.nb_bits, output_bits.bits, output_bits.nb_bits); + + //printf("Comparaison :\n"); + //printf(" Erreurs : %f\n", erreurs * 100); + + /* + double om_hat = fc_autocorrelation_1_rad(s, total_samples); + printf("Estimation de fc 1: %f\n", (om_hat * qam.Fs) / (2 * M_PI)); + + om_hat = fc_autocorrelation_multilag_rad(s, total_samples, 10); + printf("Estimation de fc multilag: %f\n", (om_hat * qam.Fs) / (2 * M_PI)); + + // Correction du signal avec Fc (en rad/sample) trouvé + signal_correction(s, total_samples, om_hat); + + printf("Demodulation...\n"); + demodulate2(&qam, s, nb_symbols, output_bits.bits); + */ + + // Avant la demodulation + printf("Test de blind carrier correction...\n"); + + // Paramètres CMA simples + double mu = 0.001; // pas d'adaptation + int num_iter = 100; // nombre d'itérations + + // Appel de la fonction blind carrier correction + blind_carrier(s, total_samples, qam.M); + + // Ensuite, utilise ton démodulateur actuel + bit_array output_bits; + output_bits.nb_bits = input_bits.nb_bits; + output_bits.bits = (uint8_t*)malloc(output_bits.nb_bits * sizeof(uint8_t)); + + demodulate2(&qam, s, nb_symbols, output_bits.bits); + + // Comparaison avec bits originaux + double erreurs = taux_erreur_bits(input_bits.bits, input_bits.nb_bits, output_bits.bits, output_bits.nb_bits); + printf("Taux d'erreur après blind carrier correction: %f %%\n", erreurs * 100); + + // Ecriture du fichier de Demodulation + printf("Ecriture...\n"); + char *output_filename = make_output_filename(input_filename); + bits_to_file(output_filename, &output_bits); + + // 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; +}