22 #include <ecl/digitization/algorithms.h>
23 #include <ecl/digitization/WrapArray2D.h>
29 int myPow(
int x,
int p)
33 return x * myPow(x, p - 1);
36 short int bit_val(
double f,
int idd,
short int bit,
int a,
int b,
int lim)
47 const int R(myPow(2, i));
48 const int val((
int)(a * R * f / idd / b + R + 0.5) - R);
50 if ((val < lim && val > -lim)) {
57 }
while (i > 8 && flag == 0);
60 if (i > bit) {i = bit;}
68 if (k > -1 && k < 8736) {
70 else if (k < 288)R = 64;
71 else if (k < 864)R = 96;
72 else if (k < 8064)R = 144;
73 else if (k < 8544)R = 96;
83 if (k > -1 && k < 8736) {
85 else if (k < 288)R = 96;
86 else if (k < 864)R = 288;
87 else if (k < 8064)R = 864;
88 else if (k < 8544)R = 8064;
98 if (k > -1 && k < 8736) {
100 else if (k < 288)R = 2;
101 else if (k < 864)R = 5;
102 else if (k < 8064)R = 11;
103 else if (k < 8544)R = 61;
116 if (k > -1 && k < 8736) {
117 const int l = Rest(k);
118 const int b = PhiC(k);
119 R = (k - l) / b + Rest1(k);
131 void bit(
int num,
char* inputPath)
144 double g1g1[192], gg[192], gg1[192];
146 double gg2[192], g1g2[192], g2g2[192];
147 double dgg1[192], dgg2[192];
171 memset(g1g1, 0,
sizeof(g1g1));
172 memset(gg, 0,
sizeof(gg));
173 memset(gg1, 0,
sizeof(gg1));
175 memset(sg1, 0,
sizeof(sg1));
176 memset(sg, 0,
sizeof(sg));
177 memset(sg2, 0,
sizeof(sg2));
178 memset(gg2, 0,
sizeof(gg2));
179 memset(g1g2, 0,
sizeof(g1g2));
180 memset(g2g2, 0,
sizeof(g2g2));
181 memset(dgg1, 0,
sizeof(dgg1));
182 memset(dgg2, 0,
sizeof(dgg2));
191 double adt, ddt, tc1;
192 double tmd, tpd, ssssj, ssssj1, ssssi, ssssi1;
198 int mxf, mxf1, mxfg31, mxfg32, mxfg33, mxfg41, mxfg43;
231 const int mapsize = 217;
233 short int bfg31[mapsize];
235 short int bfg32[mapsize];
236 short int bfg33[mapsize];
237 short int bfg41[mapsize];
239 short int bfg43[mapsize];
240 short int bf[mapsize];
241 short int bf1[mapsize];
242 short int bad[mapsize];
246 for (j = 0; j < mapsize; j++) {
295 const double ndt = 96.;
297 const int endc = 2 * ndt;
300 const double dt0 = adt * dt;
301 double t0 = -dt0 * ndt;
309 const double del = 0.;
310 const double ts0 = 0.5;
313 for (j1 = 0; j1 < endc; j1++) {
315 double t = t0 - dt - del;
319 for (j = 0; j < 16; j++) {
325 f[j1][j] = ShaperDSP(t / 0.881944444);
334 svp = ShaperDSP(tpd / 0.881944444);
335 svm = ShaperDSP(tmd / 0.881944444);
337 f1[j1][j] = (svp - svm) / 2. / ddt;
342 f1[j1][j] = ShaperDSP(tpd / 0.881944444) / (ddt + tc1);
348 if (fabs(f1[j1][j]) > 100. || (j == 2 && j1 == 156)) {
374 for (ChN = 0; ChN < mapsize; ChN++) {
375 if (ChN % (mapsize / 50) == 0) {printf(
"CnN=%d badN=%d \n", ChN, badN);}
387 if ((BMcoIN = fopen(
filename.c_str(),
"rb")) == NULL) {
388 printf(
" file %s is absent \n",
filename.c_str());
393 if (fread(ss1,
sizeof(
double), 256, BMcoIN) != 256) {
394 printf(
"Error writing file %s \n",
filename.c_str());
409 for (j1 = 0; j1 < endc; j1++) {
417 for (j = 0; j < 16; j++) {
425 for (i = 0; i < 16; i++) {
427 sg[j][j1] = sg[j][j1] + ss1[j][i] * f[j1][i];
428 sg1[j][j1] = sg1[j][j1] + ss1[j][i] * f1[j1][i];
430 sg2[j][j1] = sg2[j][j1] + ss1[j][i];
435 gg[j1] = gg[j1] + ss1[j][i] * ssssj * ssssi;
439 gg1[j1] = gg1[j1] + ss1[j][i] * ssssj * ssssi1;
440 g1g1[j1] = g1g1[j1] + ss1[j][i] * ssssi1 * ssssj1;
442 gg2[j1] = gg2[j1] + ss1[j][i] * ssssj;
443 g1g2[j1] = g1g2[j1] + ss1[j][i] * ssssj1;
444 g2g2[j1] = g2g2[j1] + ss1[j][i];
453 dgg1[j1] = gg[j1] * g2g2[j1] - gg2[j1] * gg2[j1];
454 dgg2[j1] = gg[j1] * g1g1[j1] * g2g2[j1] - gg1[j1] * gg1[j1] * g2g2[j1] + 2 * gg1[j1] * g1g2[j1] * gg2[j1] - gg2[j1] * gg2[j1] *
455 g1g1[j1] - g1g2[j1] * g1g2[j1] * gg[j1];
458 for (i = 0; i < 16; i++) {
461 fg31[j1][i] = ((g1g1[j1] * g2g2[j1] - g1g2[j1] * g1g2[j1]) * sg[i][j1] + (g1g2[j1] * gg2[j1] - gg1[j1] * g2g2[j1]) * sg1[i][j1] +
462 (gg1[j1] * g1g2[j1] - g1g1[j1] * gg2[j1]) * sg2[i][j1]) / dgg2[j1];
465 fg32[j1][i] = ((g1g2[j1] * gg2[j1] - gg1[j1] * g2g2[j1]) * sg[i][j1] + (gg[j1] * g2g2[j1] - gg2[j1] * gg2[j1]) * sg1[i][j1] +
466 (gg1[j1] * gg2[j1] - gg[j1] * g1g2[j1]) * sg2[i][j1]) / dgg2[j1];
468 fg33[j1][i] = ((gg1[j1] * g1g2[j1] - g1g1[j1] * gg2[j1]) * sg[i][j1] + (gg1[j1] * gg2[j1] - gg[j1] * g1g2[j1]) * sg1[i][j1] +
469 (gg[j1] * g1g1[j1] - gg1[j1] * gg1[j1]) * sg2[i][j1]) / dgg2[j1];
478 fg41[j1][i] = (g2g2[j1] * sg[i][j1] - gg2[j1] * sg2[i][j1]) / dgg1[j1];
479 fg43[j1][i] = (gg[j1] * sg2[i][j1] - gg2[j1] * sg[i][j1]) / dgg1[j1];
497 constexpr
int n16 = 16;
500 const int ipardsp13 = 14 + 14 * 256;
501 const int ipardsp14 = 0 * 256 + 17;
503 const int ibb = ipardsp13 / 256;
504 const int iaa = ipardsp13 - ibb * 256;
505 idd = ipardsp14 / 256;
519 for (i = 0; i < 16; i++) {
527 for (k = 0; k < 192; k++) {
530 bf[ChN] = bit_val(val_f, idd, bf[ChN], 1, 1, i16);
531 if (bf[ChN] < bitf) {bitf = bf[ChN];}
534 bf1[ChN] = bit_val(val_f, idd, bf1[ChN], 4, 3, i16);
535 if (bf1[ChN] < bitf1) {bitf1 = bf1[ChN];}
537 bfg31[ChN] = bit_val(val_f, idd, bfg31[ChN], 1, 1, i16);
538 if (bfg31[ChN] < bitfg31) {bitfg31 = bfg31[ChN];}
540 bfg32[ChN] = bit_val(val_f, idd, bfg32[ChN], 3, 4, i16);
541 if (bfg32[ChN] < bitfg32) {bitfg32 = bfg32[ChN];}
543 bfg33[ChN] = bit_val(val_f, idd, bfg33[ChN], 1, 1, i16);
544 if (bfg33[ChN] < bitfg33) {bitfg33 = bfg33[ChN];}
547 bfg41[ChN] = bit_val(val_f, idd, bfg41[ChN], 1, 1, i16);
548 if (bfg41[ChN] < bitfg41) {bitfg41 = bfg41[ChN];}
550 bfg43[ChN] = bit_val(val_f, idd, bfg43[ChN], 1, 1, i16);
551 if (bfg43[ChN] < bitfg43) {bitfg43 = bfg43[ChN];}
556 if (fabs(f[k][i] / fdd) > xf) {xf = fabs(f[k][i] / fdd);}
557 if (fabs(f1[k][i] / fdd) > xf1) {xf1 = fabs(f1[k][i] / fdd);}
558 if (fabs(fg31[k][i] / fdd) > xfg31) {xfg31 = fabs(fg31[k][i] / fdd);}
559 if (fabs(fg32[k][i] / fdd) > xfg32) {xfg32 = fabs(fg32[k][i] / fdd);}
560 if (fabs(fg33[k][i] / fdd) > xfg33) {xfg33 = fabs(fg33[k][i] / fdd);}
562 if (fabs(fg41[k][i] / fdd) > xfg41) {xfg41 = fabs(fg41[k][i] / fdd);}
563 if (fabs(fg43[k][i] / fdd) > xfg43) {xfg43 = fabs(fg43[k][i] / fdd);}
573 if (bf[ChN] < cf || bf1[ChN] < cf1 || bfg31[ChN] < cf31 || bfg32[ChN] < cf32 || bfg33[ChN] < cf33 || bfg41[ChN] < cf41
574 || bfg43[ChN] < cf43) {
577 if (bf[ChN] < cf) {bad[ChN] = bad[ChN] + 1;}
578 if (bf1[ChN] < cf1) {bad[ChN] = bad[ChN] + 2;}
579 if (bfg31[ChN] < cf31) {bad[ChN] = bad[ChN] + 4;}
580 if (bfg32[ChN] < cf32) {bad[ChN] = bad[ChN] + 8;}
581 if (bfg33[ChN] < cf33) {bad[ChN] = bad[ChN] + 16;}
582 if (bfg41[ChN] < cf41) {bad[ChN] = bad[ChN] + 32;}
583 if (bfg43[ChN] < cf43) {bad[ChN] = bad[ChN] + 64;}
591 printf(
"bxf=%d bxf1=%d bxfg31=%d bxfg32=%d bxfg33=%d bxfg41=%d bxfg43=%d \n", bitf, bitf1, bitfg31, bitfg32, bitfg33, bitfg41,
593 printf(
"xf=%lf xf1=%lf xfg31=%lf xfg32=%lf xfg33=%lf xfg41=%lf xfg43=%lf \n", xf, xf1, xfg31, xfg32, xfg33, xfg41, xfg43);
596 mxf = (int)(xf * ia / idd + ia + 0.5) - ia;
598 ia = myPow(2, bitf1);
599 mxf1 = (int)(4 * xf1 * ia / idd / 3 + ia + 0.5) - ia;
600 ia = myPow(2, bitfg31);
601 mxfg31 = (int)(xfg31 * ia / idd + ia + 0.5) - ia;
602 ia = myPow(2, bitfg32);
603 mxfg32 = (int)(3 * xfg32 * ia / idd / 4 + ia + 0.5) - ia;
604 ia = myPow(2, bitfg33);
605 mxfg33 = (int)(xfg33 * ia / idd + ia + 0.5) - ia;
606 ia = myPow(2, bitfg41);
607 mxfg41 = (int)(xfg41 * ia / idd + ia + 0.5) - ia;
608 ia = myPow(2, bitfg43);
609 mxfg43 = (int)(xfg43 * ia / idd + ia + 0.5) - ia;
611 printf(
"mf=%d mf1=%d mfg31=%d mfg32=%d mfg33=%d mfg41=%d mfg43=%d \n", mxf, mxf1, mxfg31, mxfg32, mxfg33, mxfg41, mxfg43);
614 printf(
"Rf=%lf Rf1=%lf Rfg31=%lf Rfg32=%lf Rfg33=%lf Rfg41=%lf Rfg43=%lf \n", (
double)mxf / (
double)i16, (
double)mxf1 / (
double)i16,
615 (
double)mxfg31 / (
double)i16, (
double)mxfg32 / (
double)i16, (
double)mxfg33 / (
double)i16, (
double)mxfg41 / (
double)i16,
616 (
double)mxfg43 / (
double)i16);
624 mxf = (int)(xf * ia / idd + ia + 0.5) - ia;
626 printf(
"f at limit %lf\n", (
double)mxf / (
double)i16);
627 if ((
double)mxf / (
double)i16 > 1) {printf(
"wronf bit for f\n");}
630 mxf = (int)(xf * ia / idd + ia + 0.5) - ia;
632 printf(
"f %lf \n", (
double)mxf / (
double)i16);
633 if ((
double)mxf / (
double)i16 > 1) {printf(
"wronf bit for f\n");}
636 ia = myPow(2, bitf + 1);
637 mxf = (int)(xf * ia / idd + ia + 0.5) - ia;
639 printf(
"f+1 %lf\n", (
double)mxf / (
double)i16);
641 ia = myPow(2, bitf - 1);
642 mxf = (int)(xf * ia / idd + ia + 0.5) - ia;
644 printf(
"f-1 %lf\n", (
double)mxf / (
double)i16);
653 ia = myPow(2, bitf1);
654 mxf1 = (int)(4 * xf1 * ia / idd / 3 + ia + 0.5) - ia;
656 printf(
"f at limit %lf\n", (
double)mxf1 / (
double)i16);
657 if ((
double)mxf / (
double)i16 > 1) {printf(
"wronf bit for f1\n");}
659 ia = myPow(2, bitf1);
660 mxf1 = (int)(4 * xf1 * ia / idd / 3 + ia + 0.5) - ia;
662 printf(
"f1 %lf\n", (
double)mxf1 / (
double)i16);
663 if ((
double)mxf1 / (
double)i16 > 1) {printf(
"wronf bit for f1\n");}
666 ia = myPow(2, bitf1 + 1);
667 mxf1 = (int)(4 * xf1 * ia / idd / 3 + ia + 0.5) - ia;
669 printf(
"f1+1 %lf \n", (
double)mxf1 / (
double)i16);
671 ia = myPow(2, bitf1 - 1);
672 printf(
"ia=%d i16=%d \n", ia, i16);
674 mxf1 = (int)(4 * xf1 * ia / idd / 3 + ia + 0.5) - ia;
676 printf(
"f1-1 %lf\n", (
double)mxf1 / (
double)i16);
685 ia = myPow(2, bitfg31);
686 mxfg31 = (int)(xfg31 * ia / idd + ia + 0.5) - ia;
688 printf(
"fg31 at limit %lf\n", (
double)mxfg31 / (
double)i16);
689 if ((
double)mxfg31 / (
double)i16 > 1) {printf(
"wronf bit for fg31\n");}
691 ia = myPow(2, bitfg31);
692 mxfg31 = (int)(xfg31 * ia / idd + ia + 0.5) - ia;
694 printf(
"fg31 %lf\n", (
double)mxfg31 / (
double)i16);
695 if ((
double)mxfg31 / (
double)i16 > 1) {printf(
"wronf bit for fg31\n");}
698 ia = myPow(2, bitfg31 + 1);
699 mxfg31 = (int)(xfg31 * ia / idd + ia + 0.5) - ia;
701 printf(
"fg31+1 %lf \n", (
double)mxfg31 / (
double)i16);
703 ia = myPow(2, bitfg31 - 1);
704 mxfg31 = (int)(xfg31 * ia / idd + ia + 0.5) - ia;
706 printf(
"fg31-1 %lf\n", (
double)mxfg31 / (
double)i16);
715 ia = myPow(2, bitfg32);
716 mxfg31 = (int)(3 * xfg32 * ia / idd / 4 + ia + 0.5) - ia;
718 printf(
"fg32 at limit %lf\n", (
double)mxfg32 / (
double)i16);
719 if ((
double)mxfg32 / (
double)i16 > 1) {printf(
"wronf bit for f32\n");}
721 ia = myPow(2, bitfg32);
722 mxfg32 = (int)(3 * xfg32 * ia / idd / 4 + ia + 0.5) - ia;
724 printf(
"fg32 %lf\n", (
double)mxfg32 / (
double)i16);
725 if ((
double)mxfg32 / (
double)i16 > 1) {printf(
"wronf bit for f\n");}
728 ia = myPow(2, bitfg32 + 1);
729 mxfg31 = (int)(3 * xfg32 * ia / idd / 4 + ia + 0.5) - ia;
731 printf(
"fg32+1 %lf \n", (
double)mxfg31 / (
double)i16);
733 ia = myPow(2, bitfg31 - 1);
734 mxfg31 = (int)(3 * xfg32 * ia / idd / 4 + ia + 0.5) - ia;
736 printf(
"fg32-1 %lf\n", (
double)mxfg31 / (
double)i16);
747 ia = myPow(2, bitfg33);
748 mxfg33 = (int)(xfg33 * ia / idd + ia + 0.5) - ia;
750 printf(
"fg33 at limit %lf\n", (
double)mxfg33 / (
double)i16);
751 if ((
double)mxfg33 / (
double)i16 > 1) {printf(
"wronf bit for fg33\n");}
753 ia = myPow(2, bitfg33);
754 mxfg33 = (int)(xfg33 * ia / idd + ia + 0.5) - ia;
756 printf(
"fg33 %lf\n", (
double)mxfg33 / (
double)i16);
757 if ((
double)mxfg33 / (
double)i16 > 1) {printf(
"wronf bit for fg33\n");}
760 ia = myPow(2, bitfg33 + 1);
761 mxfg33 = (int)(xfg33 * ia / idd + ia + 0.5) - ia;
763 printf(
"fg33+1 %lf \n", (
double)mxfg33 / (
double)i16);
765 ia = myPow(2, bitfg33 - 1);
766 mxfg33 = (int)(xfg33 * ia / idd + ia + 0.5) - ia;
768 printf(
"fg33-1 %lf\n", (
double)mxfg33 / (
double)i16);
778 ia = myPow(2, bitfg41);
779 mxfg41 = (int)(xfg41 * ia / idd + ia + 0.5) - ia;
781 printf(
"fg41 at limit %lf\n", (
double)mxfg41 / (
double)i16);
782 if ((
double)mxfg41 / (
double)i16 > 1) {printf(
"wronf bit for fg41\n");}
784 ia = myPow(2, bitfg41);
785 mxfg41 = (int)(xfg41 * ia / idd + ia + 0.5) - ia;
787 printf(
"fg41 %lf\n", (
double)mxfg41 / (
double)i16);
788 if ((
double)mxfg41 / (
double)i16 > 1) {printf(
"wronf bit for fg41\n");}
791 ia = myPow(2, bitfg41 + 1);
792 mxfg41 = (int)(xfg41 * ia / idd + ia + 0.5) - ia;
794 printf(
"fg41+1 %lf \n", (
double)mxfg41 / (
double)i16);
796 ia = myPow(2, bitfg41 - 1);
797 mxfg41 = (int)(xfg41 * ia / idd + ia + 0.5) - ia;
799 printf(
"fg41-1 %lf\n", (
double)mxfg41 / (
double)i16);
808 ia = myPow(2, bitfg43);
809 mxfg43 = (int)(xfg43 * ia / idd + ia + 0.5) - ia;
811 printf(
"fg43 at limit %lf\n", (
double)mxfg43 / (
double)i16);
812 if ((
double)mxfg43 / (
double)i16 > 1) {printf(
"wronf bit for fg33\n");}
814 ia = myPow(2, bitfg43);
815 mxfg43 = (int)(xfg43 * ia / idd + ia + 0.5) - ia;
817 printf(
"fg43 %lf\n", (
double)mxfg43 / (
double)i16);
818 if ((
double)mxfg33 / (
double)i16 > 1) {printf(
"wronf bit for fg33\n");}
821 ia = myPow(2, bitfg43 + 1);
822 mxfg43 = (int)(xfg43 * ia / idd + ia + 0.5) - ia;
824 printf(
"fg43+1 %lf \n", (
double)mxfg43 / (
double)i16);
826 ia = myPow(2, bitfg43 - 1);
827 mxfg43 = (int)(xfg43 * ia / idd + ia + 0.5) - ia;
829 printf(
"fg43-1 %lf\n", (
double)mxfg43 / (
double)i16);
841 if (bitf < 13) {printf(
"problem bitf=%d \n", bitf);}
842 if (bitf1 < 13) {printf(
"problem bitf1=%d \n", bitf1);}
843 if (bitfg31 < 13) {printf(
"problem bitfg31=%d \n", bitfg31);}
844 if (bitfg32 < 13) {printf(
"problem bitfg32=%d \n", bitfg32);}
845 if (bitfg33 < 16) {printf(
"problem bitfg33=%d \n", bitfg33);}
847 if (bitfg41 < 13) {printf(
"problem bitfg41=%d \n", bitfg41);}
848 if (bitfg43 < 16) {printf(
"problem bitfg43=%d \n", bitfg43);}
853 printf(
"bxf=%d bxf1=%d bxfg31=%d bxfg32=%d bxfg33=%d bxfg41=%d bxfg43=%d \n", bitf, bitf1, bitfg31, bitfg32, bitfg33, bitfg41,
856 printf(
"badN=%d \n", badN);
875 for (i = 0; i < 69; i++) {
896 for (j = 0; j < mapsize; j++) {
898 for (i = 10; i < 20; i++) {
899 if (bf[j] == i) {BF[i] = BF[i] + 1;}
900 if (bf1[j] == i) {BF1[i] = BF1[i] + 1;}
901 if (bfg31[j] == i) {BFG31[i] = BFG31[i] + 1;}
902 if (bfg32[j] == i) {BFG32[i] = BFG32[i] + 1;}
903 if (bfg33[j] == i) {BFG33[i] = BFG33[i] + 1;}
904 if (bfg41[j] == i) {BFG41[i] = BFG41[i] + 1;}
905 if (bfg43[j] == i) {BFG43[i] = BFG43[i] + 1;}
907 if (bf[j] < cf || bf1[j] < cf1 || bfg31[j] < cf31 || bfg32[j] < cf32 || bfg33[j] < cf33 || bfg41[j] < cf41 || bfg43[j] < cf43) {
908 if (bfg31[j] < cf31) {
909 RFG31[ThetR(j)] = RFG31[ThetR(j)] + 1;
911 if (bfg32[j] < cf32) {
912 RFG32[ThetR(j)] = RFG32[ThetR(j)] + 1;
914 if (bfg33[j] < cf33) {
915 RFG33[ThetR(j)] = RFG33[ThetR(j)] + 1;
917 if (bfg43[j] < cf43) {
918 RFG43[ThetR(j)] = RFG43[ThetR(j)] + 1;
924 for (i = 10; i < 20; i++) {
925 if (BF[i] > 0 || BF1[i] > 0 || BFG31[i] > 0 || BFG32[i] > 0 || BFG33[i] > 0 || BFG41[i] > 0 || BFG43[i] > 0) {
926 printf(
"bit=%d BF=%d BF1=%d BFG31=%d BFG32=%d BFG33=%d BFG41=%d BFG43=%d\n", i, BF[i], BF1[i], BFG31[i], BFG32[i], BFG33[i],
933 sprintf(BMin,
"bitst%d.dat", num);
935 if ((BMcoIN = fopen(BMin,
"w")) == NULL) {
936 printf(
" file %s is absent \n", BMin);
940 for (i = 0; i < mapsize; i++) {
941 fprintf(BMcoIN,
"%d %d %d %d %d %d %d %d \n", i, bf[i], bf1[i], bfg31[i], bfg32[i], bfg33[i], bfg41[i], bfg43[i]);
948 int main(
int argc,
char** argv)
950 assert(argc == 1 || argc == 3);
952 cout <<
"Usage:" << endl;
953 cout << argv[0] <<
" <type> <covmat_path>" << endl;
954 cout <<
"type is a numerical integer value to define the source of the calibration" << endl;
955 cout <<
"covmat_path is the path to look for the directory corr<type> that contains the covariance matrix files Binmcorr<type>_L.dat"
957 cout <<
" An output file with the translation to integers in written in the work directory" << endl;
958 }
else bit(atoi(argv[1]), argv[2]);