VE7UWU/sdrangel
1
0
Fork 0
mirror of https://github.com/f4exb/sdrangel.git synced 2024-12-23 01:55:48 -05:00
Code Issues Projects Releases Wiki Activity

WDSP: FIR: calculate on double precision as much as possible

Browse Source
This commit is contained in:
f4exb 2024-07-17 03:04:47 +02:00
parent 6d3516482e
commit 3f800dd0a9
4 changed files with 40 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
wdsp/emnr.cpp
View File

@ -24,6 +24,8 @@ The author can be reached by email at
warren@wpratt.com warren@wpratt.com
*/ */
#include <limits>
#include "comm.hpp" #include "comm.hpp"
#include "calculus.hpp" #include "calculus.hpp"
#include "emnr.hpp" #include "emnr.hpp"
@ -273,7 +275,7 @@ void EMNR::calc_emnr(EMNR *a)
float tau = -128.0 / 8000.0 / log(0.98); float tau = -128.0 / 8000.0 / log(0.98);
a->g.alpha = exp(-a->incr / a->rate / tau); a->g.alpha = exp(-a->incr / a->rate / tau);
} }
a->g.eps_floor = 1.0e-300; a->g.eps_floor = std::numeric_limits<double>::min();
a->g.gamma_max = 1000.0; a->g.gamma_max = 1000.0;
a->g.q = 0.2; a->g.q = 0.2;
for (i = 0; i < a->g.msize; i++) for (i = 0; i < a->g.msize; i++)

51
wdsp/fir.cpp
View File

@ -25,6 +25,9 @@ warren@pratt.one
*/ */
#define _CRT_SECURE_NO_WARNINGS #define _CRT_SECURE_NO_WARNINGS
#include <limits>
#include "fftw3.h" #include "fftw3.h"
#include "comm.hpp" #include "comm.hpp"
#include "fir.hpp" #include "fir.hpp"
@ -54,15 +57,15 @@ float* FIR::fftcv_mults (int NM, float* c_impulse)
float* FIR::get_fsamp_window(int N, int wintype) float* FIR::get_fsamp_window(int N, int wintype)
{ {
int i; int i;
float arg0, arg1; double arg0, arg1;
float* window = new float[N]; // (float *) malloc0 (N * sizeof(float)); float* window = new float[N]; // (float *) malloc0 (N * sizeof(float));
switch (wintype) switch (wintype)
{ {
case 0: case 0:
arg0 = 2.0 * PI / ((float)N - 1.0); arg0 = 2.0 * PI / ((double)N - 1.0);
for (i = 0; i < N; i++) for (i = 0; i < N; i++)
{ {
arg1 = cos(arg0 * (float)i); arg1 = cos(arg0 * (double)i);
window[i] = +0.21747 window[i] = +0.21747
+ arg1 * (-0.45325 + arg1 * (-0.45325
+ arg1 * (+0.28256 + arg1 * (+0.28256
@ -70,10 +73,10 @@ float* FIR::get_fsamp_window(int N, int wintype)
} }
break; break;
case 1: case 1:
arg0 = 2.0 * PI / ((float)N - 1.0); arg0 = 2.0 * PI / ((double)N - 1.0);
for (i = 0; i < N; ++i) for (i = 0; i < N; ++i)
{ {
arg1 = cos(arg0 * (float)i); arg1 = cos(arg0 * (double)i);
window[i] = +6.3964424114390378e-02 window[i] = +6.3964424114390378e-02
+ arg1 * (-2.3993864599352804e-01 + arg1 * (-2.3993864599352804e-01
+ arg1 * (+3.5015956323820469e-01 + arg1 * (+3.5015956323820469e-01
@ -90,11 +93,11 @@ float* FIR::get_fsamp_window(int N, int wintype)
return window; return window;
} }
float* FIR::fir_fsamp_odd (int N, float* A, int rtype, float scale, int wintype) float* FIR::fir_fsamp_odd (int N, float* A, int rtype, double scale, int wintype)
{ {
int i, j; int i, j;
int mid = (N - 1) / 2; int mid = (N - 1) / 2;
float mag, phs; double mag, phs;
float* window; float* window;
float *fcoef = new float[N * 2]; float *fcoef = new float[N * 2];
float *c_impulse = new float[N * 2]; float *c_impulse = new float[N * 2];
@ -105,11 +108,11 @@ float* FIR::fir_fsamp_odd (int N, float* A, int rtype, float scale, int wintype)
FFTW_BACKWARD, FFTW_BACKWARD,
FFTW_PATIENT FFTW_PATIENT
); );
float local_scale = 1.0 / (float)N; double local_scale = 1.0 / (double) N;
for (i = 0; i <= mid; i++) for (i = 0; i <= mid; i++)
{ {
mag = A[i] * local_scale; mag = A[i] * local_scale;
phs = - (float)mid * TWOPI * (float)i / (float)N; phs = - (double)mid * TWOPI * (double)i / (double)N;
fcoef[2 * i + 0] = mag * cos (phs); fcoef[2 * i + 0] = mag * cos (phs);
fcoef[2 * i + 1] = mag * sin (phs); fcoef[2 * i + 1] = mag * sin (phs);
} }
@ -140,10 +143,10 @@ float* FIR::fir_fsamp_odd (int N, float* A, int rtype, float scale, int wintype)
return c_impulse; return c_impulse;
} }
float* FIR::fir_fsamp (int N, float* A, int rtype, float scale, int wintype) float* FIR::fir_fsamp (int N, float* A, int rtype, double scale, int wintype)
{ {
int n, i, j, k; int n, i, j, k;
float sum; double sum;
float* window; float* window;
float *c_impulse = new float[N * 2]; // (float *) malloc0 (N * sizeof (complex)); float *c_impulse = new float[N * 2]; // (float *) malloc0 (N * sizeof (complex));
@ -166,7 +169,7 @@ float* FIR::fir_fsamp (int N, float* A, int rtype, float scale, int wintype)
} }
else else
{ {
float M = (float)(N - 1) / 2.0; double M = (double)(N - 1) / 2.0;
for (n = 0; n < N / 2; n++) for (n = 0; n < N / 2; n++)
{ {
sum = 0.0; sum = 0.0;
@ -200,18 +203,18 @@ float* FIR::fir_fsamp (int N, float* A, int rtype, float scale, int wintype)
return c_impulse; return c_impulse;
} }
float* FIR::fir_bandpass (int N, float f_low, float f_high, float samplerate, int wintype, int rtype, float scale) float* FIR::fir_bandpass (int N, double f_low, double f_high, double samplerate, int wintype, int rtype, double scale)
{ {
float *c_impulse = new float[N * 2]; // (float *) malloc0 (N * sizeof (complex)); float *c_impulse = new float[N * 2]; // (float *) malloc0 (N * sizeof (complex));
float ft = (f_high - f_low) / (2.0 * samplerate); double ft = (f_high - f_low) / (2.0 * samplerate);
float ft_rad = TWOPI * ft; double ft_rad = TWOPI * ft;
float w_osc = PI * (f_high + f_low) / samplerate; double w_osc = PI * (f_high + f_low) / samplerate;
int i, j; int i, j;
float m = 0.5 * (float)(N - 1); double m = 0.5 * (double)(N - 1);
float delta = PI / m; double delta = PI / m;
float cosphi; double cosphi;
float posi, posj; double posi, posj;
float sinc, window, coef; double sinc, window, coef;
if (N & 1) if (N & 1)
{ {
@ -228,8 +231,8 @@ float* FIR::fir_bandpass (int N, float f_low, float f_high, float samplerate, in
} }
for (i = (N + 1) / 2, j = N / 2 - 1; i < N; i++, j--) for (i = (N + 1) / 2, j = N / 2 - 1; i < N; i++, j--)
{ {
posi = (float)i - m; posi = (double)i - m;
posj = (float)j - m; posj = (double)j - m;
sinc = sin (ft_rad * posi) / (PI * posi); sinc = sin (ft_rad * posi) / (PI * posi);
switch (wintype) switch (wintype)
{ {
@ -384,7 +387,7 @@ void FIR::mp_imp (int N, float* fir, float* mpfir, int pfactor, int polarity)
if (mag[i] > 0.0) if (mag[i] > 0.0)
ana[2 * i + 0] = log (mag[i]); ana[2 * i + 0] = log (mag[i]);
else else
ana[2 * i + 0] = log (1.0e-300); ana[2 * i + 0] = log (std::numeric_limits<float>::min());
} }
analytic (size, ana, ana); analytic (size, ana, ana);
for (i = 0; i < size; i++) for (i = 0; i < size; i++)

6
wdsp/fir.hpp
View File

@ -35,9 +35,9 @@ class WDSP_API FIR
{ {
public: public:
static float* fftcv_mults (int NM, float* c_impulse); static float* fftcv_mults (int NM, float* c_impulse);
static float* fir_fsamp_odd (int N, float* A, int rtype, float scale, int wintype); static float* fir_fsamp_odd (int N, float* A, int rtype, double scale, int wintype);
static float* fir_fsamp (int N, float* A, int rtype, float scale, int wintype); static float* fir_fsamp (int N, float* A, int rtype, double scale, int wintype);
static float* fir_bandpass (int N, float f_low, float f_high, float samplerate, int wintype, int rtype, float scale); static float* fir_bandpass (int N, double f_low, double f_high, double samplerate, int wintype, int rtype, double scale);
static void mp_imp (int N, float* fir, float* mpfir, int pfactor, int polarity); static void mp_imp (int N, float* fir, float* mpfir, int pfactor, int polarity);
private: private:

8
wdsp/fircore.cpp
View File

@ -75,7 +75,13 @@ void FIRCORE::plan_fircore (FIRCORE *a)
FFTW_FORWARD, FFTW_FORWARD,
FFTW_PATIENT FFTW_PATIENT
); );
a->maskplan[1][i] = fftwf_plan_dft_1d(2 * a->size, (fftwf_complex *)a->maskgen, (fftwf_complex *)a->fmask[1][i], FFTW_FORWARD, FFTW_PATIENT); a->maskplan[1][i] = fftwf_plan_dft_1d(
2 * a->size,
(fftwf_complex *)a->maskgen,
(fftwf_complex *)a->fmask[1][i],
FFTW_FORWARD,
FFTW_PATIENT
);
} }
a->accum = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); a->accum = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex));
a->crev = fftwf_plan_dft_1d( a->crev = fftwf_plan_dft_1d(