From 531e96de00dafe9d0326229965d02a1e4c30c807 Mon Sep 17 00:00:00 2001 From: f4exb Date: Sat, 13 Jul 2024 03:20:26 +0200 Subject: [PATCH] WDSP: split firmin into fircore, firmin and firopt --- wdsp/CMakeLists.txt | 4 + wdsp/RXA.cpp | 2 +- wdsp/bandpass.cpp | 142 ++++++++++++++--- wdsp/bpsnba.cpp | 2 +- wdsp/cfir.cpp | 2 +- wdsp/emph.cpp | 2 +- wdsp/eq.cpp | 2 +- wdsp/fircore.cpp | 236 +++++++++++++++++++++++++++++ wdsp/fircore.hpp | 88 +++++++++++ wdsp/firmin.cpp | 360 -------------------------------------------- wdsp/firmin.hpp | 122 --------------- wdsp/firopt.cpp | 193 ++++++++++++++++++++++++ wdsp/firopt.hpp | 86 +++++++++++ wdsp/fmd.cpp | 2 +- wdsp/fmmod.cpp | 2 +- wdsp/fmsq.cpp | 2 +- wdsp/icfir.cpp | 2 +- wdsp/nbp.cpp | 2 +- wdsp/patchpanel.cpp | 12 +- wdsp/patchpanel.hpp | 12 +- 20 files changed, 759 insertions(+), 516 deletions(-) create mode 100644 wdsp/fircore.cpp create mode 100644 wdsp/fircore.hpp create mode 100644 wdsp/firopt.cpp create mode 100644 wdsp/firopt.hpp diff --git a/wdsp/CMakeLists.txt b/wdsp/CMakeLists.txt index 060b11a86..14530dea4 100644 --- a/wdsp/CMakeLists.txt +++ b/wdsp/CMakeLists.txt @@ -23,7 +23,9 @@ set(wdsp_SOURCES eq.cpp fcurve.cpp fir.cpp + fircore.cpp firmin.cpp + firopt.cpp fmd.cpp fmmod.cpp fmsq.cpp @@ -77,7 +79,9 @@ set(wdsp_HEADERS eq.hpp fcurve.hpp fir.hpp + fircore.hpp firmin.hpp + firopt.hpp fmd.hpp fmmod.hpp fmsq.hpp diff --git a/wdsp/RXA.cpp b/wdsp/RXA.cpp index fe9475ac5..405be0a73 100644 --- a/wdsp/RXA.cpp +++ b/wdsp/RXA.cpp @@ -51,7 +51,7 @@ warren@wpratt.com #include "cblock.hpp" #include "ssql.hpp" #include "iir.hpp" -#include "firmin.hpp" +#include "fircore.hpp" #include "wcpAGC.hpp" #include "anb.hpp" #include "nob.hpp" diff --git a/wdsp/bandpass.cpp b/wdsp/bandpass.cpp index 26982bc74..d8ac236c2 100644 --- a/wdsp/bandpass.cpp +++ b/wdsp/bandpass.cpp @@ -28,7 +28,7 @@ warren@wpratt.com #include "comm.hpp" #include "bandpass.hpp" #include "fir.hpp" -#include "firmin.hpp" +#include "fircore.hpp" #include "RXA.hpp" #include "TXA.hpp" @@ -40,12 +40,23 @@ namespace WDSP { * * ********************************************************************************************************/ -BANDPASS* BANDPASS::create_bandpass (int run, int position, int size, int nc, int mp, float* in, float* out, - float f_low, float f_high, int samplerate, int wintype, float gain) +BANDPASS* BANDPASS::create_bandpass ( + int run, + int position, + int size, + int nc, + int mp, + float* in, + float* out, + float f_low, + float f_high, + int samplerate, + int wintype, + float gain +) { // NOTE: 'nc' must be >= 'size' BANDPASS *a = new BANDPASS; - float* impulse; a->run = run; a->position = position; a->size = size; @@ -58,7 +69,15 @@ BANDPASS* BANDPASS::create_bandpass (int run, int position, int size, int nc, in a->samplerate = samplerate; a->wintype = wintype; a->gain = gain; - impulse = FIR::fir_bandpass (a->nc, a->f_low, a->f_high, a->samplerate, a->wintype, 1, a->gain / (float)(2 * a->size)); + float* impulse = FIR::fir_bandpass ( + a->nc, + a->f_low, + a->f_high, + a->samplerate, + a->wintype, + 1, + a->gain / (float)(2 * a->size) + ); a->p = FIRCORE::create_fircore (a->size, a->in, a->out, a->nc, a->mp, impulse); delete[] impulse; return a; @@ -92,9 +111,16 @@ void BANDPASS::setBuffers_bandpass (BANDPASS *a, float* in, float* out) void BANDPASS::setSamplerate_bandpass (BANDPASS *a, int rate) { - float* impulse; a->samplerate = rate; - impulse = FIR::fir_bandpass (a->nc, a->f_low, a->f_high, a->samplerate, a->wintype, 1, a->gain / (float)(2 * a->size)); + float* impulse = FIR::fir_bandpass ( + a->nc, + a->f_low, + a->f_high, + a->samplerate, + a->wintype, + 1, + a->gain / (float)(2 * a->size) + ); FIRCORE::setImpulse_fircore (a->p, impulse, 1); delete[] impulse; } @@ -102,33 +128,54 @@ void BANDPASS::setSamplerate_bandpass (BANDPASS *a, int rate) void BANDPASS::setSize_bandpass (BANDPASS *a, int size) { // NOTE: 'size' must be <= 'nc' - float* impulse; a->size = size; FIRCORE::setSize_fircore (a->p, a->size); // recalc impulse because scale factor is a function of size - impulse = FIR::fir_bandpass (a->nc, a->f_low, a->f_high, a->samplerate, a->wintype, 1, a->gain / (float)(2 * a->size)); + float* impulse = FIR::fir_bandpass ( + a->nc, + a->f_low, + a->f_high, + a->samplerate, + a->wintype, + 1, + a->gain / (float)(2 * a->size) + ); FIRCORE::setImpulse_fircore (a->p, impulse, 1); delete[] (impulse); } void BANDPASS::setGain_bandpass (BANDPASS *a, float gain, int update) { - float* impulse; a->gain = gain; - impulse = FIR::fir_bandpass (a->nc, a->f_low, a->f_high, a->samplerate, a->wintype, 1, a->gain / (float)(2 * a->size)); + float* impulse = FIR::fir_bandpass ( + a->nc, + a->f_low, + a->f_high, + a->samplerate, + a->wintype, + 1, + a->gain / (float)(2 * a->size) + ); FIRCORE::setImpulse_fircore (a->p, impulse, update); delete[] (impulse); } void BANDPASS::CalcBandpassFilter (BANDPASS *a, float f_low, float f_high, float gain) { - float* impulse; if ((a->f_low != f_low) || (a->f_high != f_high) || (a->gain != gain)) { a->f_low = f_low; a->f_high = f_high; a->gain = gain; - impulse = FIR::fir_bandpass (a->nc, a->f_low, a->f_high, a->samplerate, a->wintype, 1, a->gain / (float)(2 * a->size)); + float* impulse = FIR::fir_bandpass ( + a->nc, + a->f_low, + a->f_high, + a->samplerate, + a->wintype, + 1, + a->gain / (float)(2 * a->size) + ); FIRCORE::setImpulse_fircore (a->p, impulse, 1); delete[] (impulse); } @@ -142,12 +189,19 @@ void BANDPASS::CalcBandpassFilter (BANDPASS *a, float f_low, float f_high, float void BANDPASS::SetBandpassFreqs (RXA& rxa, float f_low, float f_high) { - float* impulse; BANDPASS *a = rxa.bp1.p; + if ((f_low != a->f_low) || (f_high != a->f_high)) { - impulse = FIR::fir_bandpass (a->nc, f_low, f_high, a->samplerate, - a->wintype, 1, a->gain / (float)(2 * a->size)); + float* impulse = FIR::fir_bandpass ( + a->nc, + f_low, + f_high, + a->samplerate, + a->wintype, + 1, + a->gain / (float)(2 * a->size) + ); FIRCORE::setImpulse_fircore (a->p, impulse, 0); delete[] (impulse); rxa.csDSP.lock(); @@ -161,17 +215,26 @@ void BANDPASS::SetBandpassFreqs (RXA& rxa, float f_low, float f_high) void BANDPASS::SetBandpassNC (RXA& rxa, int nc) { // NOTE: 'nc' must be >= 'size' - float* impulse; BANDPASS *a; rxa.csDSP.lock(); a = rxa.bp1.p; + if (nc != a->nc) { a->nc = nc; - impulse = FIR::fir_bandpass (a->nc, a->f_low, a->f_high, a->samplerate, a->wintype, 1, a->gain / (float)(2 * a->size)); + float* impulse = FIR::fir_bandpass ( + a->nc, + a->f_low, + a->f_high, + a->samplerate, + a->wintype, + 1, + a->gain / (float)(2 * a->size) + ); FIRCORE::setNc_fircore (a->p, a->nc, impulse); delete[] (impulse); } + rxa.csDSP.unlock(); } @@ -179,6 +242,7 @@ void BANDPASS::SetBandpassMP (RXA& rxa, int mp) { BANDPASS *a; a = rxa.bp1.p; + if (mp != a->mp) { a->mp = mp; @@ -229,33 +293,62 @@ void BANDPASS::SetBandpassMP (RXA& rxa, int mp) void BANDPASS::SetBandpassNC (TXA& txa, int nc) { // NOTE: 'nc' must be >= 'size' - float* impulse; BANDPASS *a; txa.csDSP.lock(); a = txa.bp0.p; + if (a->nc != nc) { a->nc = nc; - impulse = FIR::fir_bandpass (a->nc, a->f_low, a->f_high, a->samplerate, a->wintype, 1, a->gain / (float)(2 * a->size)); + float* impulse = FIR::fir_bandpass ( + a->nc, + a->f_low, + a->f_high, + a->samplerate, + a->wintype, + 1, + a->gain / (float)(2 * a->size) + ); FIRCORE::setNc_fircore (a->p, a->nc, impulse); delete[] (impulse); } + a = txa.bp1.p; + if (a->nc != nc) { a->nc = nc; - impulse = FIR::fir_bandpass (a->nc, a->f_low, a->f_high, a->samplerate, a->wintype, 1, a->gain / (float)(2 * a->size)); + float* impulse = FIR::fir_bandpass ( + a->nc, + a->f_low, + a->f_high, + a->samplerate, + a->wintype, + 1, + a->gain / (float)(2 * a->size) + ); FIRCORE::setNc_fircore (a->p, a->nc, impulse); delete[] (impulse); } + a = txa.bp2.p; + if (a->nc != nc) { a->nc = nc; - impulse = FIR::fir_bandpass (a->nc, a->f_low, a->f_high, a->samplerate, a->wintype, 1, a->gain / (float)(2 * a->size)); + float* impulse = FIR::fir_bandpass ( + a->nc, + a->f_low, + a->f_high, + a->samplerate, + a->wintype, + 1, + a->gain / (float)(2 * a->size) + ); FIRCORE::setNc_fircore (a->p, a->nc, impulse); delete[] (impulse); } + txa.csDSP.unlock(); } @@ -263,18 +356,23 @@ void BANDPASS::SetBandpassMP (TXA& txa, int mp) { BANDPASS *a; a = txa.bp0.p; + if (mp != a->mp) { a->mp = mp; FIRCORE::setMp_fircore (a->p, a->mp); } + a = txa.bp1.p; + if (mp != a->mp) { a->mp = mp; FIRCORE::setMp_fircore (a->p, a->mp); } + a = txa.bp2.p; + if (mp != a->mp) { a->mp = mp; diff --git a/wdsp/bpsnba.cpp b/wdsp/bpsnba.cpp index 540609ba4..87ab4dc82 100644 --- a/wdsp/bpsnba.cpp +++ b/wdsp/bpsnba.cpp @@ -28,7 +28,7 @@ warren@wpratt.com #include "comm.hpp" #include "resample.hpp" #include "lmath.hpp" -#include "firmin.hpp" +#include "fircore.hpp" #include "nbp.hpp" #include "amd.hpp" #include "anf.hpp" diff --git a/wdsp/cfir.cpp b/wdsp/cfir.cpp index 64ebb5ad6..554d055fb 100644 --- a/wdsp/cfir.cpp +++ b/wdsp/cfir.cpp @@ -27,7 +27,7 @@ warren@wpratt.com #include "comm.hpp" #include "cfir.hpp" #include "fir.hpp" -#include "firmin.hpp" +#include "fircore.hpp" #include "TXA.hpp" namespace WDSP { diff --git a/wdsp/emph.cpp b/wdsp/emph.cpp index fc31dafd1..9a8a9026d 100644 --- a/wdsp/emph.cpp +++ b/wdsp/emph.cpp @@ -28,7 +28,7 @@ warren@wpratt.com #include "comm.hpp" #include "emph.hpp" #include "fcurve.hpp" -#include "firmin.hpp" +#include "fircore.hpp" #include "TXA.hpp" namespace WDSP { diff --git a/wdsp/eq.cpp b/wdsp/eq.cpp index f627484e5..4cb591b84 100644 --- a/wdsp/eq.cpp +++ b/wdsp/eq.cpp @@ -27,7 +27,7 @@ warren@wpratt.com #include "comm.hpp" #include "eq.hpp" -#include "firmin.hpp" +#include "fircore.hpp" #include "fir.hpp" #include "RXA.hpp" #include "TXA.hpp" diff --git a/wdsp/fircore.cpp b/wdsp/fircore.cpp new file mode 100644 index 000000000..65bbe1361 --- /dev/null +++ b/wdsp/fircore.cpp @@ -0,0 +1,236 @@ +/* firmin.c + +This file is part of a program that implements a Software-Defined Radio. + +Copyright (C) 2016 Warren Pratt, NR0V +Copyright (C) 2024 Edouard Griffiths, F4EXB Adapted to SDRangel + +This program is free software; you can redistribute it and/or +modify it under the terms of the GNU General Public License +as published by the Free Software Foundation; either version 2 +of the License, or (at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + +The author can be reached by email at + +warren@wpratt.com + +*/ + +#include "comm.hpp" +#include "fir.hpp" +#include "fircore.hpp" + +namespace WDSP { + +/******************************************************************************************************** +* * +* Partitioned Overlap-Save Filter Kernel * +* * +********************************************************************************************************/ + + +void FIRCORE::plan_fircore (FIRCORE *a) +{ + // must call for change in 'nc', 'size', 'out' + int i; + a->nfor = a->nc / a->size; + a->cset = 0; + a->buffidx = 0; + a->idxmask = a->nfor - 1; + a->fftin = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); + a->fftout = new float*[a->nfor]; // (float **) malloc0 (a->nfor * sizeof (float *)); + a->fmask = new float**[2]; // (float ***) malloc0 (2 * sizeof (float **)); + a->fmask[0] = new float*[a->nfor]; // (float **) malloc0 (a->nfor * sizeof (float *)); + a->fmask[1] = new float*[a->nfor]; // (float **) malloc0 (a->nfor * sizeof (float *)); + a->maskgen = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); + a->pcfor = new fftwf_plan[a->nfor]; // (fftwf_plan *) malloc0 (a->nfor * sizeof (fftwf_plan)); + a->maskplan = new fftwf_plan*[2]; // (fftwf_plan **) malloc0 (2 * sizeof (fftwf_plan *)); + a->maskplan[0] = new fftwf_plan[a->nfor]; // (fftwf_plan *) malloc0 (a->nfor * sizeof (fftwf_plan)); + a->maskplan[1] = new fftwf_plan[a->nfor]; // (fftwf_plan *) malloc0 (a->nfor * sizeof (fftwf_plan)); + for (i = 0; i < a->nfor; i++) + { + a->fftout[i] = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); + a->fmask[0][i] = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); + a->fmask[1][i] = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); + a->pcfor[i] = fftwf_plan_dft_1d(2 * a->size, (fftwf_complex *)a->fftin, (fftwf_complex *)a->fftout[i], FFTW_FORWARD, FFTW_PATIENT); + a->maskplan[0][i] = fftwf_plan_dft_1d(2 * a->size, (fftwf_complex *)a->maskgen, (fftwf_complex *)a->fmask[0][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->crev = fftwf_plan_dft_1d(2 * a->size, (fftwf_complex *)a->accum, (fftwf_complex *)a->out, FFTW_BACKWARD, FFTW_PATIENT); + a->masks_ready = 0; +} + +void FIRCORE::calc_fircore (FIRCORE *a, int flip) +{ + // call for change in frequency, rate, wintype, gain + // must also call after a call to plan_firopt() + int i; + if (a->mp) + FIR::mp_imp (a->nc, a->impulse, a->imp, 16, 0); + else + memcpy (a->imp, a->impulse, a->nc * sizeof (wcomplex)); + for (i = 0; i < a->nfor; i++) + { + // I right-justified the impulse response => take output from left side of output buff, discard right side + // Be careful about flipping an asymmetrical impulse response. + memcpy (&(a->maskgen[2 * a->size]), &(a->imp[2 * a->size * i]), a->size * sizeof(wcomplex)); + fftwf_execute (a->maskplan[1 - a->cset][i]); + } + a->masks_ready = 1; + if (flip) + { + a->update.lock(); + a->cset = 1 - a->cset; + a->update.unlock(); + a->masks_ready = 0; + } +} + +FIRCORE* FIRCORE::create_fircore (int size, float* in, float* out, int nc, int mp, float* impulse) +{ + FIRCORE *a = new FIRCORE; + a->size = size; + a->in = in; + a->out = out; + a->nc = nc; + a->mp = mp; + // InitializeCriticalSectionAndSpinCount (&a->update, 2500); + plan_fircore (a); + a->impulse = new float[a->nc * 2]; // (float *) malloc0 (a->nc * sizeof (complex)); + a->imp = new float[a->nc * 2]; // (float *) malloc0 (a->nc * sizeof (complex)); + memcpy (a->impulse, impulse, a->nc * sizeof (wcomplex)); + calc_fircore (a, 1); + return a; +} + +void FIRCORE::deplan_fircore (FIRCORE *a) +{ + int i; + fftwf_destroy_plan (a->crev); + delete[] (a->accum); + for (i = 0; i < a->nfor; i++) + { + delete[] (a->fftout[i]); + delete[] (a->fmask[0][i]); + delete[] (a->fmask[1][i]); + fftwf_destroy_plan (a->pcfor[i]); + fftwf_destroy_plan (a->maskplan[0][i]); + fftwf_destroy_plan (a->maskplan[1][i]); + } + delete[] (a->maskplan[0]); + delete[] (a->maskplan[1]); + delete[] (a->maskplan); + delete[] (a->pcfor); + delete[] (a->maskgen); + delete[] (a->fmask[0]); + delete[] (a->fmask[1]); + delete[] (a->fmask); + delete[] (a->fftout); + delete[] (a->fftin); +} + +void FIRCORE::destroy_fircore (FIRCORE *a) +{ + deplan_fircore (a); + delete[] (a->imp); + delete[] (a->impulse); + delete (a); +} + +void FIRCORE::flush_fircore (FIRCORE *a) +{ + int i; + memset (a->fftin, 0, 2 * a->size * sizeof (wcomplex)); + for (i = 0; i < a->nfor; i++) + memset (a->fftout[i], 0, 2 * a->size * sizeof (wcomplex)); + a->buffidx = 0; +} + +void FIRCORE::xfircore (FIRCORE *a) +{ + int i, j, k; + memcpy (&(a->fftin[2 * a->size]), a->in, a->size * sizeof (wcomplex)); + fftwf_execute (a->pcfor[a->buffidx]); + k = a->buffidx; + memset (a->accum, 0, 2 * a->size * sizeof (wcomplex)); + a->update.lock(); + for (j = 0; j < a->nfor; j++) + { + for (i = 0; i < 2 * a->size; i++) + { + a->accum[2 * i + 0] += a->fftout[k][2 * i + 0] * a->fmask[a->cset][j][2 * i + 0] - a->fftout[k][2 * i + 1] * a->fmask[a->cset][j][2 * i + 1]; + a->accum[2 * i + 1] += a->fftout[k][2 * i + 0] * a->fmask[a->cset][j][2 * i + 1] + a->fftout[k][2 * i + 1] * a->fmask[a->cset][j][2 * i + 0]; + } + k = (k + a->idxmask) & a->idxmask; + } + a->update.unlock(); + a->buffidx = (a->buffidx + 1) & a->idxmask; + fftwf_execute (a->crev); + memcpy (a->fftin, &(a->fftin[2 * a->size]), a->size * sizeof(wcomplex)); +} + +void FIRCORE::setBuffers_fircore (FIRCORE *a, float* in, float* out) +{ + a->in = in; + a->out = out; + deplan_fircore (a); + plan_fircore (a); + calc_fircore (a, 1); +} + +void FIRCORE::setSize_fircore (FIRCORE *a, int size) +{ + a->size = size; + deplan_fircore (a); + plan_fircore (a); + calc_fircore (a, 1); +} + +void FIRCORE::setImpulse_fircore (FIRCORE *a, float* impulse, int update) +{ + memcpy (a->impulse, impulse, a->nc * sizeof (wcomplex)); + calc_fircore (a, update); +} + +void FIRCORE::setNc_fircore (FIRCORE *a, int nc, float* impulse) +{ + // because of FFT planning, this will probably cause a glitch in audio if done during dataflow + deplan_fircore (a); + delete[] (a->impulse); + delete[] (a->imp); + a->nc = nc; + plan_fircore (a); + a->imp = new float[a->nc * 2]; // (float *) malloc0 (a->nc * sizeof (complex)); + a->impulse = new float[a->nc * 2]; // (float *) malloc0 (a->nc * sizeof (complex)); + memcpy (a->impulse, impulse, a->nc * sizeof (wcomplex)); + calc_fircore (a, 1); +} + +void FIRCORE::setMp_fircore (FIRCORE *a, int mp) +{ + a->mp = mp; + calc_fircore (a, 1); +} + +void FIRCORE::setUpdate_fircore (FIRCORE *a) +{ + if (a->masks_ready) + { + a->update.lock(); + a->cset = 1 - a->cset; + a->update.unlock(); + a->masks_ready = 0; + } +} + +} // namespace WDSP diff --git a/wdsp/fircore.hpp b/wdsp/fircore.hpp new file mode 100644 index 000000000..b57aa2c06 --- /dev/null +++ b/wdsp/fircore.hpp @@ -0,0 +1,88 @@ +/* firmin.h + +This file is part of a program that implements a Software-Defined Radio. + +Copyright (C) 2016 Warren Pratt, NR0V +Copyright (C) 2024 Edouard Griffiths, F4EXB Adapted to SDRangel + +This program is free software; you can redistribute it and/or +modify it under the terms of the GNU General Public License +as published by the Free Software Foundation; either version 2 +of the License, or (at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + +The author can be reached by email at + +warren@wpratt.com + +*/ + +/******************************************************************************************************** +* * +* Partitioned Overlap-Save Filter Kernel * +* * +********************************************************************************************************/ + +#ifndef wdsp_fircore_h +#define wdsp_fircore_h + +#include +#include "fftw3.h" +#include "export.h" + +namespace WDSP { + +class WDSP_API FIRCORE +{ +public: + int size; // input/output buffer size, power of two + float* in; // input buffer + float* out; // output buffer, can be same as input + int nc; // number of filter coefficients, power of two, >= size + float* impulse; // impulse response of filter + float* imp; + int nfor; // number of buffers in delay line + float* fftin; // fft input buffer + float*** fmask; // frequency domain masks + float** fftout; // fftout delay line + float* accum; // frequency domain accumulator + int buffidx; // fft out buffer index + int idxmask; // mask for index computations + float* maskgen; // input for mask generation FFT + fftwf_plan* pcfor; // array of forward FFT plans + fftwf_plan crev; // reverse fft plan + fftwf_plan** maskplan; // plans for frequency domain masks + QRecursiveMutex update; + int cset; + int mp; + int masks_ready; + + static FIRCORE* create_fircore (int size, float* in, float* out, + int nc, int mp, float* impulse); + static void xfircore (FIRCORE *a); + static void destroy_fircore (FIRCORE *a); + static void flush_fircore (FIRCORE *a); + static void setBuffers_fircore (FIRCORE *a, float* in, float* out); + static void setSize_fircore (FIRCORE *a, int size); + static void setImpulse_fircore (FIRCORE *a, float* impulse, int update); + static void setNc_fircore (FIRCORE *a, int nc, float* impulse); + static void setMp_fircore (FIRCORE *a, int mp); + static void setUpdate_fircore (FIRCORE *a); + +private: + static void plan_fircore (FIRCORE *a); + static void calc_fircore (FIRCORE *a, int flip); + static void deplan_fircore (FIRCORE *a); +}; + +} // namespace WDSP + +#endif diff --git a/wdsp/firmin.cpp b/wdsp/firmin.cpp index 235d4f60b..d11c4c8e1 100644 --- a/wdsp/firmin.cpp +++ b/wdsp/firmin.cpp @@ -127,364 +127,4 @@ void FIRMIN::setFreqs_firmin (FIRMIN *a, float f_low, float f_high) calc_firmin (a); } -/******************************************************************************************************** -* * -* Standalone Partitioned Overlap-Save Bandpass * -* * -********************************************************************************************************/ - -void FIROPT::plan_firopt (FIROPT *a) -{ - // must call for change in 'nc', 'size', 'out' - int i; - a->nfor = a->nc / a->size; - a->buffidx = 0; - a->idxmask = a->nfor - 1; - a->fftin = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); - a->fftout = new float*[a->nfor]; // (float **) malloc0 (a->nfor * sizeof (float *)); - a->fmask = new float*[a->nfor]; // (float **) malloc0 (a->nfor * sizeof (float *)); - a->maskgen = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); - a->pcfor = new fftwf_plan[a->nfor]; // (fftwf_plan *) malloc0 (a->nfor * sizeof (fftwf_plan)); - a->maskplan = new fftwf_plan[a->nfor]; // (fftwf_plan *) malloc0 (a->nfor * sizeof (fftwf_plan)); - for (i = 0; i < a->nfor; i++) - { - a->fftout[i] = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); - a->fmask[i] = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); - a->pcfor[i] = fftwf_plan_dft_1d(2 * a->size, (fftwf_complex *)a->fftin, (fftwf_complex *)a->fftout[i], FFTW_FORWARD, FFTW_PATIENT); - a->maskplan[i] = fftwf_plan_dft_1d(2 * a->size, (fftwf_complex *)a->maskgen, (fftwf_complex *)a->fmask[i], FFTW_FORWARD, FFTW_PATIENT); - } - a->accum = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); - a->crev = fftwf_plan_dft_1d(2 * a->size, (fftwf_complex *)a->accum, (fftwf_complex *)a->out, FFTW_BACKWARD, FFTW_PATIENT); -} - -void FIROPT::calc_firopt (FIROPT *a) -{ - // call for change in frequency, rate, wintype, gain - // must also call after a call to plan_firopt() - int i; - float* impulse = FIR::fir_bandpass (a->nc, a->f_low, a->f_high, a->samplerate, a->wintype, 1, a->gain); - a->buffidx = 0; - for (i = 0; i < a->nfor; i++) - { - // I right-justified the impulse response => take output from left side of output buff, discard right side - // Be careful about flipping an asymmetrical impulse response. - memcpy (&(a->maskgen[2 * a->size]), &(impulse[2 * a->size * i]), a->size * sizeof(wcomplex)); - fftwf_execute (a->maskplan[i]); - } - delete[] (impulse); -} - -FIROPT* FIROPT::create_firopt (int run, int position, int size, float* in, float* out, - int nc, float f_low, float f_high, int samplerate, int wintype, float gain) -{ - FIROPT *a = new FIROPT; - a->run = run; - a->position = position; - a->size = size; - a->in = in; - a->out = out; - a->nc = nc; - a->f_low = f_low; - a->f_high = f_high; - a->samplerate = samplerate; - a->wintype = wintype; - a->gain = gain; - plan_firopt (a); - calc_firopt (a); - return a; -} - -void FIROPT::deplan_firopt (FIROPT *a) -{ - int i; - fftwf_destroy_plan (a->crev); - delete[] (a->accum); - for (i = 0; i < a->nfor; i++) - { - delete[] (a->fftout[i]); - delete[] (a->fmask[i]); - fftwf_destroy_plan (a->pcfor[i]); - fftwf_destroy_plan (a->maskplan[i]); - } - delete[] (a->maskplan); - delete[] (a->pcfor); - delete[] (a->maskgen); - delete[] (a->fmask); - delete[] (a->fftout); - delete[] (a->fftin); -} - -void FIROPT::destroy_firopt (FIROPT *a) -{ - deplan_firopt (a); - delete (a); -} - -void FIROPT::flush_firopt (FIROPT *a) -{ - int i; - memset (a->fftin, 0, 2 * a->size * sizeof (wcomplex)); - for (i = 0; i < a->nfor; i++) - memset (a->fftout[i], 0, 2 * a->size * sizeof (wcomplex)); - a->buffidx = 0; -} - -void FIROPT::xfiropt (FIROPT *a, int pos) -{ - if (a->run && (a->position == pos)) - { - int i, j, k; - memcpy (&(a->fftin[2 * a->size]), a->in, a->size * sizeof (wcomplex)); - fftwf_execute (a->pcfor[a->buffidx]); - k = a->buffidx; - memset (a->accum, 0, 2 * a->size * sizeof (wcomplex)); - for (j = 0; j < a->nfor; j++) - { - for (i = 0; i < 2 * a->size; i++) - { - a->accum[2 * i + 0] += a->fftout[k][2 * i + 0] * a->fmask[j][2 * i + 0] - a->fftout[k][2 * i + 1] * a->fmask[j][2 * i + 1]; - a->accum[2 * i + 1] += a->fftout[k][2 * i + 0] * a->fmask[j][2 * i + 1] + a->fftout[k][2 * i + 1] * a->fmask[j][2 * i + 0]; - } - k = (k + a->idxmask) & a->idxmask; - } - a->buffidx = (a->buffidx + 1) & a->idxmask; - fftwf_execute (a->crev); - memcpy (a->fftin, &(a->fftin[2 * a->size]), a->size * sizeof(wcomplex)); - } - else if (a->in != a->out) - memcpy (a->out, a->in, a->size * sizeof (wcomplex)); -} - -void FIROPT::setBuffers_firopt (FIROPT *a, float* in, float* out) -{ - a->in = in; - a->out = out; - deplan_firopt (a); - plan_firopt (a); - calc_firopt (a); -} - -void FIROPT::setSamplerate_firopt (FIROPT *a, int rate) -{ - a->samplerate = rate; - calc_firopt (a); -} - -void FIROPT::setSize_firopt (FIROPT *a, int size) -{ - a->size = size; - deplan_firopt (a); - plan_firopt (a); - calc_firopt (a); -} - -void FIROPT::setFreqs_firopt (FIROPT *a, float f_low, float f_high) -{ - a->f_low = f_low; - a->f_high = f_high; - calc_firopt (a); -} - -/******************************************************************************************************** -* * -* Partitioned Overlap-Save Filter Kernel * -* * -********************************************************************************************************/ - - -void FIRCORE::plan_fircore (FIRCORE *a) -{ - // must call for change in 'nc', 'size', 'out' - int i; - a->nfor = a->nc / a->size; - a->cset = 0; - a->buffidx = 0; - a->idxmask = a->nfor - 1; - a->fftin = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); - a->fftout = new float*[a->nfor]; // (float **) malloc0 (a->nfor * sizeof (float *)); - a->fmask = new float**[2]; // (float ***) malloc0 (2 * sizeof (float **)); - a->fmask[0] = new float*[a->nfor]; // (float **) malloc0 (a->nfor * sizeof (float *)); - a->fmask[1] = new float*[a->nfor]; // (float **) malloc0 (a->nfor * sizeof (float *)); - a->maskgen = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); - a->pcfor = new fftwf_plan[a->nfor]; // (fftwf_plan *) malloc0 (a->nfor * sizeof (fftwf_plan)); - a->maskplan = new fftwf_plan*[2]; // (fftwf_plan **) malloc0 (2 * sizeof (fftwf_plan *)); - a->maskplan[0] = new fftwf_plan[a->nfor]; // (fftwf_plan *) malloc0 (a->nfor * sizeof (fftwf_plan)); - a->maskplan[1] = new fftwf_plan[a->nfor]; // (fftwf_plan *) malloc0 (a->nfor * sizeof (fftwf_plan)); - for (i = 0; i < a->nfor; i++) - { - a->fftout[i] = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); - a->fmask[0][i] = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); - a->fmask[1][i] = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); - a->pcfor[i] = fftwf_plan_dft_1d(2 * a->size, (fftwf_complex *)a->fftin, (fftwf_complex *)a->fftout[i], FFTW_FORWARD, FFTW_PATIENT); - a->maskplan[0][i] = fftwf_plan_dft_1d(2 * a->size, (fftwf_complex *)a->maskgen, (fftwf_complex *)a->fmask[0][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->crev = fftwf_plan_dft_1d(2 * a->size, (fftwf_complex *)a->accum, (fftwf_complex *)a->out, FFTW_BACKWARD, FFTW_PATIENT); - a->masks_ready = 0; -} - -void FIRCORE::calc_fircore (FIRCORE *a, int flip) -{ - // call for change in frequency, rate, wintype, gain - // must also call after a call to plan_firopt() - int i; - if (a->mp) - FIR::mp_imp (a->nc, a->impulse, a->imp, 16, 0); - else - memcpy (a->imp, a->impulse, a->nc * sizeof (wcomplex)); - for (i = 0; i < a->nfor; i++) - { - // I right-justified the impulse response => take output from left side of output buff, discard right side - // Be careful about flipping an asymmetrical impulse response. - memcpy (&(a->maskgen[2 * a->size]), &(a->imp[2 * a->size * i]), a->size * sizeof(wcomplex)); - fftwf_execute (a->maskplan[1 - a->cset][i]); - } - a->masks_ready = 1; - if (flip) - { - a->update.lock(); - a->cset = 1 - a->cset; - a->update.unlock(); - a->masks_ready = 0; - } -} - -FIRCORE* FIRCORE::create_fircore (int size, float* in, float* out, int nc, int mp, float* impulse) -{ - FIRCORE *a = new FIRCORE; - a->size = size; - a->in = in; - a->out = out; - a->nc = nc; - a->mp = mp; - // InitializeCriticalSectionAndSpinCount (&a->update, 2500); - plan_fircore (a); - a->impulse = new float[a->nc * 2]; // (float *) malloc0 (a->nc * sizeof (complex)); - a->imp = new float[a->nc * 2]; // (float *) malloc0 (a->nc * sizeof (complex)); - memcpy (a->impulse, impulse, a->nc * sizeof (wcomplex)); - calc_fircore (a, 1); - return a; -} - -void FIRCORE::deplan_fircore (FIRCORE *a) -{ - int i; - fftwf_destroy_plan (a->crev); - delete[] (a->accum); - for (i = 0; i < a->nfor; i++) - { - delete[] (a->fftout[i]); - delete[] (a->fmask[0][i]); - delete[] (a->fmask[1][i]); - fftwf_destroy_plan (a->pcfor[i]); - fftwf_destroy_plan (a->maskplan[0][i]); - fftwf_destroy_plan (a->maskplan[1][i]); - } - delete[] (a->maskplan[0]); - delete[] (a->maskplan[1]); - delete[] (a->maskplan); - delete[] (a->pcfor); - delete[] (a->maskgen); - delete[] (a->fmask[0]); - delete[] (a->fmask[1]); - delete[] (a->fmask); - delete[] (a->fftout); - delete[] (a->fftin); -} - -void FIRCORE::destroy_fircore (FIRCORE *a) -{ - deplan_fircore (a); - delete[] (a->imp); - delete[] (a->impulse); - delete (a); -} - -void FIRCORE::flush_fircore (FIRCORE *a) -{ - int i; - memset (a->fftin, 0, 2 * a->size * sizeof (wcomplex)); - for (i = 0; i < a->nfor; i++) - memset (a->fftout[i], 0, 2 * a->size * sizeof (wcomplex)); - a->buffidx = 0; -} - -void FIRCORE::xfircore (FIRCORE *a) -{ - int i, j, k; - memcpy (&(a->fftin[2 * a->size]), a->in, a->size * sizeof (wcomplex)); - fftwf_execute (a->pcfor[a->buffidx]); - k = a->buffidx; - memset (a->accum, 0, 2 * a->size * sizeof (wcomplex)); - a->update.lock(); - for (j = 0; j < a->nfor; j++) - { - for (i = 0; i < 2 * a->size; i++) - { - a->accum[2 * i + 0] += a->fftout[k][2 * i + 0] * a->fmask[a->cset][j][2 * i + 0] - a->fftout[k][2 * i + 1] * a->fmask[a->cset][j][2 * i + 1]; - a->accum[2 * i + 1] += a->fftout[k][2 * i + 0] * a->fmask[a->cset][j][2 * i + 1] + a->fftout[k][2 * i + 1] * a->fmask[a->cset][j][2 * i + 0]; - } - k = (k + a->idxmask) & a->idxmask; - } - a->update.unlock(); - a->buffidx = (a->buffidx + 1) & a->idxmask; - fftwf_execute (a->crev); - memcpy (a->fftin, &(a->fftin[2 * a->size]), a->size * sizeof(wcomplex)); -} - -void FIRCORE::setBuffers_fircore (FIRCORE *a, float* in, float* out) -{ - a->in = in; - a->out = out; - deplan_fircore (a); - plan_fircore (a); - calc_fircore (a, 1); -} - -void FIRCORE::setSize_fircore (FIRCORE *a, int size) -{ - a->size = size; - deplan_fircore (a); - plan_fircore (a); - calc_fircore (a, 1); -} - -void FIRCORE::setImpulse_fircore (FIRCORE *a, float* impulse, int update) -{ - memcpy (a->impulse, impulse, a->nc * sizeof (wcomplex)); - calc_fircore (a, update); -} - -void FIRCORE::setNc_fircore (FIRCORE *a, int nc, float* impulse) -{ - // because of FFT planning, this will probably cause a glitch in audio if done during dataflow - deplan_fircore (a); - delete[] (a->impulse); - delete[] (a->imp); - a->nc = nc; - plan_fircore (a); - a->imp = new float[a->nc * 2]; // (float *) malloc0 (a->nc * sizeof (complex)); - a->impulse = new float[a->nc * 2]; // (float *) malloc0 (a->nc * sizeof (complex)); - memcpy (a->impulse, impulse, a->nc * sizeof (wcomplex)); - calc_fircore (a, 1); -} - -void FIRCORE::setMp_fircore (FIRCORE *a, int mp) -{ - a->mp = mp; - calc_fircore (a, 1); -} - -void FIRCORE::setUpdate_fircore (FIRCORE *a) -{ - if (a->masks_ready) - { - a->update.lock(); - a->cset = 1 - a->cset; - a->update.unlock(); - a->masks_ready = 0; - } -} - } // namespace WDSP diff --git a/wdsp/firmin.hpp b/wdsp/firmin.hpp index effc44738..0e0e1b21b 100644 --- a/wdsp/firmin.hpp +++ b/wdsp/firmin.hpp @@ -76,125 +76,3 @@ private: } // namespace WDSP #endif - - -/******************************************************************************************************** -* * -* Standalone Partitioned Overlap-Save Bandpass * -* * -********************************************************************************************************/ - -#ifndef wdsp_firopt_h -#define wdsp_firopt_h - -#include "fftw3.h" -#include "export.h" - -namespace WDSP { - -class WDSP_API FIROPT -{ - int run; // run control - int position; // position at which to execute - int size; // input/output buffer size, power of two - float* in; // input buffer - float* out; // output buffer, can be same as input - int nc; // number of filter coefficients, power of two, >= size - float f_low; // low cutoff frequency - float f_high; // high cutoff frequency - float samplerate; // sample rate - int wintype; // filter window type - float gain; // filter gain - int nfor; // number of buffers in delay line - float* fftin; // fft input buffer - float** fmask; // frequency domain masks - float** fftout; // fftout delay line - float* accum; // frequency domain accumulator - int buffidx; // fft out buffer index - int idxmask; // mask for index computations - float* maskgen; // input for mask generation FFT - fftwf_plan* pcfor; // array of forward FFT plans - fftwf_plan crev; // reverse fft plan - fftwf_plan* maskplan; // plans for frequency domain masks - - static FIROPT* create_firopt (int run, int position, int size, float* in, float* out, - int nc, float f_low, float f_high, int samplerate, int wintype, float gain); - static void xfiropt (FIROPT *a, int pos); - static void destroy_firopt (FIROPT *a); - static void flush_firopt (FIROPT *a); - static void setBuffers_firopt (FIROPT *a, float* in, float* out); - static void setSamplerate_firopt (FIROPT *a, int rate); - static void setSize_firopt (FIROPT *a, int size); - static void setFreqs_firopt (FIROPT *a, float f_low, float f_high); - -private: - static void plan_firopt (FIROPT *a); - static void calc_firopt (FIROPT *a); - static void deplan_firopt (FIROPT *a); -}; - -} // namespace WDSP - -#endif - - -/******************************************************************************************************** -* * -* Partitioned Overlap-Save Filter Kernel * -* * -********************************************************************************************************/ - -#ifndef wdsp_fircore_h -#define wdsp_fircore_h - -#include -#include "export.h" - -namespace WDSP { - -class WDSP_API FIRCORE -{ -public: - int size; // input/output buffer size, power of two - float* in; // input buffer - float* out; // output buffer, can be same as input - int nc; // number of filter coefficients, power of two, >= size - float* impulse; // impulse response of filter - float* imp; - int nfor; // number of buffers in delay line - float* fftin; // fft input buffer - float*** fmask; // frequency domain masks - float** fftout; // fftout delay line - float* accum; // frequency domain accumulator - int buffidx; // fft out buffer index - int idxmask; // mask for index computations - float* maskgen; // input for mask generation FFT - fftwf_plan* pcfor; // array of forward FFT plans - fftwf_plan crev; // reverse fft plan - fftwf_plan** maskplan; // plans for frequency domain masks - QRecursiveMutex update; - int cset; - int mp; - int masks_ready; - - static FIRCORE* create_fircore (int size, float* in, float* out, - int nc, int mp, float* impulse); - static void xfircore (FIRCORE *a); - static void destroy_fircore (FIRCORE *a); - static void flush_fircore (FIRCORE *a); - static void setBuffers_fircore (FIRCORE *a, float* in, float* out); - static void setSize_fircore (FIRCORE *a, int size); - static void setImpulse_fircore (FIRCORE *a, float* impulse, int update); - static void setNc_fircore (FIRCORE *a, int nc, float* impulse); - static void setMp_fircore (FIRCORE *a, int mp); - static void setUpdate_fircore (FIRCORE *a); - -private: - static void plan_fircore (FIRCORE *a); - static void calc_fircore (FIRCORE *a, int flip); - static void deplan_fircore (FIRCORE *a); -}; - -} // namespace WDSP - -#endif diff --git a/wdsp/firopt.cpp b/wdsp/firopt.cpp new file mode 100644 index 000000000..73e451a3a --- /dev/null +++ b/wdsp/firopt.cpp @@ -0,0 +1,193 @@ +/* firmin.c + +This file is part of a program that implements a Software-Defined Radio. + +Copyright (C) 2016 Warren Pratt, NR0V +Copyright (C) 2024 Edouard Griffiths, F4EXB Adapted to SDRangel + +This program is free software; you can redistribute it and/or +modify it under the terms of the GNU General Public License +as published by the Free Software Foundation; either version 2 +of the License, or (at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + +The author can be reached by email at + +warren@wpratt.com + +*/ + +#include "comm.hpp" +#include "fir.hpp" +#include "firopt.hpp" + +namespace WDSP { + +/******************************************************************************************************** +* * +* Standalone Partitioned Overlap-Save Bandpass * +* * +********************************************************************************************************/ + +void FIROPT::plan_firopt (FIROPT *a) +{ + // must call for change in 'nc', 'size', 'out' + int i; + a->nfor = a->nc / a->size; + a->buffidx = 0; + a->idxmask = a->nfor - 1; + a->fftin = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); + a->fftout = new float*[a->nfor]; // (float **) malloc0 (a->nfor * sizeof (float *)); + a->fmask = new float*[a->nfor]; // (float **) malloc0 (a->nfor * sizeof (float *)); + a->maskgen = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); + a->pcfor = new fftwf_plan[a->nfor]; // (fftwf_plan *) malloc0 (a->nfor * sizeof (fftwf_plan)); + a->maskplan = new fftwf_plan[a->nfor]; // (fftwf_plan *) malloc0 (a->nfor * sizeof (fftwf_plan)); + for (i = 0; i < a->nfor; i++) + { + a->fftout[i] = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); + a->fmask[i] = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); + a->pcfor[i] = fftwf_plan_dft_1d(2 * a->size, (fftwf_complex *)a->fftin, (fftwf_complex *)a->fftout[i], FFTW_FORWARD, FFTW_PATIENT); + a->maskplan[i] = fftwf_plan_dft_1d(2 * a->size, (fftwf_complex *)a->maskgen, (fftwf_complex *)a->fmask[i], FFTW_FORWARD, FFTW_PATIENT); + } + a->accum = new float[2 * a->size * 2]; // (float *) malloc0 (2 * a->size * sizeof (complex)); + a->crev = fftwf_plan_dft_1d(2 * a->size, (fftwf_complex *)a->accum, (fftwf_complex *)a->out, FFTW_BACKWARD, FFTW_PATIENT); +} + +void FIROPT::calc_firopt (FIROPT *a) +{ + // call for change in frequency, rate, wintype, gain + // must also call after a call to plan_firopt() + int i; + float* impulse = FIR::fir_bandpass (a->nc, a->f_low, a->f_high, a->samplerate, a->wintype, 1, a->gain); + a->buffidx = 0; + for (i = 0; i < a->nfor; i++) + { + // I right-justified the impulse response => take output from left side of output buff, discard right side + // Be careful about flipping an asymmetrical impulse response. + memcpy (&(a->maskgen[2 * a->size]), &(impulse[2 * a->size * i]), a->size * sizeof(wcomplex)); + fftwf_execute (a->maskplan[i]); + } + delete[] (impulse); +} + +FIROPT* FIROPT::create_firopt (int run, int position, int size, float* in, float* out, + int nc, float f_low, float f_high, int samplerate, int wintype, float gain) +{ + FIROPT *a = new FIROPT; + a->run = run; + a->position = position; + a->size = size; + a->in = in; + a->out = out; + a->nc = nc; + a->f_low = f_low; + a->f_high = f_high; + a->samplerate = samplerate; + a->wintype = wintype; + a->gain = gain; + plan_firopt (a); + calc_firopt (a); + return a; +} + +void FIROPT::deplan_firopt (FIROPT *a) +{ + int i; + fftwf_destroy_plan (a->crev); + delete[] (a->accum); + for (i = 0; i < a->nfor; i++) + { + delete[] (a->fftout[i]); + delete[] (a->fmask[i]); + fftwf_destroy_plan (a->pcfor[i]); + fftwf_destroy_plan (a->maskplan[i]); + } + delete[] (a->maskplan); + delete[] (a->pcfor); + delete[] (a->maskgen); + delete[] (a->fmask); + delete[] (a->fftout); + delete[] (a->fftin); +} + +void FIROPT::destroy_firopt (FIROPT *a) +{ + deplan_firopt (a); + delete (a); +} + +void FIROPT::flush_firopt (FIROPT *a) +{ + int i; + memset (a->fftin, 0, 2 * a->size * sizeof (wcomplex)); + for (i = 0; i < a->nfor; i++) + memset (a->fftout[i], 0, 2 * a->size * sizeof (wcomplex)); + a->buffidx = 0; +} + +void FIROPT::xfiropt (FIROPT *a, int pos) +{ + if (a->run && (a->position == pos)) + { + int i, j, k; + memcpy (&(a->fftin[2 * a->size]), a->in, a->size * sizeof (wcomplex)); + fftwf_execute (a->pcfor[a->buffidx]); + k = a->buffidx; + memset (a->accum, 0, 2 * a->size * sizeof (wcomplex)); + for (j = 0; j < a->nfor; j++) + { + for (i = 0; i < 2 * a->size; i++) + { + a->accum[2 * i + 0] += a->fftout[k][2 * i + 0] * a->fmask[j][2 * i + 0] - a->fftout[k][2 * i + 1] * a->fmask[j][2 * i + 1]; + a->accum[2 * i + 1] += a->fftout[k][2 * i + 0] * a->fmask[j][2 * i + 1] + a->fftout[k][2 * i + 1] * a->fmask[j][2 * i + 0]; + } + k = (k + a->idxmask) & a->idxmask; + } + a->buffidx = (a->buffidx + 1) & a->idxmask; + fftwf_execute (a->crev); + memcpy (a->fftin, &(a->fftin[2 * a->size]), a->size * sizeof(wcomplex)); + } + else if (a->in != a->out) + memcpy (a->out, a->in, a->size * sizeof (wcomplex)); +} + +void FIROPT::setBuffers_firopt (FIROPT *a, float* in, float* out) +{ + a->in = in; + a->out = out; + deplan_firopt (a); + plan_firopt (a); + calc_firopt (a); +} + +void FIROPT::setSamplerate_firopt (FIROPT *a, int rate) +{ + a->samplerate = rate; + calc_firopt (a); +} + +void FIROPT::setSize_firopt (FIROPT *a, int size) +{ + a->size = size; + deplan_firopt (a); + plan_firopt (a); + calc_firopt (a); +} + +void FIROPT::setFreqs_firopt (FIROPT *a, float f_low, float f_high) +{ + a->f_low = f_low; + a->f_high = f_high; + calc_firopt (a); +} + + +} // namespace WDSP diff --git a/wdsp/firopt.hpp b/wdsp/firopt.hpp new file mode 100644 index 000000000..05ea82cfe --- /dev/null +++ b/wdsp/firopt.hpp @@ -0,0 +1,86 @@ +/* firmin.h + +This file is part of a program that implements a Software-Defined Radio. + +Copyright (C) 2016 Warren Pratt, NR0V +Copyright (C) 2024 Edouard Griffiths, F4EXB Adapted to SDRangel + +This program is free software; you can redistribute it and/or +modify it under the terms of the GNU General Public License +as published by the Free Software Foundation; either version 2 +of the License, or (at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + +The author can be reached by email at + +warren@wpratt.com + +*/ + +/******************************************************************************************************** +* * +* Standalone Partitioned Overlap-Save Bandpass * +* * +********************************************************************************************************/ + +#ifndef wdsp_firopt_h +#define wdsp_firopt_h + +#include "fftw3.h" +#include "export.h" + +namespace WDSP { + +class WDSP_API FIROPT +{ + int run; // run control + int position; // position at which to execute + int size; // input/output buffer size, power of two + float* in; // input buffer + float* out; // output buffer, can be same as input + int nc; // number of filter coefficients, power of two, >= size + float f_low; // low cutoff frequency + float f_high; // high cutoff frequency + float samplerate; // sample rate + int wintype; // filter window type + float gain; // filter gain + int nfor; // number of buffers in delay line + float* fftin; // fft input buffer + float** fmask; // frequency domain masks + float** fftout; // fftout delay line + float* accum; // frequency domain accumulator + int buffidx; // fft out buffer index + int idxmask; // mask for index computations + float* maskgen; // input for mask generation FFT + fftwf_plan* pcfor; // array of forward FFT plans + fftwf_plan crev; // reverse fft plan + fftwf_plan* maskplan; // plans for frequency domain masks + + static FIROPT* create_firopt (int run, int position, int size, float* in, float* out, + int nc, float f_low, float f_high, int samplerate, int wintype, float gain); + static void xfiropt (FIROPT *a, int pos); + static void destroy_firopt (FIROPT *a); + static void flush_firopt (FIROPT *a); + static void setBuffers_firopt (FIROPT *a, float* in, float* out); + static void setSamplerate_firopt (FIROPT *a, int rate); + static void setSize_firopt (FIROPT *a, int size); + static void setFreqs_firopt (FIROPT *a, float f_low, float f_high); + +private: + static void plan_firopt (FIROPT *a); + static void calc_firopt (FIROPT *a); + static void deplan_firopt (FIROPT *a); +}; + +} // namespace WDSP + +#endif + diff --git a/wdsp/fmd.cpp b/wdsp/fmd.cpp index 2b4ca47ff..aa2d8435c 100644 --- a/wdsp/fmd.cpp +++ b/wdsp/fmd.cpp @@ -27,7 +27,7 @@ warren@wpratt.com #include "comm.hpp" #include "iir.hpp" -#include "firmin.hpp" +#include "fircore.hpp" #include "fcurve.hpp" #include "fir.hpp" #include "wcpAGC.hpp" diff --git a/wdsp/fmmod.cpp b/wdsp/fmmod.cpp index 2d68cc96e..d8e023ba4 100644 --- a/wdsp/fmmod.cpp +++ b/wdsp/fmmod.cpp @@ -26,7 +26,7 @@ warren@wpratt.com */ #include "comm.hpp" -#include "firmin.hpp" +#include "fircore.hpp" #include "fir.hpp" #include "fmmod.hpp" #include "TXA.hpp" diff --git a/wdsp/fmsq.cpp b/wdsp/fmsq.cpp index 8bc2756e6..2cd146b27 100644 --- a/wdsp/fmsq.cpp +++ b/wdsp/fmsq.cpp @@ -26,7 +26,7 @@ warren@wpratt.com */ #include "comm.hpp" -#include "firmin.hpp" +#include "fircore.hpp" #include "eq.hpp" #include "fmsq.hpp" #include "RXA.hpp" diff --git a/wdsp/icfir.cpp b/wdsp/icfir.cpp index 63f2747ba..5f52153e3 100644 --- a/wdsp/icfir.cpp +++ b/wdsp/icfir.cpp @@ -26,7 +26,7 @@ warren@pratt.one */ #include "comm.hpp" -#include "firmin.hpp" +#include "fircore.hpp" #include "fir.hpp" #include "icfir.hpp" diff --git a/wdsp/nbp.cpp b/wdsp/nbp.cpp index edb55c17c..f8fc5ad17 100644 --- a/wdsp/nbp.cpp +++ b/wdsp/nbp.cpp @@ -27,7 +27,7 @@ warren@wpratt.com #include "comm.hpp" #include "fir.hpp" -#include "firmin.hpp" +#include "fircore.hpp" #include "bpsnba.hpp" #include "nbp.hpp" #include "RXA.hpp" diff --git a/wdsp/patchpanel.cpp b/wdsp/patchpanel.cpp index e5f41ce73..5e685bfa1 100644 --- a/wdsp/patchpanel.cpp +++ b/wdsp/patchpanel.cpp @@ -32,7 +32,17 @@ warren@wpratt.com namespace WDSP { -PANEL* PANEL::create_panel (int run, int size, float* in, float* out, float gain1, float gain2I, float gain2Q, int inselect, int copy) +PANEL* PANEL::create_panel ( + int run, + int size, + float* in, + float* out, + float gain1, + float gain2I, + float gain2Q, + int inselect, + int copy +) { PANEL* a = new PANEL; a->run = run; diff --git a/wdsp/patchpanel.hpp b/wdsp/patchpanel.hpp index 524a8f6f4..34b3ea712 100644 --- a/wdsp/patchpanel.hpp +++ b/wdsp/patchpanel.hpp @@ -48,7 +48,17 @@ public: int inselect; int copy; - static PANEL* create_panel (int run, int size, float* in, float* out, float gain1, float gain2I, float gain2Q, int inselect, int copy); + static PANEL* create_panel ( + int run, + int size, + float* in, + float* out, + float gain1, + float gain2I, + float gain2Q, + int inselect, + int copy + ); static void destroy_panel (PANEL *a); static void flush_panel (PANEL *a); static void xpanel (PANEL *a);