From 7cb15bbd95ab08f71efaac1ba296ead0a4152996 Mon Sep 17 00:00:00 2001 From: f4exb Date: Mon, 29 Jul 2024 20:12:44 +0200 Subject: [PATCH] WDSP: EMNR rework --- wdsp/calculus.cpp | 4 +- wdsp/calculus.hpp | 6 +- wdsp/emnr.cpp | 390 ++++++++++++++++++++-------------------------- wdsp/emnr.hpp | 115 ++++++++------ 4 files changed, 241 insertions(+), 274 deletions(-) diff --git a/wdsp/calculus.cpp b/wdsp/calculus.cpp index 1a08b51e8..2d6a4c5b6 100644 --- a/wdsp/calculus.cpp +++ b/wdsp/calculus.cpp @@ -29,7 +29,7 @@ warren@wpratt.com namespace WDSP { -const double Calculus::GG[241 * 241] = { +const std::array Calculus::GG = { 7.25654181154076983e-01, 7.05038822098223439e-01, 6.85008217584843870e-01, 6.65545775927326222e-01, 6.46635376294157682e-01, 6.28261355371665386e-01, 6.10408494407843394e-01, 5.93062006626410732e-01, 5.76207525000389742e-01, 5.59831090374464435e-01, 5.43919139925240769e-01, 5.28458495948192608e-01, @@ -14552,7 +14552,7 @@ const double Calculus::GG[241 * 241] = { 1.00000000000000000e+00, 1.00000000000000000e+00, 1.00000000000000000e+00, 1.00000000000000000e+00, 1.00000000000000000e+00 }; -const double Calculus::GGS[241 * 241] = { +const std::array Calculus::GGS = { 8.00014908335353492e-01, 8.00020707540703313e-01, 8.00026700706648830e-01, 8.00032894400760863e-01, 8.00039295417528384e-01, 8.00045910786425396e-01, 8.00052747780268358e-01, 8.00059813923879481e-01, 8.00067117003061101e-01, 8.00074665073896907e-01, 8.00082466472385456e-01, 8.00090529824419749e-01, diff --git a/wdsp/calculus.hpp b/wdsp/calculus.hpp index 93fe9306e..c5d3e3888 100644 --- a/wdsp/calculus.hpp +++ b/wdsp/calculus.hpp @@ -29,6 +29,8 @@ warren@wpratt.com #ifndef wdsp_calculus_h #define wdsp_calculus_h +#include + #include "export.h" namespace WDSP { @@ -36,8 +38,8 @@ namespace WDSP { class WDSP_API Calculus { public: - static const double GG[]; - static const double GGS[]; + static const std::array GG; + static const std::array GGS; }; } // namespace WDSP diff --git a/wdsp/emnr.cpp b/wdsp/emnr.cpp index 35909775b..11ecd1f55 100644 --- a/wdsp/emnr.cpp +++ b/wdsp/emnr.cpp @@ -40,7 +40,7 @@ namespace WDSP { EMNR::AE::AE( int _msize, - double* _lambda_y, + const std::vector& _lambda_y, double _zetaThresh, double _psi ) : @@ -49,20 +49,15 @@ EMNR::AE::AE( zetaThresh(_zetaThresh), psi(_psi) { - nmask = new double[msize]; -} - -EMNR::AE::~AE() -{ - delete[] nmask; + nmask.resize(msize); } EMNR::NPS::NPS( int _incr, double _rate, int _msize, - double* _lambda_y, - double* _lambda_d, + const std::vector& _lambda_y, + std::vector& _lambda_d, double _alpha_pow, double _alpha_Pbar, @@ -78,10 +73,10 @@ EMNR::NPS::NPS( epsH1(_epsH1) { epsH1r = epsH1 / (1.0 + epsH1); - sigma2N = new double[msize]; // (float *)malloc0(nps.msize * sizeof(float)); - PH1y = new double[msize]; // (float *)malloc0(nps.msize * sizeof(float)); - Pbar = new double[msize]; // (float *)malloc0(nps.msize * sizeof(float)); - EN2y = new double[msize]; // (float *)malloc0(nps.msize * sizeof(float)); + sigma2N.resize(msize); // (float *)malloc0(nps.msize * sizeof(float)); + PH1y.resize(msize); // (float *)malloc0(nps.msize * sizeof(float)); + Pbar.resize(msize); // (float *)malloc0(nps.msize * sizeof(float)); + EN2y.resize(msize); // (float *)malloc0(nps.msize * sizeof(float)); for (int i = 0; i < msize; i++) { @@ -90,14 +85,6 @@ EMNR::NPS::NPS( } } -EMNR::NPS::~NPS() -{ - delete[] sigma2N; - delete[] PH1y; - delete[] Pbar; - delete[] EN2y; -} - void EMNR::NPS::LambdaDs() { int k; @@ -114,7 +101,7 @@ void EMNR::NPS::LambdaDs() sigma2N[k] = alpha_pow * sigma2N[k] + (1.0 - alpha_pow) * EN2y[k]; } - std::copy(sigma2N, sigma2N + msize, lambda_d); + std::copy(sigma2N.begin(), sigma2N.end(), lambda_d.begin()); } const std::array EMNR::NP::DVals = { 1.0, 2.0, 5.0, 8.0, 10.0, 15.0, 20.0, 30.0, 40.0, @@ -126,15 +113,15 @@ EMNR::NP::NP( int _incr, double _rate, int _msize, - double* _lambda_y, - double* _lambda_d -) + std::vector& _lambda_y, + std::vector& _lambda_d +) : + incr(_incr), + rate(_rate), + msize(_msize), + lambda_y(_lambda_y), + lambda_d(_lambda_d) { - incr = _incr; - rate = _rate; - msize = _msize; - lambda_y = _lambda_y; - lambda_d = _lambda_d; { double tau = -128.0 / 8000.0 / log(0.7); @@ -194,24 +181,24 @@ EMNR::NP::NP( nsmax[3] = pow(10.0, db / 10.0 * V * incr / rate); } - p = new double[msize]; // (float *)malloc0(np.msize * sizeof(float)); - alphaOptHat = new double[msize]; // (float *)malloc0(np.msize * sizeof(float)); - alphaHat = new double[msize]; // (float *)malloc0(np.msize * sizeof(float)); - sigma2N = new double[msize]; // (float *)malloc0(np.msize * sizeof(float)); - pbar = new double[msize]; // (float *)malloc0(np.msize * sizeof(float)); - p2bar = new double[msize]; // (float *)malloc0(np.msize * sizeof(float)); - Qeq = new double[msize]; // (float *)malloc0(np.msize * sizeof(float)); - bmin = new double[msize]; // (float *)malloc0(np.msize * sizeof(float)); - bmin_sub = new double[msize]; // (float *)malloc0(np.msize * sizeof(float)); - k_mod = new int[msize]; // (int *)malloc0(np.msize * sizeof(int)); - actmin = new double[msize]; // (float *)malloc0(np.msize * sizeof(float)); - actmin_sub = new double[msize]; // (float *)malloc0(np.msize * sizeof(float)); - lmin_flag = new int[msize]; // (int *)malloc0(np.msize * sizeof(int)); - pmin_u = new double[msize]; // (float *)malloc0(np.msize * sizeof(float)); - actminbuff = new double*[U]; // (float**)malloc0(np.U * sizeof(float*)); + p.resize(msize); // (float *)malloc0(np.msize * sizeof(float)); + alphaOptHat.resize(msize); // (float *)malloc0(np.msize * sizeof(float)); + alphaHat.resize(msize); // (float *)malloc0(np.msize * sizeof(float)); + sigma2N.resize(msize); // (float *)malloc0(np.msize * sizeof(float)); + pbar.resize(msize); // (float *)malloc0(np.msize * sizeof(float)); + p2bar.resize(msize); // (float *)malloc0(np.msize * sizeof(float)); + Qeq.resize(msize); // (float *)malloc0(np.msize * sizeof(float)); + bmin.resize(msize); // (float *)malloc0(np.msize * sizeof(float)); + bmin_sub.resize(msize); // (float *)malloc0(np.msize * sizeof(float)); + k_mod.resize(msize); // (int *)malloc0(np.msize * sizeof(int)); + actmin.resize(msize); // (float *)malloc0(np.msize * sizeof(float)); + actmin_sub.resize(msize); // (float *)malloc0(np.msize * sizeof(float)); + lmin_flag.resize(msize); // (int *)malloc0(np.msize * sizeof(int)); + pmin_u.resize(msize); // (float *)malloc0(np.msize * sizeof(float)); + actminbuff.resize(U); // (float**)malloc0(np.U * sizeof(float*)); for (int i = 0; i < U; i++) { - actminbuff[i] = new double[msize]; // (float *)malloc0(np.msize * sizeof(float)); + actminbuff[i].resize(msize); // (float *)malloc0(np.msize * sizeof(float)); } { @@ -224,10 +211,10 @@ EMNR::NP::NP( lambda_y[k] = 0.5; } - std::copy(lambda_y, lambda_y + msize, p); - std::copy(lambda_y, lambda_y + msize, sigma2N); - std::copy(lambda_y, lambda_y + msize, pbar); - std::copy(lambda_y, lambda_y + msize, pmin_u); + std::copy(lambda_y.begin(), lambda_y.end(), p.begin()); + std::copy(lambda_y.begin(), lambda_y.end(), sigma2N.begin()); + std::copy(lambda_y.begin(), lambda_y.end(), pbar.begin()); + std::copy(lambda_y.begin(), lambda_y.end(), pmin_u.begin()); for (k = 0; k < msize; k++) { @@ -240,32 +227,10 @@ EMNR::NP::NP( } } - std::fill(lmin_flag, lmin_flag + msize, 0); + std::fill(lmin_flag.begin(), lmin_flag.end(), 0); } } -EMNR::NP::~NP() -{ - for (int i = 0; i < U; i++) - delete[] (actminbuff[i]); - - delete[] (actminbuff); - delete[] (pmin_u); - delete[] (lmin_flag); - delete[] (actmin_sub); - delete[] (actmin); - delete[] (k_mod); - delete[] (bmin_sub); - delete[] (bmin); - delete[] (Qeq); - delete[] (p2bar); - delete[] (pbar); - delete[] (sigma2N); - delete[] (alphaHat); - delete[] (alphaOptHat); - delete[] (p); -} - void EMNR::NP::interpM ( double* res, double x, @@ -371,7 +336,7 @@ void EMNR::NP::LambdaD() bmin_sub[k] = 1.0 + 2.0 * (V - 1.0) / QeqTildaSub; } - std::fill(k_mod, k_mod + msize, 0); + std::fill(k_mod.begin(), k_mod.end(), 0); for (k = 0; k < msize; k++) { @@ -452,27 +417,27 @@ void EMNR::NP::LambdaD() ++subwc; } - std::copy(sigma2N, sigma2N + msize, lambda_d); + std::copy(sigma2N.begin(), sigma2N.end(), lambda_d.begin()); } EMNR::G::G( int _incr, double _rate, int _msize, - double* _mask, - float* _y -) + std::vector& _mask, + const std::vector& _y +) : + incr(_incr), + rate(_rate), + msize(_msize), + mask(_mask), + y(_y) { - incr = _incr; - rate = _rate; - msize = _msize; - mask = _mask; - y = _y; - lambda_y = new double[msize]; // (float *)malloc0(msize * sizeof(float)); - lambda_d = new double[msize]; // (float *)malloc0(msize * sizeof(float)); - prev_gamma = new double[msize]; // (float *)malloc0(msize * sizeof(float)); - prev_mask = new double[msize]; // (float *)malloc0(msize * sizeof(float)); + lambda_y.resize(msize); // (float *)malloc0(msize * sizeof(float)); + lambda_d.resize(msize); // (float *)malloc0(msize * sizeof(float)); + prev_gamma.resize(msize); // (float *)malloc0(msize * sizeof(float)); + prev_mask.resize(msize); // (float *)malloc0(msize * sizeof(float)); gf1p5 = sqrt(PI) / 2.0; @@ -485,41 +450,113 @@ EMNR::G::G( gamma_max = 1000.0; q = 0.2; - std::fill(prev_mask, prev_mask + msize, 1.0); - std::fill(prev_gamma, prev_gamma + msize, 1.0); + std::fill(prev_mask.begin(), prev_mask.end(), 1.0); + std::fill(prev_gamma.begin(), prev_gamma.end(), 1.0); gmax = 10000.0; - // - GG = new double[241 * 241]; // (float *)malloc0(241 * 241 * sizeof(float)); - GGS = new double[241 * 241]; // (float *)malloc0(241 * 241 * sizeof(float)); + std::copy(Calculus::GG.begin(), Calculus::GG.end(), GG.begin()); + std::copy(Calculus::GGS.begin(), Calculus::GGS.end(), GGS.begin()); + + // We keep this pretty useless part just in case... if ((fileb = fopen("calculus", "rb"))) { - std::size_t lgg = fread(GG, sizeof(float), 241 * 241, fileb); + std::array gg; + std::size_t lgg = fread(&gg[0], sizeof(double), 241 * 241, fileb); if (lgg != 241 * 241) { fprintf(stderr, "GG file has an invalid size\n"); + } else { + std::copy(gg.begin(), gg.end(), GG.begin()); } - std::size_t lggs =fread(GGS, sizeof(float), 241 * 241, fileb); + std::array ggs; + std::size_t lggs =fread(&ggs[0], sizeof(double), 241 * 241, fileb); if (lggs != 241 * 241) { fprintf(stderr, "GGS file has an invalid size\n"); + } else { + std::copy(ggs.begin(), ggs.end(), GGS.begin()); } fclose(fileb); } - else +} + +void EMNR::G::calc_gamma0() +{ + double gamma, eps_hat, v; + + for (int k = 0; k < msize; k++) { - std::copy(Calculus::GG, Calculus::GG + (241 * 241), GG); - std::copy(Calculus::GGS, Calculus::GGS + (241 * 241), GGS); + gamma = std::min (lambda_y[k] / lambda_d[k], gamma_max); + eps_hat = alpha * prev_mask[k] * prev_mask[k] * prev_gamma[k] + + (1.0 - alpha) * std::max (gamma - 1.0f, eps_floor); + v = (eps_hat / (1.0 + eps_hat)) * gamma; + mask[k] = gf1p5 * sqrt (v) / gamma * exp (- 0.5 * v) + * ((1.0 + v) * bessI0 (0.5 * v) + v * bessI1 (0.5 * v)); + { + double v2 = std::min (v, 700.0); + double eta = mask[k] * mask[k] * lambda_y[k] / lambda_d[k]; + double eps = eta / (1.0 - q); + double witchHat = (1.0 - q) / q * exp (v2) / (1.0 + eps); + mask[k] *= witchHat / (1.0 + witchHat); + } + + if (mask[k] > gmax) + mask[k] = gmax; + + if (mask[k] != mask[k]) + mask[k] = 0.01; + + prev_gamma[k] = gamma; + prev_mask[k] = mask[k]; } } -EMNR::G::~G() +void EMNR::G::calc_gamma1() { - delete[] (GGS); - delete[] (GG); - delete[] (prev_mask); - delete[] (prev_gamma); - delete[] (lambda_d); - delete[] (lambda_y); + double gamma, eps_hat, v, ehr; + + for (int k = 0; k < msize; k++) + { + gamma = std::min (lambda_y[k] / lambda_d[k], gamma_max); + eps_hat = alpha * prev_mask[k] * prev_mask[k] * prev_gamma[k] + + (1.0 - alpha) * std::max (gamma - 1.0f, eps_floor); + ehr = eps_hat / (1.0 + eps_hat); + v = ehr * gamma; + + if ((mask[k] = ehr * exp (std::min (700.0, 0.5 * e1xb(v)))) > gmax) + mask[k] = gmax; + + if (mask[k] != mask[k]) + mask[k] = 0.01; + + prev_gamma[k] = gamma; + prev_mask[k] = mask[k]; + } +} + +void EMNR::G::calc_gamma2() +{ + double gamma, eps_hat, eps_p; + + for (int k = 0; k < msize; k++) + { + gamma = std::min(lambda_y[k] / lambda_d[k], gamma_max); + eps_hat = alpha * prev_mask[k] * prev_mask[k] * prev_gamma[k] + + (1.0 - alpha) * std::max(gamma - 1.0f, eps_floor); + eps_p = eps_hat / (1.0 - q); + mask[k] = getKey(GG, gamma, eps_hat) * getKey(GGS, gamma, eps_p); + prev_gamma[k] = gamma; + prev_mask[k] = mask[k]; + } +} + +void EMNR::G::calc_lambda_y() +{ + for (int k = 0; k < msize; k++) + { + double y0 = y[2 * k + 0]; + double y1 = y[2 * k + 1]; + lambda_y[k] = y0 * y0 + y1 * y1; + } } /******************************************************************************************************** @@ -534,7 +571,7 @@ EMNR::G::~G() // Shanjie Zhang and Jianming Jin, "Computation of Special Functions." New York, NY, John Wiley and Sons, // Inc., 1996. [Sample code given in FORTRAN] -double EMNR::bessI0 (double x) +double EMNR::G::bessI0 (double x) { double res, p; @@ -577,7 +614,7 @@ double EMNR::bessI0 (double x) return res; } -double EMNR::bessI1 (double x) +double EMNR::G::bessI1 (double x) { double res, p; @@ -629,7 +666,7 @@ double EMNR::bessI1 (double x) // Shanjie Zhang and Jianming Jin, "Computation of Special Functions." New York, NY, John Wiley and Sons, // Inc., 1996. [Sample code given in FORTRAN] -double EMNR::e1xb (double x) +double EMNR::G::e1xb (double x) { double e1, ga, r, t, t0; int k, m; @@ -702,31 +739,6 @@ void EMNR::calc_window() } } -void EMNR::interpM (double* res, double x, int nvals, double* xvals, double* yvals) -{ - if (x <= xvals[0]) - { - *res = yvals[0]; - } - else if (x >= xvals[nvals - 1]) - { - *res = yvals[nvals - 1]; - } - else - { - int idx = 1; - double xllow, xlhigh, frac; - - while ((x >= xvals[idx]) && (idx < nvals - 1)) - idx++; - - xllow = log10 (xvals[idx - 1]); - xlhigh = log10(xvals[idx]); - frac = (log10 (x) - xllow) / (xlhigh - xllow); - *res = yvals[idx - 1] + frac * (yvals[idx] - yvals[idx - 1]); - } -} - void EMNR::calc() { // float Hvals[18] = { 0.000, 0.150, 0.480, 0.780, 0.980, 1.550, 2.000, 2.300, 2.520, @@ -797,8 +809,8 @@ void EMNR::calc() incr, rate, msize, - mask.data(), - forfftout.data() + mask, + forfftout ); // NP @@ -941,10 +953,10 @@ void EMNR::aepf() ae->nmask[k] /= (float)N; } - std::copy(ae->nmask, ae->nmask + (ae->msize - 2 * n), &mask[n]); + std::copy(ae->nmask.begin(), ae->nmask.end() - 2*n, &mask[n]); } -double EMNR::getKey(double* type, double gamma, double xi) +double EMNR::G::getKey(const std::array& type, double gamma, double xi) { int ngamma1, ngamma2, nxi1 = 0, nxi2 = 0; double tg, tx, dg, dx; @@ -987,22 +999,26 @@ double EMNR::getKey(double* type, double gamma, double xi) dg = (tg - 0.25 * ngamma1) / 0.25; dx = (tx - 0.25 * nxi1) / 0.25; - return (1.0 - dg) * (1.0 - dx) * type[241 * nxi1 + ngamma1] - + (1.0 - dg) * dx * type[241 * nxi2 + ngamma1] - + dg * (1.0 - dx) * type[241 * nxi1 + ngamma2] - + dg * dx * type[241 * nxi2 + ngamma2]; + + std::array ix; + ix[0] = 241 * nxi1 + ngamma1; + ix[1] = 241 * nxi2 + ngamma1; + ix[2] = 241 * nxi1 + ngamma2; + ix[3] = 241 * nxi2 + ngamma2; + + for (auto& ixi : ix) { + ixi = ixi < 0 ? 0 : (ixi >= 241*241 ? 241*241 - 1 : ixi); + } + + return (1.0 - dg) * (1.0 - dx) * type[ix[0]] + + (1.0 - dg) * dx * type[ix[1]] + + dg * (1.0 - dx) * type[ix[2]] + + dg * dx * type[ix[3]]; } void EMNR::calc_gain() { - int k; - - for (k = 0; k < msize; k++) - { - double y0 = g->y[2 * k + 0]; - double y1 = g->y[2 * k + 1]; - g->lambda_y[k] = y0 * y0 + y1 * y1; - } + g->calc_lambda_y(); switch (g->npe_method) { @@ -1017,80 +1033,14 @@ void EMNR::calc_gain() switch (g->gain_method) { case 0: - { - double gamma, eps_hat, v; - - for (k = 0; k < msize; k++) - { - gamma = std::min (g->lambda_y[k] / g->lambda_d[k], g->gamma_max); - eps_hat = g->alpha * g->prev_mask[k] * g->prev_mask[k] * g->prev_gamma[k] - + (1.0 - g->alpha) * std::max (gamma - 1.0f, g->eps_floor); - v = (eps_hat / (1.0 + eps_hat)) * gamma; - g->mask[k] = g->gf1p5 * sqrt (v) / gamma * exp (- 0.5 * v) - * ((1.0 + v) * bessI0 (0.5 * v) + v * bessI1 (0.5 * v)); - { - double v2 = std::min (v, 700.0); - double eta = g->mask[k] * g->mask[k] * g->lambda_y[k] / g->lambda_d[k]; - double eps = eta / (1.0 - g->q); - double witchHat = (1.0 - g->q) / g->q * exp (v2) / (1.0 + eps); - g->mask[k] *= witchHat / (1.0 + witchHat); - } - - if (g->mask[k] > g->gmax) - g->mask[k] = g->gmax; - - if (g->mask[k] != g->mask[k]) - g->mask[k] = 0.01; - - g->prev_gamma[k] = gamma; - g->prev_mask[k] = g->mask[k]; - } - - break; - } - + g->calc_gamma0(); + break; case 1: - { - double gamma, eps_hat, v, ehr; - - for (k = 0; k < msize; k++) - { - gamma = std::min (g->lambda_y[k] / g->lambda_d[k], g->gamma_max); - eps_hat = g->alpha * g->prev_mask[k] * g->prev_mask[k] * g->prev_gamma[k] - + (1.0 - g->alpha) * std::max (gamma - 1.0f, g->eps_floor); - ehr = eps_hat / (1.0 + eps_hat); - v = ehr * gamma; - - if ((g->mask[k] = ehr * exp (std::min (700.0, 0.5 * e1xb(v)))) > g->gmax) - g->mask[k] = g->gmax; - - if (g->mask[k] != g->mask[k]) - g->mask[k] = 0.01; - - g->prev_gamma[k] = gamma; - g->prev_mask[k] = g->mask[k]; - } - - break; - } - + g->calc_gamma1(); + break; case 2: - { - double gamma, eps_hat, eps_p; - - for (k = 0; k < msize; k++) - { - gamma = std::min(g->lambda_y[k] / g->lambda_d[k], g->gamma_max); - eps_hat = g->alpha * g->prev_mask[k] * g->prev_mask[k] * g->prev_gamma[k] - + (1.0 - g->alpha) * std::max(gamma - 1.0f, g->eps_floor); - eps_p = eps_hat / (1.0 - g->q); - g->mask[k] = getKey(g->GG, gamma, eps_hat) * getKey(g->GGS, gamma, eps_p); - g->prev_gamma[k] = gamma; - g->prev_mask[k] = g->mask[k]; - } - - break; - } + g->calc_gamma2(); + break; } if (g->ae_run) diff --git a/wdsp/emnr.hpp b/wdsp/emnr.hpp index fc4a6efb9..58f99390e 100644 --- a/wdsp/emnr.hpp +++ b/wdsp/emnr.hpp @@ -74,20 +74,21 @@ public: int saveidx; fftwf_plan Rfor; fftwf_plan Rrev; + struct G { int incr; double rate; int msize; - double* mask; - float* y; + std::vector& mask; + const std::vector& y; int gain_method; int npe_method; int ae_run; - double* lambda_y; - double* lambda_d; - double* prev_mask; - double* prev_gamma; + std::vector lambda_y; + std::vector lambda_d; + std::vector prev_mask; + std::vector prev_gamma; double gf1p5; double alpha; double eps_floor; @@ -95,27 +96,40 @@ public: double q; double gmax; // - double* GG; - double* GGS; + std::array GG; + std::array GGS; FILE* fileb; + G( int incr, double rate, int msize, - double* mask, - float* y + std::vector& mask, + const std::vector& y ); G(const G&) = delete; G& operator=(const G& other) = delete; - ~G(); + ~G() = default; + + void calc_gamma0(); + void calc_gamma1(); + void calc_gamma2(); + void calc_lambda_y(); + + private: + static double getKey(const std::array& type, double gamma, double xi); + static double e1xb (double x); + static double bessI0 (double x); + static double bessI1 (double x); } *g; + struct NP { int incr; double rate; int msize; - double* lambda_y; - double* lambda_d; + std::vector& lambda_y; + std::vector& lambda_d; double alphaCsmooth; double alphaMax; double alphaCmin; @@ -128,38 +142,40 @@ public: int U; int V; int D; - double* p; - double* alphaOptHat; + std::vector p; + std::vector alphaOptHat; double alphaC; - double* alphaHat; - double* sigma2N; - double* pbar; - double* p2bar; - double* Qeq; + std::vector alphaHat; + std::vector sigma2N; + std::vector pbar; + std::vector p2bar; + std::vector Qeq; double MofD; double MofV; std::array invQbar_points; std::array nsmax; - double* bmin; - double* bmin_sub; - int* k_mod; - double* actmin; - double* actmin_sub; + std::vector bmin; + std::vector bmin_sub; + std::vector k_mod; + std::vector actmin; + std::vector actmin_sub; int subwc; - int* lmin_flag; - double* pmin_u; - double** actminbuff; + std::vector lmin_flag; + std::vector pmin_u; + std::vector> actminbuff; int amb_idx; + NP( int incr, double rate, int msize, - double* lambda_y, - double* lambda_d + std::vector& lambda_y, + std::vector& lambda_d ); NP(const NP&) = delete; NP& operator=(const NP& other) = delete; - ~NP(); + ~NP() = default; + void LambdaD(); private: @@ -173,55 +189,59 @@ public: const std::array& yvals ); } *np; + struct NPS { int incr; double rate; int msize; - double* lambda_y; - double* lambda_d; + const std::vector& lambda_y; + std::vector& lambda_d; double alpha_pow; double alpha_Pbar; double epsH1; double epsH1r; - double* sigma2N; - double* PH1y; - double* Pbar; - double* EN2y; + std::vector sigma2N; + std::vector PH1y; + std::vector Pbar; + std::vector EN2y; + NPS( int incr, double rate, int msize, - double* lambda_y, - double* lambda_d, - + const std::vector& lambda_y, + std::vector& lambda_d, double alpha_pow, double alpha_Pbar, double epsH1 ); NPS(const NPS&) = delete; NPS& operator=(const NPS& other) = delete; - ~NPS(); + ~NPS() = default; + void LambdaDs(); } *nps; + struct AE { int msize; - double* lambda_y; + const std::vector& lambda_y; double zetaThresh; double psi; - double* nmask; + std::vector nmask; + AE( int msize, - double* lambda_y, + const std::vector& lambda_y, double zetaThresh, double psi ); AE(const AE&) = delete; AE& operator=(const AE& other) = delete; - ~AE(); + ~AE() = default; } *ae; EMNR( @@ -256,11 +276,6 @@ public: void setAePsi(double psi); private: - static double bessI0 (double x); - static double bessI1 (double x); - static double e1xb (double x); - static void interpM (double* res, double x, int nvals, double* xvals, double* yvals); - static double getKey(double* type, double gamma, double xi); void calc_window(); void calc(); void decalc();