1
0
mirror of https://github.com/f4exb/sdrangel.git synced 2024-12-23 01:55:48 -05:00

EMNR: make internal calculations in double precision

This commit is contained in:
f4exb 2024-07-06 22:40:25 +02:00
parent 043d9da47e
commit ed3a57c92c
4 changed files with 167 additions and 138 deletions

View File

@ -29,7 +29,7 @@ warren@wpratt.com
namespace WDSP {
const float Calculus::GG[241 * 241] = {
const double Calculus::GG[241 * 241] = {
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 float Calculus::GG[241 * 241] = {
1.00000000000000000e+00, 1.00000000000000000e+00, 1.00000000000000000e+00, 1.00000000000000000e+00,
1.00000000000000000e+00 };
const float Calculus::GGS[241 * 241] = {
const double Calculus::GGS[241 * 241] = {
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,

View File

@ -36,8 +36,8 @@ namespace WDSP {
class WDSP_API Calculus
{
public:
static const float GG[];
static const float GGS[];
static const double GG[];
static const double GGS[];
};
} // namespace WDSP

View File

@ -48,9 +48,9 @@ namespace WDSP {
// Shanjie Zhang and Jianming Jin, "Computation of Special Functions." New York, NY, John Wiley and Sons,
// Inc., 1996. [Sample code given in FORTRAN]
float EMNR::bessI0 (float x)
double EMNR::bessI0 (double x)
{
float res, p;
double res, p;
if (x == 0.0)
res = 1.0;
else
@ -86,10 +86,10 @@ float EMNR::bessI0 (float x)
return res;
}
float EMNR::bessI1 (float x)
double EMNR::bessI1 (double x)
{
float res, p;
double res, p;
if (x == 0.0)
res = 0.0;
else
@ -132,9 +132,9 @@ float EMNR::bessI1 (float x)
// Shanjie Zhang and Jianming Jin, "Computation of Special Functions." New York, NY, John Wiley and Sons,
// Inc., 1996. [Sample code given in FORTRAN]
float EMNR::e1xb (float x)
double EMNR::e1xb (double x)
{
float e1, ga, r, t, t0;
double e1, ga, r, t, t0;
int k, m;
if (x == 0.0)
e1 = 1.0e300;
@ -193,7 +193,7 @@ void EMNR::calc_window (EMNR *a)
}
}
void EMNR::interpM (float* res, float x, int nvals, float* xvals, float* yvals)
void EMNR::interpM (double* res, double x, int nvals, double* xvals, double* yvals)
{
if (x <= xvals[0])
*res = yvals[0];
@ -202,7 +202,7 @@ void EMNR::interpM (float* res, float x, int nvals, float* xvals, float* yvals)
else
{
int idx = 0;
float xllow, xlhigh, frac;
double xllow, xlhigh, frac;
while (x >= xvals[idx]) idx++;
xllow = log10 (xvals[idx - 1]);
xlhigh = log10(xvals[idx]);
@ -214,9 +214,9 @@ void EMNR::interpM (float* res, float x, int nvals, float* xvals, float* yvals)
void EMNR::calc_emnr(EMNR *a)
{
int i;
float Dvals[18] = { 1.0, 2.0, 5.0, 8.0, 10.0, 15.0, 20.0, 30.0, 40.0,
double Dvals[18] = { 1.0, 2.0, 5.0, 8.0, 10.0, 15.0, 20.0, 30.0, 40.0,
60.0, 80.0, 120.0, 140.0, 160.0, 180.0, 220.0, 260.0, 300.0 };
float Mvals[18] = { 0.000, 0.260, 0.480, 0.580, 0.610, 0.668, 0.705, 0.762, 0.800,
double Mvals[18] = { 0.000, 0.260, 0.480, 0.580, 0.610, 0.668, 0.705, 0.762, 0.800,
0.841, 0.865, 0.890, 0.900, 0.910, 0.920, 0.930, 0.935, 0.940 };
// float Hvals[18] = { 0.000, 0.150, 0.480, 0.780, 0.980, 1.550, 2.000, 2.300, 2.520,
// 3.100, 3.380, 4.150, 4.350, 4.250, 3.900, 4.100, 4.700, 5.000 };
@ -246,7 +246,8 @@ void EMNR::calc_emnr(EMNR *a)
a->inaccum = new float[a->iasize]; // (float *)malloc0(a->iasize * sizeof(float));
a->forfftin = new float[a->fsize]; // (float *)malloc0(a->fsize * sizeof(float));
a->forfftout = new float[a->msize * 2]; // (float *)malloc0(a->msize * sizeof(complex));
a->mask = new float[a->msize]; // (float *)malloc0(a->msize * sizeof(float));
a->mask = new double[a->msize]; // (float *)malloc0(a->msize * sizeof(float));
std::fill(a->mask, a->mask + a->msize, 1.0);
a->revfftin = new float[a->msize * 2]; // (float *)malloc0(a->msize * sizeof(complex));
a->revfftout = new float[a->fsize]; // (float *)malloc0(a->fsize * sizeof(float));
a->save = new float*[a->ovrlp]; // (float **)malloc0(a->ovrlp * sizeof(float *));
@ -262,10 +263,10 @@ void EMNR::calc_emnr(EMNR *a)
a->g.msize = a->msize;
a->g.mask = a->mask;
a->g.y = a->forfftout;
a->g.lambda_y = new float[a->msize]; // (float *)malloc0(a->msize * sizeof(float));
a->g.lambda_d = new float[a->msize]; // (float *)malloc0(a->msize * sizeof(float));
a->g.prev_gamma = new float[a->msize]; // (float *)malloc0(a->msize * sizeof(float));
a->g.prev_mask = new float[a->msize]; // (float *)malloc0(a->msize * sizeof(float));
a->g.lambda_y = new double[a->msize]; // (float *)malloc0(a->msize * sizeof(float));
a->g.lambda_d = new double[a->msize]; // (float *)malloc0(a->msize * sizeof(float));
a->g.prev_gamma = new double[a->msize]; // (float *)malloc0(a->msize * sizeof(float));
a->g.prev_mask = new double[a->msize]; // (float *)malloc0(a->msize * sizeof(float));
a->g.gf1p5 = sqrt(PI) / 2.0;
{
@ -282,8 +283,8 @@ void EMNR::calc_emnr(EMNR *a)
}
a->g.gmax = 10000.0;
//
a->g.GG = new float[241 * 241]; // (float *)malloc0(241 * 241 * sizeof(float));
a->g.GGS = new float[241 * 241]; // (float *)malloc0(241 * 241 * sizeof(float));
a->g.GG = new double[241 * 241]; // (float *)malloc0(241 * 241 * sizeof(float));
a->g.GGS = new double[241 * 241]; // (float *)malloc0(241 * 241 * sizeof(float));
if ((a->g.fileb = fopen("calculus", "rb")))
{
std::size_t lgg = fread(a->g.GG, sizeof(float), 241 * 241, a->g.fileb);
@ -356,23 +357,23 @@ void EMNR::calc_emnr(EMNR *a)
a->np.nsmax[3] = pow(10.0, db / 10.0 * a->np.V * a->np.incr / a->np.rate);
}
a->np.p = new float[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.alphaOptHat = new float[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.alphaHat = new float[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.sigma2N = new float[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.pbar = new float[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.p2bar = new float[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.Qeq = new float[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.bmin = new float[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.bmin_sub = new float[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.p = new double[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.alphaOptHat = new double[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.alphaHat = new double[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.sigma2N = new double[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.pbar = new double[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.p2bar = new double[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.Qeq = new double[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.bmin = new double[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.bmin_sub = new double[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.k_mod = new int[a->np.msize]; // (int *)malloc0(a->np.msize * sizeof(int));
a->np.actmin = new float[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.actmin_sub = new float[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.actmin = new double[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.actmin_sub = new double[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.lmin_flag = new int[a->np.msize]; // (int *)malloc0(a->np.msize * sizeof(int));
a->np.pmin_u = new float[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.actminbuff = new float*[a->np.U]; // (float**)malloc0(a->np.U * sizeof(float*));
a->np.pmin_u = new double[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.actminbuff = new double*[a->np.U]; // (float**)malloc0(a->np.U * sizeof(float*));
for (i = 0; i < a->np.U; i++)
a->np.actminbuff[i] = new float[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
a->np.actminbuff[i] = new double[a->np.msize]; // (float *)malloc0(a->np.msize * sizeof(float));
{
int k, ku;
@ -412,10 +413,10 @@ void EMNR::calc_emnr(EMNR *a)
a->nps.epsH1 = pow(10.0, 15.0 / 10.0);
a->nps.epsH1r = a->nps.epsH1 / (1.0 + a->nps.epsH1);
a->nps.sigma2N = new float[a->nps.msize]; // (float *)malloc0(a->nps.msize * sizeof(float));
a->nps.PH1y = new float[a->nps.msize]; // (float *)malloc0(a->nps.msize * sizeof(float));
a->nps.Pbar = new float[a->nps.msize]; // (float *)malloc0(a->nps.msize * sizeof(float));
a->nps.EN2y = new float[a->nps.msize]; // (float *)malloc0(a->nps.msize * sizeof(float));
a->nps.sigma2N = new double[a->nps.msize]; // (float *)malloc0(a->nps.msize * sizeof(float));
a->nps.PH1y = new double[a->nps.msize]; // (float *)malloc0(a->nps.msize * sizeof(float));
a->nps.Pbar = new double[a->nps.msize]; // (float *)malloc0(a->nps.msize * sizeof(float));
a->nps.EN2y = new double[a->nps.msize]; // (float *)malloc0(a->nps.msize * sizeof(float));
for (i = 0; i < a->nps.msize; i++)
{
@ -429,7 +430,7 @@ void EMNR::calc_emnr(EMNR *a)
a->ae.zetaThresh = 0.75;
a->ae.psi = 10.0;
a->ae.nmask = new float[a->ae.msize]; // (float *)malloc0(a->ae.msize * sizeof(float));
a->ae.nmask = new double[a->ae.msize]; // (float *)malloc0(a->ae.msize * sizeof(float));
}
void EMNR::decalc_emnr(EMNR *a)
@ -482,8 +483,21 @@ void EMNR::decalc_emnr(EMNR *a)
delete[] (a->window);
}
EMNR* EMNR::create_emnr (int run, int position, int size, float* in, float* out, int fsize, int ovrlp,
int rate, int wintype, float gain, int gain_method, int npe_method, int ae_run)
EMNR* EMNR::create_emnr (
int run,
int position,
int size,
float* in,
float* out,
int fsize,
int ovrlp,
int rate,
int wintype,
float gain,
int gain_method,
int npe_method,
int ae_run
)
{
EMNR *a = new EMNR;
@ -528,17 +542,17 @@ void EMNR::destroy_emnr (EMNR *a)
void EMNR::LambdaD(EMNR *a)
{
int k;
float f0, f1, f2, f3;
float sum_prev_p;
float sum_lambda_y;
float alphaCtilda;
float sum_prev_sigma2N;
float alphaMin, SNR;
float beta, varHat, invQeq;
float invQbar;
float bc;
float QeqTilda, QeqTildaSub;
float noise_slope_max;
double f0, f1, f2, f3;
double sum_prev_p;
double sum_lambda_y;
double alphaCtilda;
double sum_prev_sigma2N;
double alphaMin, SNR;
double beta, varHat, invQeq;
double invQbar;
double bc;
double QeqTilda, QeqTildaSub;
double noise_slope_max;
sum_prev_p = 0.0;
sum_lambda_y = 0.0;
@ -658,7 +672,7 @@ void EMNR::LambdaDs (EMNR *a)
a->nps.PH1y[k] = 1.0 / (1.0 + (1.0 + a->nps.epsH1) * exp (- a->nps.epsH1r * a->nps.lambda_y[k] / a->nps.sigma2N[k]));
a->nps.Pbar[k] = a->nps.alpha_Pbar * a->nps.Pbar[k] + (1.0 - a->nps.alpha_Pbar) * a->nps.PH1y[k];
if (a->nps.Pbar[k] > 0.99)
a->nps.PH1y[k] = std::min (a->nps.PH1y[k], 0.99f);
a->nps.PH1y[k] = std::min (a->nps.PH1y[k], 0.99);
a->nps.EN2y[k] = (1.0 - a->nps.PH1y[k]) * a->nps.lambda_y[k] + a->nps.PH1y[k] * a->nps.sigma2N[k];
a->nps.sigma2N[k] = a->nps.alpha_pow * a->nps.sigma2N[k] + (1.0 - a->nps.alpha_pow) * a->nps.EN2y[k];
}
@ -697,12 +711,12 @@ void EMNR::aepf(EMNR *a)
std::copy(a->ae.nmask, a->ae.nmask + (a->ae.msize - 2 * n), a->mask + n);
}
float EMNR::getKey(float* type, float gamma, float xi)
double EMNR::getKey(double* type, double gamma, double xi)
{
int ngamma1, ngamma2, nxi1, nxi2;
float tg, tx, dg, dx;
const float dmin = 0.001;
const float dmax = 1000.0;
int ngamma1, ngamma2, nxi1 = 0, nxi2 = 0;
double tg, tx, dg, dx;
const double dmin = 0.001;
const double dmax = 1000.0;
if (gamma <= dmin)
{
ngamma1 = ngamma2 = 0;
@ -748,7 +762,9 @@ void EMNR::calc_gain (EMNR *a)
int k;
for (k = 0; k < a->g.msize; k++)
{
a->g.lambda_y[k] = a->g.y[2 * k + 0] * a->g.y[2 * k + 0] + a->g.y[2 * k + 1] * a->g.y[2 * k + 1];
double y0 = a->g.y[2 * k + 0];
double y1 = a->g.y[2 * k + 1];
a->g.lambda_y[k] = y0 * y0 + y1 * y1;
}
switch (a->g.npe_method)
{
@ -763,7 +779,7 @@ void EMNR::calc_gain (EMNR *a)
{
case 0:
{
float gamma, eps_hat, v;
double gamma, eps_hat, v;
for (k = 0; k < a->msize; k++)
{
gamma = std::min (a->g.lambda_y[k] / a->g.lambda_d[k], a->g.gamma_max);
@ -773,10 +789,10 @@ void EMNR::calc_gain (EMNR *a)
a->g.mask[k] = a->g.gf1p5 * sqrt (v) / gamma * exp (- 0.5 * v)
* ((1.0 + v) * bessI0 (0.5 * v) + v * bessI1 (0.5 * v));
{
float v2 = std::min (v, 700.0f);
float eta = a->g.mask[k] * a->g.mask[k] * a->g.lambda_y[k] / a->g.lambda_d[k];
float eps = eta / (1.0 - a->g.q);
float witchHat = (1.0 - a->g.q) / a->g.q * exp (v2) / (1.0 + eps);
double v2 = std::min (v, 700.0);
double eta = a->g.mask[k] * a->g.mask[k] * a->g.lambda_y[k] / a->g.lambda_d[k];
double eps = eta / (1.0 - a->g.q);
double witchHat = (1.0 - a->g.q) / a->g.q * exp (v2) / (1.0 + eps);
a->g.mask[k] *= witchHat / (1.0 + witchHat);
}
if (a->g.mask[k] > a->g.gmax) a->g.mask[k] = a->g.gmax;
@ -788,7 +804,7 @@ void EMNR::calc_gain (EMNR *a)
}
case 1:
{
float gamma, eps_hat, v, ehr;
double gamma, eps_hat, v, ehr;
for (k = 0; k < a->g.msize; k++)
{
gamma = std::min (a->g.lambda_y[k] / a->g.lambda_d[k], a->g.gamma_max);
@ -805,7 +821,7 @@ void EMNR::calc_gain (EMNR *a)
}
case 2:
{
float gamma, eps_hat, eps_p;
double gamma, eps_hat, eps_p;
for (k = 0; k < a->msize; k++)
{
gamma = std::min(a->g.lambda_y[k] / a->g.lambda_d[k], a->g.gamma_max);
@ -952,14 +968,14 @@ void EMNR::SetEMNRPosition (RXA& rxa, int position)
rxa.csDSP.unlock();
}
void EMNR::SetEMNRaeZetaThresh (RXA& rxa, float zetathresh)
void EMNR::SetEMNRaeZetaThresh (RXA& rxa, double zetathresh)
{
rxa.csDSP.lock();
rxa.emnr.p->ae.zetaThresh = zetathresh;
rxa.csDSP.unlock();
}
void EMNR::SetEMNRaePsi (RXA& rxa, float psi)
void EMNR::SetEMNRaePsi (RXA& rxa, double psi)
{
rxa.csDSP.lock();
rxa.emnr.p->ae.psi = psi;

View File

@ -52,7 +52,7 @@ public:
float* forfftin;
float* forfftout;
int msize;
float* mask;
double* mask;
float* revfftin;
float* revfftout;
float** save;
@ -76,95 +76,108 @@ public:
int gain_method;
int npe_method;
int ae_run;
float msize;
float* mask;
double msize;
double* mask;
float* y;
float* lambda_y;
float* lambda_d;
float* prev_mask;
float* prev_gamma;
float gf1p5;
float alpha;
float eps_floor;
float gamma_max;
float q;
float gmax;
double* lambda_y;
double* lambda_d;
double* prev_mask;
double* prev_gamma;
double gf1p5;
double alpha;
double eps_floor;
double gamma_max;
double q;
double gmax;
//
float* GG;
float* GGS;
double* GG;
double* GGS;
FILE* fileb;
} g;
struct _npest
{
int incr;
float rate;
double rate;
int msize;
float* lambda_y;
float* lambda_d;
float* p;
float* alphaOptHat;
float alphaC;
float alphaCsmooth;
float alphaCmin;
float* alphaHat;
float alphaMax;
float* sigma2N;
float alphaMin_max_value;
float snrq;
float betamax;
float* pbar;
float* p2bar;
float invQeqMax;
float av;
float* Qeq;
double* lambda_y;
double* lambda_d;
double* p;
double* alphaOptHat;
double alphaC;
double alphaCsmooth;
double alphaCmin;
double* alphaHat;
double alphaMax;
double* sigma2N;
double alphaMin_max_value;
double snrq;
double betamax;
double* pbar;
double* p2bar;
double invQeqMax;
double av;
double* Qeq;
int U;
float Dtime;
double Dtime;
int V;
int D;
float MofD;
float MofV;
float* bmin;
float* bmin_sub;
double MofD;
double MofV;
double* bmin;
double* bmin_sub;
int* k_mod;
float* actmin;
float* actmin_sub;
double* actmin;
double* actmin_sub;
int subwc;
int* lmin_flag;
float* pmin_u;
float invQbar_points[4];
float nsmax[4];
float** actminbuff;
double* pmin_u;
double invQbar_points[4];
double nsmax[4];
double** actminbuff;
int amb_idx;
} np;
struct _npests
{
int incr;
float rate;
double rate;
int msize;
float* lambda_y;
float* lambda_d;
double* lambda_y;
double* lambda_d;
float alpha_pow;
float alpha_Pbar;
float epsH1;
float epsH1r;
double alpha_pow;
double alpha_Pbar;
double epsH1;
double epsH1r;
float* sigma2N;
float* PH1y;
float* Pbar;
float* EN2y;
double* sigma2N;
double* PH1y;
double* Pbar;
double* EN2y;
} nps;
struct _ae
{
int msize;
float* lambda_y;
float zetaThresh;
float psi;
float* nmask;
double* lambda_y;
double zetaThresh;
double psi;
double* nmask;
} ae;
static EMNR* create_emnr (int run, int position, int size, float* in, float* out, int fsize, int ovrlp,
int rate, int wintype, float gain, int gain_method, int npe_method, int ae_run);
static EMNR* create_emnr (
int run,
int position,
int size,
float* in,
float* out,
int fsize,
int ovrlp,
int rate,
int wintype,
float gain,
int gain_method,
int npe_method,
int ae_run
);
static void destroy_emnr (EMNR *a);
static void flush_emnr (EMNR *a);
static void xemnr (EMNR *a, int pos);
@ -177,21 +190,21 @@ public:
static void SetEMNRnpeMethod (RXA& rxa, int method);
static void SetEMNRaeRun (RXA& rxa, int run);
static void SetEMNRPosition (RXA& rxa, int position);
static void SetEMNRaeZetaThresh (RXA& rxa, float zetathresh);
static void SetEMNRaePsi (RXA& rxa, float psi);
static void SetEMNRaeZetaThresh (RXA& rxa, double zetathresh);
static void SetEMNRaePsi (RXA& rxa, double psi);
private:
static float bessI0 (float x);
static float bessI1 (float x);
static float e1xb (float x);
static double bessI0 (double x);
static double bessI1 (double x);
static double e1xb (double x);
static void calc_window (EMNR *a);
static void interpM (float* res, float x, int nvals, float* xvals, float* yvals);
static void interpM (double* res, double x, int nvals, double* xvals, double* yvals);
static void calc_emnr(EMNR *a);
static void decalc_emnr(EMNR *a);
static void LambdaD(EMNR *a);
static void LambdaDs (EMNR *a);
static void aepf(EMNR *a);
static float getKey(float* type, float gamma, float xi);
static double getKey(double* type, double gamma, double xi);
static void calc_gain (EMNR *a);
};