12 Initialize(model_order, exponent_number, gate_beg, gate_end);
36 uint16_t uZeroLevel) {
44 float front_time_end) {
45 return std::ceil((3 * front_time_beg_03 - front_time_end) / 2.);
57 std::vector<float> a_f;
58 std::vector<float> a_b;
64 printf(
"LSerr %e, SLE roots forward ", rho_f);
65 for (uint8_t
i = 0;
i < a_f.size();
i++)
66 printf(
" %e ", a_f.at(
i));
68 printf(
"LSerr %e, SLE roots backward ", rho_b);
69 for (uint8_t
i = 0;
i < a_b.size();
i++)
70 printf(
" %e ", a_b.at(
i));
79 for (uint8_t
i = 0;
i < a_f.size();
i++)
80 a_arr[
i + 1] = a_f.at(
i);
94 printf(
"forward polinom roots ");
96 printf(
"%.5f%+.5fi ", zr[
i], zi[
i]);
99 printf(
"forward freqs ");
101 printf(
"%.5f ",
log(zr[
i]));
105 std::complex<float>* z_arr =
new std::complex<float>[
fExpNumber + 1];
107 if (std::isfinite(zr[
i]))
108 z_arr[
i + 1] = {zr[
i], zi[
i]};
123 std::vector<float>& a_f,
125 std::vector<float>& ) {
126 std::vector<int32_t> kiWfmSignal;
130 int n = kiWfmSignal.size();
146 for (
int sample_curr =
fModelOrder; sample_curr < n; sample_curr++)
147 Rkm_arr[
i - 1][j - 1] += (
float) (kiWfmSignal.at(sample_curr -
i)
148 * kiWfmSignal.at(sample_curr - j));
151 for (
int sample_curr =
fModelOrder; sample_curr < n; sample_curr++)
152 R0k_arr[
i - 1] -= (
float) (kiWfmSignal.at(sample_curr)
153 * kiWfmSignal.at(sample_curr -
i));
156 printf(
"system forward\n");
159 printf(
"%e ", Rkm_arr[
i][j]);
160 printf(
"%e\n", R0k_arr[
i]);
171 printf(
"SLE roots ");
173 printf(
" %e ", a[
i]);
189 std::vector<float>& a_f,
191 std::vector<float>& a_b) {
229 std::vector<int32_t> kiWfmSignal;
233 int n = kiWfmSignal.size();
236 printf(
"ERROR: Order too high; will make solution singular\n");
241 for (
int k = 1; k <= n - 2; k++)
242 r1 += std::pow(kiWfmSignal.at(k), 2);
244 float r2 = std::pow(kiWfmSignal.front(), 2);
245 float r3 = std::pow(kiWfmSignal.back(), 2);
249 r1 = 1. / (rho_b + r3);
251 std::vector<float> c,
d;
252 c.push_back(kiWfmSignal.back() * r1);
253 d.push_back(kiWfmSignal.front() * r1);
257 std::vector<float> ef, eb, ec, ed;
258 std::vector<float> coeffs;
261 for (
int v_iter = 0; v_iter < n; v_iter++) {
262 ef.push_back(kiWfmSignal.at(v_iter));
263 eb.push_back(kiWfmSignal.at(v_iter));
264 ec.push_back(c.at(0) * kiWfmSignal.at(v_iter));
265 ed.push_back(
d.at(0) * kiWfmSignal.at(v_iter));
271 if (rho_f <= 0 || rho_b <= 0) {
273 printf(
"PsdPronyFitter::ERROR: prediction squared error was less "
274 "than or equal to zero\n");
278 gam = 1. - ec.at(n - k);
279 del = 1. - ed.front();
280 if (gam <= 0 || gam > 1 || del <= 0 || del > 1) {
282 printf(
"PsdPronyFitter::ERROR: GAM or DEL gain factor not in "
283 "expected range 0 to 1\n");
288 std::vector<float> eff, ebb;
292 for (
int v_iter = 0; v_iter < n - 1; v_iter++) {
293 eff.at(v_iter) = ef.at(v_iter + 1);
294 ebb.at(v_iter) = eb.at(v_iter);
295 delta += eff.at(v_iter) * ebb.at(v_iter);
298 float k_f = -delta / rho_b;
299 float k_b = -delta / rho_f;
302 rho_f = rho_f * (1 - k_f * k_b);
303 rho_b = rho_b * (1 - k_f * k_b);
306 std::vector<float> temp_af = a_f;
307 std::reverse(std::begin(a_b), std::end(a_b));
308 for (uint8_t
i = 0;
i < a_f.size();
i++)
309 a_f.at(
i) += k_f * a_b.at(
i);
312 std::reverse(std::begin(a_b), std::end(a_b));
313 std::reverse(std::begin(temp_af), std::end(temp_af));
314 for (uint8_t
i = 0;
i < a_b.size();
i++)
315 a_b.at(
i) += k_b * temp_af.at(
i);
326 for (
int v_iter = 0; v_iter < n; v_iter++) {
327 eb.at(v_iter) = ebb.at(v_iter) + k_b * eff.at(v_iter);
328 ef.at(v_iter) = eff.at(v_iter) + k_f * ebb.at(v_iter);
332 coeffs.at(0) = ec.front();
333 coeffs.at(1) = coeffs.at(0) / del;
334 coeffs.at(2) = coeffs.at(0) / gam;
337 std::vector<float> temp_c = c;
338 for (uint8_t v_iter = 0; v_iter < c.size(); v_iter++) {
339 c.at(v_iter) += coeffs.at(1) *
d.at(v_iter);
340 d.at(v_iter) += coeffs.at(2) * temp_c.at(v_iter);
344 std::vector<float> temp_ec = ec;
345 for (
int v_iter = 0; v_iter < n; v_iter++) {
346 ec.at(v_iter) += coeffs.at(1) * ed.at(v_iter);
347 ed.at(v_iter) += coeffs.at(2) * temp_ec.at(v_iter);
350 std::vector<float> ecc, edd;
353 for (
int v_iter = 0; v_iter < n - 1; v_iter++) {
354 ecc.at(v_iter) = ec.at(v_iter + 1);
355 edd.at(v_iter) = ed.at(v_iter);
358 if (rho_f <= 0 || rho_b <= 0) {
360 printf(
"PsdPronyFitter::ERROR2: prediction squared error was less "
361 "than or equal to zero\n");
364 gam = 1 - ecc.at(n - k - 1);
365 del = 1 - edd.front();
366 if (gam <= 0 || gam > 1 || del <= 0 || del > 1) {
368 printf(
"PsdPronyFitter::ERROR2: GAM or DEL gain factor not in "
369 "expected range 0 to 1\n");
375 coeffs.at(0) = ef.front();
376 coeffs.at(1) = eb.at(n - k - 1);
377 coeffs.at(2) = coeffs.at(1) / rho_b;
378 coeffs.at(3) = coeffs.at(0) / rho_f;
379 coeffs.at(4) = coeffs.at(0) / del;
380 coeffs.at(5) = coeffs.at(1) / gam;
383 std::vector<float> temp_ab = a_b;
384 std::reverse(std::begin(temp_ab), std::end(temp_ab));
385 std::reverse(std::begin(c), std::end(c));
387 for (uint8_t
i = 0;
i < a_b.size();
i++)
388 a_b.at(
i) += coeffs.at(5) * c.at(
i);
389 std::reverse(std::begin(c), std::end(c));
391 for (uint8_t
i = 0;
i < c.size();
i++)
392 c.at(
i) += coeffs.at(2) * temp_ab.at(
i);
393 c.push_back(coeffs.at(2));
395 std::vector<float> temp_af2 = a_f;
396 for (uint8_t
i = 0;
i < a_f.size();
i++)
397 a_f.at(
i) += coeffs.at(4) *
d.at(
i);
399 for (uint8_t
i = 0;
i <
d.size();
i++)
400 d.at(
i) += coeffs.at(3) * temp_af2.at(
i);
401 d.insert(
d.begin(), coeffs.at(3));
404 rho_f = rho_f - coeffs.at(4) * coeffs.at(0);
405 rho_b = rho_b - coeffs.at(5) * coeffs.at(1);
407 if (rho_f <= 0 || rho_b <= 0) {
409 printf(
"PsdPronyFitter::ERROR3: prediction squared error was less "
410 "than or equal to zero\n");
415 for (
int v_iter = 0; v_iter < n; v_iter++) {
416 ec.at(v_iter) = ecc.at(v_iter) + coeffs.at(2) * eb.at(v_iter);
417 eb.at(v_iter) = eb.at(v_iter) + coeffs.at(5) * ecc.at(v_iter);
418 ed.at(v_iter) = edd.at(v_iter) + coeffs.at(3) * ef.at(v_iter);
419 ef.at(v_iter) = ef.at(v_iter) + coeffs.at(4) * edd.at(v_iter);
425 std::memcpy(
fz, z, (
fExpNumber + 1) *
sizeof(std::complex<float>));
429 std::complex<float>
z2) {
430 std::complex<float>* z_arr =
new std::complex<float>[
fExpNumber + 1];
447 std::complex<float>** Zik_arr =
new std::complex<float>*[
fExpNumber + 1];
449 Zik_arr[
i] =
new std::complex<float>[
fExpNumber + 1];
451 Zik_arr[
i][j] = {0., 0.};
454 std::complex<float>* Zyk_arr =
new std::complex<float>[
fExpNumber + 1];
456 Zyk_arr[
i] = {0., 0.};
459 const std::complex<float> unit = {1., 0.};
463 std::complex<float> temp = std::conj(
fz[
i]) *
fz[j];
464 if (std::abs(temp - unit) > 1e-3) {
465 Zik_arr[
i][j] =
static_cast<std::complex<float>
>(
466 (std::pow(temp,
static_cast<float>(samples_in_gate)) - unit)
469 Zik_arr[
i][j] =
static_cast<std::complex<float>
>(samples_in_gate);
474 std::complex<float>* z_power =
new std::complex<float>[
fExpNumber + 1];
481 Zyk_arr[
i] += (std::complex<float>) (std::conj(z_power[
i])
482 * (float)
fuWfm.at(sample_curr));
488 printf(
"\nampl calculation\n");
492 "%e%+ei ", std::real(Zik_arr[
i][j]), std::imag(Zik_arr[
i][j]));
494 " %e%+ei\n", std::real(Zyk_arr[
i]), std::imag(Zyk_arr[
i]));
501 printf(
"amplitudes\n%.0f%+.0fi ", std::real(
fh[0]), std::imag(
fh[0]));
503 printf(
"%e%+ei ", std::real(
fh[
i]), std::imag(
fh[
i]));
510 std::complex<float> fit_ampl_in_sample = {0., 0.};
512 for (
int sample_curr = 0; sample_curr <
fSampleTotal; sample_curr++) {
513 fit_ampl_in_sample = {0., 0.};
516 fit_ampl_in_sample +=
fh[
i] * z_power[
i];
519 fuFitWfm.at(sample_curr) = (uint16_t) std::real(fit_ampl_in_sample);
525 printf(
"waveform:\n");
526 for (uint8_t
i = 0;
i <
fuWfm.size();
i++)
527 printf(
"%u ",
fuWfm.at(
i));
529 printf(
"\nfit waveform:\n");
548 for (
int sample_curr = gate_beg; sample_curr <= gate_end; sample_curr++)
551 if (std::isfinite(integral))
return integral;
556 if (std::isfinite(
fuFitWfm.at(sample_number)))
562 std::complex<float> amplitude = {0., 0.};
568 if (std::isfinite(std::real(amplitude)))
return std::real(amplitude);
582 if (std::real(
fh[2] /
fh[1]) < 0)
594 if (first_sample < last_sample)
598 float result_sample = 0.;
599 int sample_to_check = first_sample;
600 float amplitude = 0.;
601 float amplitude_prev =
GetFitValue(sample_to_check - step);
602 while ((first_sample - sample_to_check) * (last_sample - sample_to_check)
605 if ((level - amplitude) * (level - amplitude_prev) <= 0) {
608 sample_to_check - step,
611 return result_sample;
613 amplitude_prev = amplitude;
614 sample_to_check += step;
624 float result_sample = 0.;
625 float sample_to_check = (float) first_sample;
626 std::complex<float> amplitude = {0., 0.};
627 std::complex<float> amplitude_prev =
GetFitValue(sample_to_check - step);
628 while ((first_sample - sample_to_check) * (last_sample - sample_to_check)
631 if ((level - std::real(amplitude)) * (level - std::real(amplitude_prev))
633 if (amplitude != amplitude_prev)
635 std::real(amplitude),
636 sample_to_check - step,
637 std::real(amplitude_prev),
639 return result_sample;
641 amplitude_prev = amplitude;
642 sample_to_check += step;
653 return (X1 * Y0 - X1 * Y2 - X2 * Y0 + X2 * Y1) / (Y1 - Y2);
660 int m = gate_end - gate_beg + 1;
662 if (
m <= params_number)
return 999;
664 for (
int sample_curr = gate_beg; sample_curr <= gate_end; sample_curr++)
665 average +=
fuWfm.at(sample_curr);
668 for (
int sample_curr = gate_beg; sample_curr <= gate_end; sample_curr++) {
669 RSS += std::pow(
fuFitWfm.at(sample_curr) -
fuWfm.at(sample_curr), 2);
670 TSS += std::pow(
fuWfm.at(sample_curr) - average, 2);
672 if (TSS == 0)
return 999;
677 float R2_adj = R2 * (
m - 1) / (
m - params_number);
687 int freedom_counter = 0;
688 int regions_number = 10;
690 if (amplitude_max == 0)
return 999;
692 int* probability_exp =
new int[regions_number];
693 int* probability_theor =
new int[regions_number];
694 float* amplitude_regions =
new float[regions_number + 1];
695 amplitude_regions[0] = 0.;
696 for (
int i = 0;
i < regions_number;
i++) {
697 probability_exp[
i] = 0;
698 probability_theor[
i] = 0;
699 amplitude_regions[
i + 1] = (
i + 1) * amplitude_max / regions_number;
702 for (
int sample_curr = gate_beg; sample_curr <= gate_end; sample_curr++) {
703 for (
int i = 0;
i < regions_number;
i++) {
705 > amplitude_regions[
i])
707 <= amplitude_regions[
i + 1]))
708 probability_exp[
i]++;
710 > amplitude_regions[
i])
712 <= amplitude_regions[
i + 1]))
713 probability_theor[
i]++;
717 for (
int i = 0;
i < regions_number;
i++) {
718 if (probability_exp[
i] > 0) {
719 chi2 += std::pow(probability_exp[
i] - probability_theor[
i], 2.)
720 / (probability_exp[
i]);
725 if (freedom_counter > 0) chi2 /= freedom_counter;
726 delete[] probability_exp;
727 delete[] probability_theor;
728 delete[] amplitude_regions;
763 int best_signal_begin = 0;
764 bool IsReasonableRoot;
765 bool IsGoodFit =
false;
766 int good_fit_counter = 0;
768 for (
int signal_begin = first_sample; signal_begin <= last_sample;
772 IsReasonableRoot =
true;
774 IsReasonableRoot = IsReasonableRoot && (std::abs(
fz[j + 1]) > 1e-6)
775 && (std::abs(
fz[j + 1]) < 1e1);
780 printf(
"good fit candidate at signal begin %i\n", signal_begin);
784 if (good_fit_counter == 1) {
786 best_signal_begin = signal_begin;
790 best_signal_begin = signal_begin;
795 return best_signal_begin;
800 int best_signal_begin = first_sample;
802 for (
int signal_begin = first_sample; signal_begin <= last_sample;
807 if (signal_begin == first_sample) best_R2 = R2;
810 best_signal_begin = signal_begin;
814 return best_signal_begin;
818 bool solvable =
true;
821 float** a =
new float*[n];
822 for (
int i = 0;
i < n;
i++) {
823 a[
i] =
new float[n + 1];
824 for (
int j = 0; j < n + 1; j++)
828 for (
int i = 0;
i < n;
i++) {
829 for (
int j = 0; j < n; j++)
834 for (
int i = 0;
i < n;
i++) {
835 maxEl = std::abs(a[
i][
i]);
837 for (
int k =
i + 1; k < n; k++)
838 if (abs(a[k][
i]) > maxEl) {
839 maxEl = std::abs(a[k][
i]);
845 if (
fIsDebug) printf(
"SLE has not solution\n");
848 for (
int k =
i; k < n + 1; k++) {
850 a[maxRow][k] = a[
i][k];
854 for (
int k =
i + 1; k < n; k++) {
855 c = -a[k][
i] / a[
i][
i];
856 for (
int j =
i; j < n + 1; j++) {
860 a[k][j] += c * a[
i][j];
865 for (
int i = n - 1;
i >= 0;
i--) {
866 x[
i] = a[
i][n] / a[
i][
i];
867 for (
int k =
i - 1; k >= 0; k--)
868 a[k][n] -= a[k][
i] *
x[
i];
872 for (
int i = n - 1;
i >= 0;
i--)
876 for (
int i = 0;
i < n;
i++)
882 std::complex<float>** r,
883 std::complex<float>* b,
885 bool solvable =
true;
888 std::complex<float> tmp, c;
889 std::complex<float>** a =
new std::complex<float>*[n];
890 for (
int i = 0;
i < n;
i++) {
891 a[
i] =
new std::complex<float>[n + 1];
892 for (
int j = 0; j < n + 1; j++)
896 for (
int i = 0;
i < n;
i++) {
897 for (
int j = 0; j < n; j++)
902 for (
int i = 0;
i < n;
i++) {
903 maxEl = std::abs(a[
i][
i]);
905 for (
int k =
i + 1; k < n; k++)
906 if (std::abs(a[k][
i]) > maxEl) {
907 maxEl = std::abs(a[k][
i]);
913 if (
fIsDebug) printf(
"PsdPronyFitter:: SLE has not solution\n");
916 for (
int k =
i; k < n + 1; k++) {
918 a[maxRow][k] = a[
i][k];
922 for (
int k =
i + 1; k < n; k++) {
923 c = -a[k][
i] / a[
i][
i];
924 for (
int j =
i; j < n + 1; j++) {
928 a[k][j] += c * a[
i][j];
933 for (
int i = n - 1;
i >= 0;
i--) {
934 x[
i] = a[
i][n] / a[
i][
i];
935 for (
int k =
i - 1; k >= 0; k--)
936 a[k][n] -= a[k][
i] *
x[
i];
940 for (
int i = n - 1;
i >= 0;
i--)
944 for (
int i = 0;
i < n;
i++)
951 float** u =
new float*[n];
952 for (
int i = 0;
i < n;
i++) {
954 for (
int j = 0; j < n; j++)
958 float*
y =
new float[n];
959 for (
int i = 0;
i < n;
i++)
962 for (
int i = 0;
i < n;
i++) {
964 for (
int k = 0; k <
i; k++)
965 temp = temp + u[k][
i] * u[k][
i];
967 for (
int j =
i; j < n; j++) {
969 for (
int k = 0; k <
i; k++)
970 temp = temp + u[k][
i] * u[k][j];
971 u[
i][j] = (a[
i][j] - temp) / u[
i][
i];
975 for (
int i = 0;
i < n;
i++) {
977 for (
int k = 0; k <
i; k++)
978 temp = temp + u[k][
i] *
y[k];
979 y[
i] = (b[
i] - temp) / u[
i][
i];
982 for (
int i = n - 1;
i >= 0;
i--) {
984 for (
int k =
i + 1; k < n; k++)
985 temp = temp + u[
i][k] *
x[k];
986 x[
i] = (
y[
i] - temp) / u[
i][
i];
989 for (
int i = 0;
i < n;
i++)