From 876d37d1ab139673a74f81b558ef711145af7f31 Mon Sep 17 00:00:00 2001 From: zeefaad Date: Tue, 7 Oct 2025 12:58:52 +0200 Subject: [PATCH] qam edits --- FFT/ffttp.c | 79 +++++++ FFT/fftwav.c | 2 +- FFT/signal.csv | 476 ++++++++++++++++++++++++++++++++++++++++++ QAM/qam | Bin 21408 -> 20480 bytes QAM/qam.c | 273 ++++++++++++++++++++++++ QAM/qamold/qamclean.c | 184 ++++++++++++++++ QAM/qamold/qamtest.c | 328 +++++++++++++++++++++++++++++ 7 files changed, 1341 insertions(+), 1 deletion(-) create mode 100644 FFT/ffttp.c create mode 100644 FFT/signal.csv create mode 100644 QAM/qam.c create mode 100644 QAM/qamold/qamclean.c create mode 100644 QAM/qamold/qamtest.c 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 470d3b124bb1ac21d085f85183bae8f1213cae4f..836a1e423b1ce2fbb2663c4a514c2a2d2c2f0734 100755 GIT binary patch delta 5885 zcmZ`d4R}-6b@$1Ze}v7GjBSZu3mRv}*cSc+9LGimzwpQqHVz)7qX@;8u!?OtmQ36f zEGM2AuP;u48#cCPtD7%h#$g#oEh=p`IZ2D-B*id?w{z}$ z4?FbR{l0h4JwNxHbMHC#y!YNtx&PFqof zdpnQjaAix6r>Lk7iCM)tdXg*HGHOtr{T+^<_BCsrWQU-{6krbX2P7 zJftm#{j1lPtK@rONF|dHeaFzWI;z=-YBqG})O>Q%h7@;=kYRfA72BnoXWswwi&oFW zOTYWSjX!gA`XAd+<>+iL?eur|^_TW<+EBWosw@~Nt03qgd+6g?4(_h(~huWQ5~unDQ;V+~go8kj;L5j_|gh(MdxGH_B0mQEkh=snG< zoKDyAXcN$ph*_0Fk7iowk<6^@IFg|uHH*HS>EJB%^UQ6Wk=B{4w9}N8ZPyg0gG4&* z0_yWXUB7^uN&m)VspF!^fJ zX8K^xa)Ya$KAN+#_-&X$`YKz+p}JxmiqZ$6Y*DWJM?jQ38}!LGMNxiYe36diG;x>c zbdJ>!xkEFT)EFMVLmQXepL=B`Av65`F`vi}M4{`BCLMHO$*~nqV2H91_3=elef*7| z!*f~Xm+{n|`_1arUnt5@lmmKnhA)Q4oc*^eRAn=~MVL;KKNL7hX?{B+J%YS?O{uUM@suH6I^e9~W?9bEP~w~-!R?g;Ps9F{K3a@cy_hhI0gVZB2tF#2g*mY;&gOA({xvhhCjn20+$8@R*)$$X zMGHg199WsF72Z4YmbTNvY#3t!d+vGX&e@JRs>u&!>%6408yPa2COr%ZlEL2neiz< zklmMyCije~`B(H7hG;zLfHMsCi{d~=)|DB4uM4tcuYh^tn#Rn7Jh33nA|!{v1Ms3& zi1N>PAkIU=B;b>^!U+G205}+iDUe!fCaw8Q%RW&Ot^!?v%y$U!!})BC8so4qDJ+b_ z!H|TDL$FF2bitD4))&`~*R)RX&;IBR+KB=qS}Q!qKl=jUVCT5pIx3lAV3u1iN@hJp zp87Ze`_?YXsf9%_0K>E8QHYJ1?IxVQ8$X9g$HBh=>>%BLH?v&N&k69$=w5<{&LZIU z8NW(zuPhB$X{3ZWEJQIWcfi7++U(0g3y}kCi@USo&kU~Sh;-o9b%|j}0=n~pVANbA? zGs4M{FnUqT*1|t?7{)SlT2HrFGEH^42Uwy!j6gok79_ul%g}di+wXxPT$NAblKFGv zQluG0SwFAr^-tLP`*ET?j3BJe5am~y1uX_2_EU_qA5Mw!FbpBUYOv|F-#5fFTfZm8Hq(Ufdt0pqsrtBAWGv> zd~CL0ELp&f4_)Whj!WTpWG19mxDkk$8ZMj&$3oXz!-S6fv#eWTbt0BLq*^*VFDXIs zaGe$9z7ah%=fGF_H2SK>7EDC;G7S1KW~QdBKvi}=7I{#2TaejUm(dE-_?M|zlg_I( z$pl-wB?AJ4k27P4&4x5Xg|*E0BQM7yPF;Po9A?M@{gEDEick=TQA%Ae_oedp|CvI0 zFomF&CGn9LSVrPnSw`+-3bEJi#ggHY8!T3%uxq38X_SS%JFoBEGc3+UK*cD=A|=`g z6gLj%k{!bl$8Z-_53fT5LDGo=uMyJE}+EiReZ`=vo$Vz>UIpvgpxuHq&1M zjyC@#i+*?AALSiGHfWyVk8Oudfce#rWYKrm*_W3ihrgkqa4U?buE84xZEwn=`6bnr z7NlJP8tw$p>!#tG6II=QViXM;6Ntld$DcThcy@dd9Vih>{YVB^rQBHp zC^3E$k9`tkq zx_UJsk9?!+u|TJNx4+BZebCGKYMqo_6gA z5k5y)xT>nqoApHNTbSGbcjQdwZQ?>c>{ zaus)x{A!zvDg>70h z4nbz^CZXS|uQzPs=#~0at21D^*zk)0o2U&mGHA}+>@ZwSqYrFeWjLQkzqh&q6q6Ow!cSD!k0AetfrHtK5(x9-3*Z!yk|dtrpe!e?;mUvnzn~Dx_)KF z4kYleYxHkhnK=1GJH)>QI_NL{gAKo>O|ZgWFE)0T0Nz1&SMo3&mwIoQmIIxIpR1;0 zjaK7BfUjZc;DJ$&0ZJ%l+k$owuEk#z%KH zxVCSl`?eL^eZ5|%u@=jouu?GA=KZHL@#9f#aLkH7OSvr4;+(g@*- zGR?cY5^x*U76{IxOdA?UmAxTy(9;(T`aRw5?tnk&1xp9vC8F}&Zh{MUw*|U`A#Z1= zC*%)wgG75U;O_8rw|9d0_9l>R_jkMdg5Gw*Kjsew-JyV6U3QSJ7jvC$fvz4;uh)${ zSa^3kfmMd?(c13q3bglu33TDHU{Ld}%j55+4{l#Y|5IFM4nE!+TDTCJy?uqV~$>n^_yV0@qK#N05?dv7DuLDt~zK<%SMZS2K$|=(G0dxBP}+rpT6F* qIq!@{Z^qXfG<<`m)9gpqnUal+K%O?xEst!-{Q5Q(&8K6og8u_MmXU1$ literal 21408 zcmeHPe{@vUoqv-b63{TgiUk!JUC@X!@_R+33_l(`S|SLEf<=Z&GJ%o&>dXYtg_;^c zh9NGttHrL{vToahd#F-H&04DyP!cJ%0nrAvYoOT92uj3Xsj|*~KKIAu&5-Ho*+2H| zIm^v??|Z*L@ArQ1efQVPynD=fvvV8{L#JHh8wRnCL6YK84A(rB0q__%7z1&?+!$t@ z2Ru*EB;O+exH1*xro#e7Cz*P?6+aVRrsOo#9uiV7TdKc6ax^rGBu>42W)+?Z3x`QL z4K+Vi-WbU@GGL=HY8V>Y^%$F7Q>EUnl5{vm)!U<>*2{qy^pD?jpb`f@qIX9OXIgLh) zG|hJy4uSUvjVVXrnorq_aIrb4|JbUpJ-zO}%G1YI%=^dNKkJ^jcJ3Q3*H;J1Cr!A% zx^i@NptfP{=(RUY8a-)(J5=W$C%T0Ty0j+{7+knAQ#lBImHUG%^qpDg*MaWBm6^Ig zxH6Ln{KvTZD)(Di=v7(h+p^H_%|ib)3;j$M`W;#5xmnulAF{|lorT_+g+4nA{iQ7Q z1zG5ev(OKMKFAnqbWssWk28jMZfS|P(jWA%2!z7^U}?$B>bhEgsjs}+Z+N{cYU*me zp|CF)_IjmgpMn+swe?1YuRLVb*R3%s>Ow~7JHfC~UAJN!Z${9+CQw^xglg-9f!c7D zQQr`Ts6QC23!)I%uUO@+Sh>nuHy;$j!C=M8pixy<@2@3R4G{FfRE@8? zx~>As7#k{-hc&(cg-{u2sFJ+iinTs(RiM^a9k?52!~W`Ocq){uDyr*3psfi8!T=$^ zf0Y4qH5^$cF}G;eym6z)xyO(&hJ*^qnlo?ijG5kX?s1ciRMtH)l^;JLm7Ty?=ZHUt zf;kxL92hrAm9$)u&lh(|CodP{Ssw0UC=-%J=z|=uIo?Jdd>F%Bo-t9;8ZHT7h|M)T zz?mG#U#ob;cMJLRjCqQ7Q6NyrQTw(u{hXA(p1*ZJiK&OTMa|EeUa#n#8U0G*vH2`r zj?H?GBHnJJ4Ew<5%QdGlOZlkLuvEZ(>(e2}Sy^a1go1L{bI_Gzt8g2B!nh3bbM!(ob-)f_O z-A3PLqYtssx7+ALZS);B`Xx5{E*t$)8@u%CtPQbnPO`&ldbwb)a56b@Kgh|E2XJ=|_aH}nC3%UiB+|$| z#5r{&It2baaZVYDn82SR&Z#1?OW^-boKs(7yTBhK&M70YRp37&&Z#2NDDdwS=M<5s z7x-r4oEj3#1%5AaP6>$x0>7I$r_h8);6dV?8WLj!ekXBG2?>|L%ZYO;NE8YDcH*3B z6NbPS66X|~=>7_Xk>4aI z#5wdQwhH_b;vDi5jRL=rIEVT~y}$<$=MbM*F7O=U9NH5LfU`ZvelRz3aG4c3W;Gq} zUQlX|jg5WPYMxq$AS4H$m}(fGRJonQ?}G{&ym{)8ZxG}kBADD09*nN=1ccZX`X{@b z!y75LS3$@>F61XZMsUqnR;1f%{cxJq+LLQJ_E`tN3SR^R4az`4va8BDd{(bM{cD_B z16`wGe3{iWbs;1TEAn2r(27o-4SI602c&ofV!bba76^_z_V)S&x$N4ZZmhV(kI+nD9?(P!Ro@|Xvq$%;#I5lXs*?K*#N5EC{%Ns%pH!X zxnqjCd)TZ0{IMzyS*^#=)*V*Vh?eb=cph*Fwx4z`E@_d*7w^z2 z3@OkZkDV>>Gn+sUyzZ9u*!Lxyy{Pa3!Cqo-gI~djXah=ioOUiLLreyhbwRwOhtVmQ zI_(tV4ZUcEFt^0)O1bPouj`C79HwJv$&;cvHsKB*l+tiiR1}29?(QWQ*vKr;>|yxf zM?2zq1o6D1J8GU3JvC4E&`!h}f(UKVf`UGIi+=+>d-zWx+$X1)hn-uhMYz-bwjvE3 zk+K6;q$CEvI-~V@QS-nfpMNF&>=X~uA#~#oy3ZRzaGwYPAEZy(%;TcEPvYKY#+g?U zI59SA#+xr-5IUotJm-_<@yH%{g1%FM?qhbKCeCLK=d&?Kd&G=4fulqxn&LYV-dV2vmL=~pJMKGZh0G{n*E7hWPR@>M<&E{6JpvJ zd1xIptVlEG^EUHQ5tCnF#UQJ1AD!CF?c{0(sw#)ASmF@KaAK3J-=nPG_VoHSqG7D) z5Y%uvo8)IJStv9_)N#vy5R)h?Xpj7MC2@A-S3$bZ$lu;0B7a1hDr#mqjNsn*a;y^4 zBe=LEh9}c(hd-)MFb5h&jV4q;)o8N0-MQs**nykknPfFp=Ya4sp^r>a333ucrigYS z-)mf$WAr}2(;Mg0%Nw5FeV^-Tkex99hge?FhYmIG>8Xgp*K-S(3$!O6zP~aJz+e$B zs{l7Sr<`^!MF$q{EluyWwV(BliR|C(hu52(t1o6A1h|BGB4Yj&&soIWk(x(ZoLks= zQu7G#Qs7bZr|9f0>8K;Je^XJDxmz_&-ek=0Pc@m3qPIUfWdl31bIUBO0PG!}Qy4mX9LD2f z@WgOFicgiEKjxs?ausr>&Fm2#j-arNZTi5*CC1l3WnFtvX3rG!%(}x`mb@=NLE&mV zSA&)nDciwOYhDDanWptEMz|lq4x`~JfNe(jVk~BM!NzuL@*c}MgX08NJw@zb2@4|} z*nAFl3HIzgn6EJ8+DJ1=>ax_J2-qR2Fb+a3(37J&b`6xm{r_&^zbKHXnZXnx= zqMrHF$AjNoSdP7f6?qM9-NG=(-}*S26zjvp``A^rnK6tqn3Pb*Q*Wc6b+VuFhR3`P zqOza8OqXNWh0r?=6ZyYV2!=oLWL`%DqN9=kDr>^|ozG?dLrS{KYQ}VmZr6-%hfTwu zF-Z^|e41V>=RE~hWWxec4+$VJ2|GAccDv}Ccs!$JFIyG+pJ$D%&p(3vV#+t??^Fg~ z{zw=UR54|>#SfAc*}z>N+HFUq&5AalpFSiA!mWdmFX)U5qZ>@!qfFfeJ**lcn5{)A z&L{UoHq3&RsJQDnu`OON3y>3q5M>z;N&K0!iuoOiLEPmsUq>!KgyvG+K#U@c-U>k* zv|%j@1+HP?d_iA^Q7C^Fvg}T>;Vwov*^UUc*+>nJM;PkOMnRQA4x7*>2$k?AY@?!7 zu$l)q%PFXLTDKyE^u}x(-0r7}x?2_YC&`R-{+)vk&P$2Bqe~Uw|drgI2k1iytAg8ScTU0#>XU@j?C! zAAzXszpS$?K|~7Zc%tgzdJc;@m`51e;+KLhyhY7>Dy?fT65~cVKVA+-Y8}7awqVT| zT;Ss=x-VKxwCgrH$i;or`fkH{-v}(dFv~W@E_Cp5uWB=2a~QaHI54`rGYh*?=W}Mq zwm@vtVuyi&pO_<76)_J%=x`BrsdyYEbmMqngNRk@^LdEHiy|OBB3fl&VpNYDe4fM? z$!7iDmZ7BPg0nPV>}B3?7V}@h#C`(}hsp`-@w2S2?=2c5s4AL=Q!ILUFSFw;=FDE^ zJ1M4i`l?=L*ICT|z05-?rY`g32W(p5?cuYSzY|Om#hY;gWku{!e5zOTcuG?{^bnYO z(=}c037ao^1@&aR_ij-z$VbdKRr%-(83CCOE8K*M&En2=&xsaApuP_J3CC8T?16IT zLJ#aZY3qR=^gv9|4Y9Mk-~BjUgD$wD`wfK}x?dN%UzgD^6pt!jy|F>w1xNJt#7ENe zuICxI7ezh5f9ZsHE)epzcoXy3G{SI~m_|}kk0>cgy^fZ&Jl9xF@$0x?XnlVgHwXuN zBIY|$^Bpk@h@+r8rQz;2vy>aVL=|(|mPk18)Qcw5-@r1XT-69xj*@$O* z+65xG7>z}6UwB`Hr4zyJG#V~LaIc^Min~`9CrU|tM=bNq&h+?tsm$6y{O_cpE{(rV zoGo_GP)#!I7KclVqGnf1UBr5QoGq+Oc1ZKmgHvo5MmAc~34yp6Dp>H{Qc5`SW0~)g z$8=^71bZm6FCl6POL$l|#2T)k1u{M$LkEXUgfwv?p1vt8C?vm^K; z?a)#>q#Q}LrAI$22vE!rOWqhx(^!)MEtqkViR)!FLB(RKDb~+wjiYm--=OcLB2_P= zpu(&)Kb4i0eHYInmTO#1ay`wnmiuDiZwgErWJEigy81;-Y%0najvgU5M#bW=7_1SH zV{*5Y66yaMtB*=|auF~6oQ#G11L%^^;4~gY{k&fjSB=H#rr}C-$pq zi@Kkfgf78`V`Y20#nV!@n|k-#^uGMF8Z%PS9fRnyOLe~JW-K88Ib_6G!TuN?D2l|% zA_IEjd@dD)KQg+Me~$7^u>uIlE(Y2nX`=^C)?+mE6T{SrQNqpUsrxz2M;j=N{&~<6 z>rf#QKk+9~$KKP;BDg25gt6XdpB6{$SLr;;x|(H?)y8heSR(5|8HRzaUhm$@;as8!gma`Zm&XUv2@C4{6kS|ZqaxU4++pC~_i0#U$s_ERfOm{8@#r=UZIvbr0 zxmTUhn~;NO(cnT9Mk}4X9mlP7*oBf*I-jj3{UIP)b~( zBv_J(eBv!Si+~nuYpW+SNY2eXL_rQg`4^cxEh=Ae!JK*TE=XX(@~ zyw#&uU#w|ih4VXzqn>-#>>yqkQ}}-8#=THG@qJOh(xZtD!bl@ZK(VE*iB&dSwtVRk z=SDB6jMfsGLF5vbfB^ZCYe7id%4(y~RYayqWbiT~vn6sQaXIAq{NTI>gRFkz$x_dy za}oIe5dr?k9aC92P#p;Q!hyQl8(kwquA8a?6)OY&px0j;#{VCFV}LQwUlDEy`dyU` zF2!@Z-9||ua8>)p-VocNTtZN9Agt@q5SVV@B! zH)i91N?qYP7ypIS(#q*KZd?GN78%cSTj zg+`roT$($DZ+LLLWyJoFOm+a)zblt)fv#k-2l(h?$>e32&AI@&0@?P*WO6Cs6UURu zb%4qDlF99W-B|X$47d{O-s6B@VRw*+?zs1IM2H50rs8zuE*?jpOzm zfal{3svU4SU>uMMhqfX)?p|m())qM~y`W#gR@f#TFMf?(|Cvl)MTuEFPcn=VxF+C= zzmrVvf}An1X!gJ%-*le8rl8T7cFB}c<40Txs^P-56juOsrO(Fnfkhj0W)_}jHs$6% z??5Tiej`}@k0z5(1Eyqho=EYxgZwD?Vlk<0 zc)+R>t4kh>Teb?uqLzCe8o!t7IMV4z&it zfs05*IEZ7y8Iedmimx|GKINbO4bKoIH&{Up^wKImUsAYKWU5ygy`ERa8L8jMsDd=T z7)QQL+V4(@8j5V_^6_tuqCQ#_M@9l~QHt392weMt<)^@dijOzB0@r@4V*b}b>i6Rr zB~McQV77v{D##c0Osf^#sNlm2{#?Ol72K!bn+pC(!Os*t5C2-gbg_b0DL6^N*$UpO zV3mTa72K%c!wUXf!DkiRr=TnazmHoWb8}|SywO#RmG?ESu@l|n++$ti#*CXd=7zE3 zwPfFr-7(MjD4EG+|LKd)=cq3__4}gBeqChf-?3AUsTh|7y;k>cLBtLJ2al7QLL+vK zr06}y9H0(ku#wrmzeELG=*1l>4pIx-Sad}0#WNGd`l6R-p+5ur7s4O?od~fe z86;ib50ei64=~_JV$TPdz}1&OXR^=-syM7t?P5WWX_%tNu8_Fc(;|Ddpi7$l{`a7} z>;#ox1o|MuDSzY&n%*~-XR*T-+t@(IxU`u+PEWoO0+X<#|pgDF4DxG^LP{n{+_n?Ps3>dL6s z%q;S^WufyIWqq~#hAi}d67)i&SuGI6c^m5WZ&~E|JsS4sW%z9jlQ_2m{14J|uu3Ts z`_uLv>}vc*8qu$<5PC=FXHq%=g%TO|L%|kVdT3r`~6MoPeJcNE0o`k+(No zhu{9=O9tY{{b_8H$)=dW4+I48O@UQ@uRIJETD86!zxvQXx+;&7( z5FZrq@~9b~LP!(*t6y7~J0Bl1s0`G48$$j{em9|4b+Y&>LNAW54alz{q>+U6K0nal z_3}7eekZ~3RaSay>+m82?O^{pLYgi0>4h{JzwnSo@OtsR2(NkbEPNe8e&Zp7T^XwL zuEh5m&;%9YrGyAmNPgZSL$*8=QVwQ*5+Wry!9Y8Oy%lw}_*O-=eA9v>eA&?QMNCM2 z-65qj4za>pok~M~?jeoNPfcW4<2NgMWhSfdN%T@+ji06vtt-DAk;aoBwMZlJdlhK} zuV1{IP@kwsWAQ^8X@s*raFOETD--zGMheBZC|39yLZN`KRy8g^!O_cx-{SiiUPjO` z+=y1#R}L5sO1M(zYU{#&_lnvEcYP2YJs7^raHmM+4FP8%Rl zVf3!7gs#SvThnLE9Ub1^>hkqIH`D(`(C403>+5~5hMYf{xG!eeT3_!!zYQ8iw7%X) zYpD0l5D_WU{|ZpKN7edzU#j5@rKkPZavIiy&pofk^}bg_@|oyAWwigKG4aJ9+eGW@ zeYb|46hhMFYkgh+2zb<=K_Mo6V4$Gh_fwxs&rJOvfX0~UbBENtt>$mFT!Ms5nf>Ph zknQ>|rKh2I2dI?+6*~nV$~5{GV~3-?^ZCgeBFN9;48@3pa0D7&AGIof|=#N z4mx&!;u)fdKG)LUb;~SI;~KsT9qhhR`ug0fWh$-p_BV<3;#&U%80>rMF;DH&_4oaB z!|490^>jOY0$sLWN`F2Oq}VhGizP_iwZ4Xb0mrVd&lTcIKc3Q*8Lg*b4wUWs`aGe- zBjvT9+Kkpq;ZP<``>)siT^=cssqd2H)V)xZKPE+y8Lgl3_y}JqL7$ljGrG8-8`CAQ wRN#4z*&mnoUynbeY2T%Wyjyc5d?|ye%hj?xmdX$jtbX5;5_%t~ +#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; +}