From 65cb3c386e204880bc97022bbdf8e2d73f6a1a20 Mon Sep 17 00:00:00 2001 From: Uladzimir Karpenka Date: Fri, 10 Jul 2026 07:43:07 +0300 Subject: [PATCH 1/3] wfm: add a wideband FM modulator and demodulator Model the pair on fmd.c / fmmod.c, but with the parts that a 75 kHz deviation forces: - the demodulator discriminates with arg(x[n] * conj(x[n-1])) rather than a PLL; an omegaN in the tens of kHz cannot track 75 kHz. - emphasis is a one-pole RC (tau 75 us) on both ends, not fmd's 1/f fc_impulse FIR. A 1/f FIR from f_low = 20 Hz would sit ~+57 dB at 20 Hz, where broadcast FM specifies flat below the corner. - TXA_WFM leaves the shared preemph block off; wfmmod carries its own. - wfmmod clamps bp_fc = deviation + f_high to 0.45 * samplerate, since +/-90 kHz exceeds Nyquist at the rates the narrowband modes use. Scope is mono: no 19 kHz pilot, no 38 kHz stereo subcarrier, no RDS. Both mode enums are appended to so the ABI stays stable for clients. The JNI bindings are deliberately left alone. Verified against a wfmmod -> wfmd loopback: the discriminator is exact, the dc-removal one-pole tracks |H_lp(f)| * ain * sdelta to ratio 1.000, and a 700 + 1900 Hz two-tone comes back with THD+N ~ 3e-5 % at +/-37.5 kHz deviation. Note the demodulator aliases unless samplerate exceeds 2 * deviation * |aud|max, so 192 kHz is the practical floor. Co-Authored-By: Claude Opus 4.8 --- Makefile | 2 + Makefile.android | 2 + RXA.c | 34 ++++- RXA.h | 7 +- TXA.c | 40 +++++- TXA.h | 7 +- comm.h | 2 + wdsp.h | 24 ++++ wfmd.c | 322 +++++++++++++++++++++++++++++++++++++++++++++++ wfmd.h | 103 +++++++++++++++ wfmmod.c | 250 ++++++++++++++++++++++++++++++++++++ wfmmod.h | 86 +++++++++++++ 12 files changed, 871 insertions(+), 8 deletions(-) create mode 100644 wfmd.c create mode 100644 wfmd.h create mode 100644 wfmmod.c create mode 100644 wfmmod.h diff --git a/Makefile b/Makefile index 618486f..f882c96 100644 --- a/Makefile +++ b/Makefile @@ -144,6 +144,8 @@ utilities.c \ varsamp.c \ version.c \ wcpAGC.c \ +wfmd.c \ +wfmmod.c \ wisdom.c \ zetaHat.c diff --git a/Makefile.android b/Makefile.android index a540ea5..79185b0 100644 --- a/Makefile.android +++ b/Makefile.android @@ -94,6 +94,8 @@ utilities.c \ varsamp.c \ version.c \ wcpAGC.c \ +wfmd.c \ +wfmmod.c \ wisdom.c \ zetaHat.c diff --git a/RXA.c b/RXA.c index e264471..ac29df2 100644 --- a/RXA.c +++ b/RXA.c @@ -233,6 +233,23 @@ void create_rxa (int channel) max(2048, ch[channel].dsp_size), // number of coefficients for noise filter 0); // minimum phase flag + // WFM demod + rxa[channel].wfmd.p = create_wfmd ( + 0, // run + ch[channel].dsp_size, // buffer size + rxa[channel].midbuff, // pointer to input buffer + rxa[channel].midbuff, // pointer to output buffer + ch[channel].dsp_rate, // sample rate + 75000.0, // deviation + 20.0, // f_low + 15000.0, // f_high + 0.02, // tau - for dc removal + 1, // run de-emphasis + 75.0e-6, // de-emphasis time constant + 0.5, // audio gain + max(2048, ch[channel].dsp_size), // # coefs for audio cutoff filter + 0); // min phase flag for audio cutoff filter + // snba rxa[channel].snba.p = create_snba ( 0, // run @@ -579,6 +596,7 @@ void destroy_rxa (int channel) destroy_anf (rxa[channel].anf.p); destroy_eqp (rxa[channel].eqp.p); destroy_snba (rxa[channel].snba.p); + destroy_wfmd (rxa[channel].wfmd.p); destroy_fmsq (rxa[channel].fmsq.p); destroy_fmd (rxa[channel].fmd.p); destroy_amd (rxa[channel].amd.p); @@ -614,6 +632,7 @@ void flush_rxa (int channel) flush_amd (rxa[channel].amd.p); flush_fmd (rxa[channel].fmd.p); flush_fmsq (rxa[channel].fmsq.p); + flush_wfmd (rxa[channel].wfmd.p); flush_snba (rxa[channel].snba.p); flush_eqp (rxa[channel].eqp.p); flush_anf (rxa[channel].anf.p); @@ -649,6 +668,7 @@ void xrxa (int channel) xamd (rxa[channel].amd.p); xfmd (rxa[channel].fmd.p); xfmsq (rxa[channel].fmsq.p); + xwfmd (rxa[channel].wfmd.p); xbpsnbain (rxa[channel].bpsnba.p, 1); xbpsnbaout (rxa[channel].bpsnba.p, 1); xsnba (rxa[channel].snba.p); @@ -733,6 +753,7 @@ void setDSPSamplerate_rxa (int channel) setSamplerate_fmd (rxa[channel].fmd.p, ch[channel].dsp_rate); setBuffers_fmsq (rxa[channel].fmsq.p, rxa[channel].midbuff, rxa[channel].midbuff, rxa[channel].fmd.p->audio); setSamplerate_fmsq (rxa[channel].fmsq.p, ch[channel].dsp_rate); + setSamplerate_wfmd (rxa[channel].wfmd.p, ch[channel].dsp_rate); setSamplerate_snba (rxa[channel].snba.p, ch[channel].dsp_rate); setSamplerate_eqp (rxa[channel].eqp.p, ch[channel].dsp_rate); setSamplerate_anf (rxa[channel].anf.p, ch[channel].dsp_rate); @@ -794,6 +815,8 @@ void setDSPBuffsize_rxa (int channel) setSize_fmd (rxa[channel].fmd.p, ch[channel].dsp_size); setBuffers_fmsq (rxa[channel].fmsq.p, rxa[channel].midbuff, rxa[channel].midbuff, rxa[channel].fmd.p->audio); setSize_fmsq (rxa[channel].fmsq.p, ch[channel].dsp_size); + setBuffers_wfmd (rxa[channel].wfmd.p, rxa[channel].midbuff, rxa[channel].midbuff); + setSize_wfmd (rxa[channel].wfmd.p, ch[channel].dsp_size); setBuffers_snba (rxa[channel].snba.p, rxa[channel].midbuff, rxa[channel].midbuff); setSize_snba (rxa[channel].snba.p, ch[channel].dsp_size); setBuffers_eqp (rxa[channel].eqp.p, rxa[channel].midbuff, rxa[channel].midbuff); @@ -857,6 +880,7 @@ void SetRXAMode (int channel, int mode) rxa[channel].mode = mode; rxa[channel].amd.p->run = 0; rxa[channel].fmd.p->run = 0; + rxa[channel].wfmd.p->run = 0; rxa[channel].agc.p->run = 1; switch (mode) { @@ -875,6 +899,10 @@ void SetRXAMode (int channel, int mode) rxa[channel].fmd.p->run = 1; rxa[channel].agc.p->run = 0; break; + case RXA_WFM: + rxa[channel].wfmd.p->run = 1; + rxa[channel].agc.p->run = 0; + break; default: break; @@ -961,13 +989,14 @@ void RXAbpsnbaCheck (int channel, int mode, int notch_run) run_notches = 0; break; case RXA_FM: + case RXA_WFM: f_low = +a->abs_low_freq; f_high = +a->abs_high_freq; run_notches = 0; break; case RXA_DRM: case RXA_SPEC: - + break; } // 'run' and 'position' are examined at run time; no filter changes required. @@ -1010,6 +1039,7 @@ void RXAbpsnbaSet (int channel) a->position = 1; break; case RXA_FM: + case RXA_WFM: a->run = rxa[channel].snba.p->run; a->position = 1; break; @@ -1046,6 +1076,7 @@ void RXASetNC (int channel, int nc) SetRXAFMSQNC (channel, nc); SetRXAFMNCde (channel, nc); SetRXAFMNCaud (channel, nc); + SetRXAWFMNCaud (channel, nc); SetChannelState (channel, oldstate, 0); } @@ -1059,4 +1090,5 @@ void RXASetMP (int channel, int mp) SetRXAFMSQMP (channel, mp); SetRXAFMMPde (channel, mp); SetRXAFMMPaud (channel, mp); + SetRXAWFMMPaud (channel, mp); } diff --git a/RXA.h b/RXA.h index 0f0a378..36c3cb5 100644 --- a/RXA.h +++ b/RXA.h @@ -41,7 +41,8 @@ enum rxaMode RXA_SPEC, RXA_DIGL, RXA_SAM, - RXA_DRM + RXA_DRM, + RXA_WFM }; enum rxaMeterType @@ -121,6 +122,10 @@ struct _rxa FMSQ p; } fmsq; struct + { + WFMD p; + } wfmd; + struct { EQP p; } eqp; diff --git a/TXA.c b/TXA.c index 7ea9c10..6b8d4e6 100644 --- a/TXA.c +++ b/TXA.c @@ -357,7 +357,22 @@ void create_txa (int channel) 1, // run bandpass filter max(2048, ch[channel].dsp_size), // number coefficients for bandpass filter 0); // minimum phase flag - + + txa[channel].wfmmod.p = create_wfmmod ( + 0, // run - OFF by default + ch[channel].dsp_size, // size + txa[channel].midbuff, // pointer to input buffer + txa[channel].midbuff, // pointer to output buffer + ch[channel].dsp_rate, // samplerate + 75000.0, // deviation + 20.0, // low cutoff frequency + 15000.0, // high cutoff frequency + 1, // run pre-emphasis + 75.0e-6, // pre-emphasis time constant + 1, // run bandpass filter + max(2048, ch[channel].dsp_size), // number coefficients for bandpass filter + 0); // minimum phase flag + txa[channel].gen1.p = create_gen ( 0, // run ch[channel].dsp_size, // buffer size @@ -490,6 +505,7 @@ void destroy_txa (int channel) destroy_meter (txa[channel].alcmeter.p); destroy_uslew (txa[channel].uslew.p); destroy_gen (txa[channel].gen1.p); + destroy_wfmmod (txa[channel].wfmmod.p); destroy_fmmod (txa[channel].fmmod.p); destroy_ammod (txa[channel].ammod.p); destroy_wcpagc (txa[channel].alc.p); @@ -544,6 +560,7 @@ void flush_txa (int channel) flush_wcpagc (txa[channel].alc.p); flush_ammod (txa[channel].ammod.p); flush_fmmod (txa[channel].fmmod.p); + flush_wfmmod (txa[channel].wfmmod.p); flush_gen (txa[channel].gen1.p); flush_uslew (txa[channel].uslew.p); flush_meter (txa[channel].alcmeter.p); @@ -580,6 +597,7 @@ void xtxa (int channel) xammod (txa[channel].ammod.p); // AM Modulator xemphp (txa[channel].preemph.p, 1); // FM pre-emphasis (second option) xfmmod (txa[channel].fmmod.p); // FM Modulator + xwfmmod (txa[channel].wfmmod.p); // WFM Modulator (pre-emphasis is internal) xgen (txa[channel].gen1.p); // output signal generator (TUN and Two-tone) xuslew (txa[channel].uslew.p); // up-slew for AM, FM, and gens xmeter (txa[channel].alcmeter.p); // ALC Meter @@ -653,6 +671,7 @@ void setDSPSamplerate_txa (int channel) setSamplerate_wcpagc (txa[channel].alc.p, ch[channel].dsp_rate); setSamplerate_ammod (txa[channel].ammod.p, ch[channel].dsp_rate); setSamplerate_fmmod (txa[channel].fmmod.p, ch[channel].dsp_rate); + setSamplerate_wfmmod (txa[channel].wfmmod.p, ch[channel].dsp_rate); setSamplerate_gen (txa[channel].gen1.p, ch[channel].dsp_rate); setSamplerate_uslew (txa[channel].uslew.p, ch[channel].dsp_rate); setSamplerate_meter (txa[channel].alcmeter.p, ch[channel].dsp_rate); @@ -723,6 +742,8 @@ void setDSPBuffsize_txa (int channel) setSize_ammod (txa[channel].ammod.p, ch[channel].dsp_size); setBuffers_fmmod (txa[channel].fmmod.p, txa[channel].midbuff, txa[channel].midbuff); setSize_fmmod (txa[channel].fmmod.p, ch[channel].dsp_size); + setBuffers_wfmmod (txa[channel].wfmmod.p, txa[channel].midbuff, txa[channel].midbuff); + setSize_wfmmod (txa[channel].wfmmod.p, ch[channel].dsp_size); setBuffers_gen (txa[channel].gen1.p, txa[channel].midbuff, txa[channel].midbuff); setSize_gen (txa[channel].gen1.p, ch[channel].dsp_size); setBuffers_uslew (txa[channel].uslew.p, txa[channel].midbuff, txa[channel].midbuff); @@ -758,6 +779,7 @@ void SetTXAMode (int channel, int mode) txa[channel].mode = mode; txa[channel].ammod.p->run = 0; txa[channel].fmmod.p->run = 0; + txa[channel].wfmmod.p->run = 0; txa[channel].preemph.p->run = 0; switch (mode) { @@ -779,6 +801,10 @@ void SetTXAMode (int channel, int mode) txa[channel].fmmod.p->run = 1; txa[channel].preemph.p->run = 1; break; + case TXA_WFM: + // wfmmod carries its own RC pre-emphasis; the shared emphp stays off + txa[channel].wfmmod.p->run = 1; + break; default: break; @@ -818,10 +844,11 @@ void TXAResCheck (int channel) int TXAUslewCheck (int channel) { - return (txa[channel].ammod.p->run == 1) || - (txa[channel].fmmod.p->run == 1) || - (txa[channel].gen0.p->run == 1) || - (txa[channel].gen1.p->run == 1); + return (txa[channel].ammod.p->run == 1) || + (txa[channel].fmmod.p->run == 1) || + (txa[channel].wfmmod.p->run == 1) || + (txa[channel].gen0.p->run == 1) || + (txa[channel].gen1.p->run == 1); } void TXASetupBPFilters (int channel) @@ -855,6 +882,7 @@ void TXASetupBPFilters (int channel) case TXA_AM: case TXA_SAM: case TXA_FM: + case TXA_WFM: if (txa[channel].compressor.p->run) { CalcBandpassFilter (txa[channel].bp0.p, 0.0, txa[channel].f_high, 2.0); @@ -914,6 +942,7 @@ void TXASetNC (int channel, int nc) SetTXAFMEmphNC (channel, nc); SetTXAEQNC (channel, nc); SetTXAFMNC (channel, nc); + SetTXAWFMNC (channel, nc); SetTXACFIRNC (channel, nc); SetChannelState (channel, oldstate, 0); } @@ -925,6 +954,7 @@ void TXASetMP (int channel, int mp) SetTXAFMEmphMP (channel, mp); SetTXAEQMP (channel, mp); SetTXAFMMP (channel, mp); + SetTXAWFMMP (channel, mp); } PORT diff --git a/TXA.h b/TXA.h index 4e94810..7ecb996 100644 --- a/TXA.h +++ b/TXA.h @@ -43,7 +43,8 @@ enum txaMode TXA_SAM, TXA_DRM, TXA_AM_LSB, - TXA_AM_USB + TXA_AM_USB, + TXA_WFM }; enum txaMeterType @@ -135,6 +136,10 @@ struct _txa FMMOD p; } fmmod; struct + { + WFMMOD p; + } wfmmod; + struct { SIPHON p; } sip1; diff --git a/comm.h b/comm.h index 2477e36..c9fd6a4 100644 --- a/comm.h +++ b/comm.h @@ -77,6 +77,8 @@ warren@wpratt.com #include "fmd.h" #include "fmmod.h" #include "fmsq.h" +#include "wfmd.h" +#include "wfmmod.h" #include "gain.h" #include "gaussian.h" #include "gen.h" diff --git a/wdsp.h b/wdsp.h index f88dc0d..0826681 100644 --- a/wdsp.h +++ b/wdsp.h @@ -517,6 +517,30 @@ extern void SetTXAFMNC (int channel, int nc); extern void SetTXAFMMP (int channel, int mp); extern void SetTXAFMAFFreqs (int channel, double low, double high); +// +// Interfaces from wfmd.c +// + +extern void SetRXAWFMDeviation (int channel, double deviation); +extern void SetRXAWFMNCaud (int channel, int nc); +extern void SetRXAWFMMPaud (int channel, int mp); +extern void SetRXAWFMAFFilter (int channel, double low, double high); +extern void SetRXAWFMDeemphRun (int channel, int run); +extern void SetRXAWFMDeemphTau (int channel, double tau); +extern void SetRXAWFMLimRun (int channel, int run); +extern void SetRXAWFMLimGain (int channel, double gaindB); + +// +// Interfaces from wfmmod.c +// + +extern void SetTXAWFMDeviation (int channel, double deviation); +extern void SetTXAWFMNC (int channel, int nc); +extern void SetTXAWFMMP (int channel, int mp); +extern void SetTXAWFMAFFreqs (int channel, double low, double high); +extern void SetTXAWFMPreEmphRun (int channel, int run); +extern void SetTXAWFMPreEmphTau (int channel, double tau); + // // Interfaces from fmsq.c // diff --git a/wfmd.c b/wfmd.c new file mode 100644 index 0000000..da0f896 --- /dev/null +++ b/wfmd.c @@ -0,0 +1,322 @@ +/* wfmd.c + +This file is part of a program that implements a Software-Defined Radio. + +Copyright (C) 2013, 2023 Warren Pratt, NR0V + +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.h" + +void calc_wfmd (WFMD a) +{ + // discriminator + a->pre_i = 0.0; + a->pre_q = 0.0; + a->again = a->rate / (a->deviation * TWOPI); + // dc removal + a->mtau = exp(-1.0 / (a->rate * a->tau)); + a->onem_mtau = 1.0 - a->mtau; + a->fmdc = 0.0; + // de-emphasis + if (a->tau_de > 0.0) a->mde = exp(-1.0 / (a->rate * a->tau_de)); + else a->mde = 0.0; + a->onem_mde = 1.0 - a->mde; + a->deemph_z = 0.0; + // detector limiter + a->plim = create_wcpagc ( + 1, // run - always ON + 5, // mode + 1, // 0 for max(I,Q), 1 for envelope + a->out, // input buff pointer + a->out, // output buff pointer + a->size, // io_buffsize + (int)a->rate, // sample rate + 0.001, // tau_attack + 0.008, // tau_decay + 4, // n_tau + a->lim_gain, // max_gain (sets threshold, initial value) + 1.0, // var_gain / slope + 1.0, // fixed_gain + 1.0, // max_input + 0.9, // out_targ + 0.250, // tau_fast_backaverage + 0.004, // tau_fast_decay + 4.0, // pop_ratio + 0, // hang_enable + 0.500, // tau_hang_backmult + 0.500, // hangtime + 2.000, // hang_thresh + 0.100); // tau_hang_decay +} + +void decalc_wfmd (WFMD a) +{ + destroy_wcpagc(a->plim); +} + +WFMD create_wfmd (int run, int size, double* in, double* out, int rate, double deviation, double f_low, double f_high, + double tau, int deemph_run, double tau_de, double afgain, int nc_aud, int mp_aud) +{ + WFMD a = (WFMD) malloc0 (sizeof (wfmd)); + double* impulse; + a->run = run; + a->size = size; + a->in = in; + a->out = out; + a->rate = (double)rate; + a->deviation = deviation; + a->f_low = f_low; + a->f_high = f_high; + a->tau = tau; + a->deemph_run = deemph_run; + a->tau_de = tau_de; + a->afgain = afgain; + a->nc_aud = nc_aud; + a->mp_aud = mp_aud; + a->lim_run = 0; + a->lim_pre_gain = 0.4; + a->lim_gain = 2.5; + calc_wfmd (a); + a->audio = (double *) malloc0 (a->size * sizeof (complex)); + // audio filter + impulse = fir_bandpass(a->nc_aud, 0.8 * a->f_low, 1.1 * a->f_high, a->rate, 0, 1, a->afgain / (2.0 * a->size)); + a->paud = create_fircore (a->size, a->audio, a->out, a->nc_aud, a->mp_aud, impulse); + _aligned_free (impulse); + return a; +} + +void destroy_wfmd (WFMD a) +{ + destroy_fircore (a->paud); + _aligned_free (a->audio); + decalc_wfmd (a); + _aligned_free (a); +} + +void flush_wfmd (WFMD a) +{ + memset (a->audio, 0, a->size * sizeof (complex)); + flush_fircore (a->paud); + a->pre_i = 0.0; + a->pre_q = 0.0; + a->fmdc = 0.0; + a->deemph_z = 0.0; + flush_wcpagc (a->plim); +} + +void xwfmd (WFMD a) +{ + if (a->run) + { + int i; + double si, sq, cr, ci, det, aud; + for (i = 0; i < a->size; i++) + { + // quadrature discriminator: det = arg (x[n] * conj (x[n-1])) + si = a->in[2 * i + 0]; + sq = a->in[2 * i + 1]; + cr = + si * a->pre_i + sq * a->pre_q; + ci = - si * a->pre_q + sq * a->pre_i; + a->pre_i = si; + a->pre_q = sq; + det = atan2 (ci, cr); + // dc removal, gain, & demod output + a->fmdc = a->mtau * a->fmdc + a->onem_mtau * det; + aud = a->again * (det - a->fmdc); + // de-emphasis + if (a->deemph_run) + { + a->deemph_z = a->mde * a->deemph_z + a->onem_mde * aud; + aud = a->deemph_z; + } + a->audio[2 * i + 0] = aud; + a->audio[2 * i + 1] = aud; + } + // audio filter + xfircore (a->paud); + if (a->lim_run) + { + for (i = 0; i < 2 * a->size; i++) + a->out[i] *= a->lim_pre_gain; + xwcpagc (a->plim); + } + } + else if (a->in != a->out) + memcpy (a->out, a->in, a->size * sizeof (complex)); +} + +void setBuffers_wfmd (WFMD a, double* in, double* out) +{ + decalc_wfmd (a); + a->in = in; + a->out = out; + calc_wfmd (a); + setBuffers_fircore (a->paud, a->audio, a->out); + setBuffers_wcpagc (a->plim, a->out, a->out); +} + +void setSamplerate_wfmd (WFMD a, int rate) +{ + double* impulse; + decalc_wfmd (a); + a->rate = rate; + calc_wfmd (a); + // audio filter + impulse = fir_bandpass(a->nc_aud, 0.8 * a->f_low, 1.1 * a->f_high, a->rate, 0, 1, a->afgain / (2.0 * a->size)); + setImpulse_fircore (a->paud, impulse, 1); + _aligned_free (impulse); + setSamplerate_wcpagc (a->plim, (int)a->rate); +} + +void setSize_wfmd (WFMD a, int size) +{ + double* impulse; + decalc_wfmd (a); + _aligned_free (a->audio); + a->size = size; + calc_wfmd (a); + a->audio = (double *) malloc0 (a->size * sizeof (complex)); + // audio filter + destroy_fircore (a->paud); + impulse = fir_bandpass(a->nc_aud, 0.8 * a->f_low, 1.1 * a->f_high, a->rate, 0, 1, a->afgain / (2.0 * a->size)); + a->paud = create_fircore (a->size, a->audio, a->out, a->nc_aud, a->mp_aud, impulse); + _aligned_free (impulse); + setSize_wcpagc (a->plim, a->size); +} + +/******************************************************************************************************** +* * +* RXA Properties * +* * +********************************************************************************************************/ + +PORT +void SetRXAWFMDeviation (int channel, double deviation) +{ + WFMD a; + EnterCriticalSection (&ch[channel].csDSP); + a = rxa[channel].wfmd.p; + a->deviation = deviation; + a->again = a->rate / (a->deviation * TWOPI); + LeaveCriticalSection (&ch[channel].csDSP); +} + +PORT +void SetRXAWFMNCaud (int channel, int nc) +{ + WFMD a; + double* impulse; + EnterCriticalSection (&ch[channel].csDSP); + a = rxa[channel].wfmd.p; + if (a->nc_aud != nc) + { + a->nc_aud = nc; + impulse = fir_bandpass(a->nc_aud, 0.8 * a->f_low, 1.1 * a->f_high, a->rate, 0, 1, a->afgain / (2.0 * a->size)); + setNc_fircore (a->paud, a->nc_aud, impulse); + _aligned_free (impulse); + } + LeaveCriticalSection (&ch[channel].csDSP); +} + +PORT +void SetRXAWFMMPaud (int channel, int mp) +{ + WFMD a; + a = rxa[channel].wfmd.p; + if (a->mp_aud != mp) + { + a->mp_aud = mp; + setMp_fircore (a->paud, a->mp_aud); + } +} + +PORT +void SetRXAWFMAFFilter (int channel, double low, double high) +{ + WFMD a = rxa[channel].wfmd.p; + double* impulse; + EnterCriticalSection (&ch[channel].csDSP); + if (a->f_low != low || a->f_high != high) + { + a->f_low = low; + a->f_high = high; + impulse = fir_bandpass (a->nc_aud, 0.8 * a->f_low, 1.1 * a->f_high, a->rate, 0, 1, a->afgain / (2.0 * a->size)); + setImpulse_fircore (a->paud, impulse, 1); + _aligned_free (impulse); + } + LeaveCriticalSection (&ch[channel].csDSP); +} + +PORT +void SetRXAWFMDeemphRun (int channel, int run) +{ + WFMD a = rxa[channel].wfmd.p; + EnterCriticalSection (&ch[channel].csDSP); + if (a->deemph_run != run) + { + a->deemph_run = run; + a->deemph_z = 0.0; + } + LeaveCriticalSection (&ch[channel].csDSP); +} + +PORT +void SetRXAWFMDeemphTau (int channel, double tau) +{ + WFMD a = rxa[channel].wfmd.p; + EnterCriticalSection (&ch[channel].csDSP); + if (a->tau_de != tau && tau > 0.0) + { + a->tau_de = tau; + a->mde = exp(-1.0 / (a->rate * a->tau_de)); + a->onem_mde = 1.0 - a->mde; + a->deemph_z = 0.0; + } + LeaveCriticalSection (&ch[channel].csDSP); +} + +PORT +void SetRXAWFMLimRun (int channel, int run) +{ + WFMD a = rxa[channel].wfmd.p; + EnterCriticalSection (&ch[channel].csDSP); + if (a->lim_run != run) + { + a->lim_run = run; + } + LeaveCriticalSection (&ch[channel].csDSP); +} + +PORT +void SetRXAWFMLimGain (int channel, double gaindB) +{ + double gain = pow(10.0, gaindB / 20.0); + WFMD a = rxa[channel].wfmd.p; + EnterCriticalSection (&ch[channel].csDSP); + if (a->lim_gain != gain) + { + decalc_wfmd (a); + a->lim_gain = gain; + calc_wfmd (a); + } + LeaveCriticalSection (&ch[channel].csDSP); +} diff --git a/wfmd.h b/wfmd.h new file mode 100644 index 0000000..78343e5 --- /dev/null +++ b/wfmd.h @@ -0,0 +1,103 @@ +/* wfmd.h + +This file is part of a program that implements a Software-Defined Radio. + +Copyright (C) 2013 Warren Pratt, NR0V + +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 + +*/ + +#ifndef _wfmd_h +#define _wfmd_h +#include "firmin.h" +#include "wcpAGC.h" +typedef struct _wfmd +{ + int run; + int size; + double* in; + double* out; + double rate; + double f_low; // audio low cutoff + double f_high; // audio high cutoff + // quadrature discriminator + double deviation; // peak deviation, Hz + double again; // discriminator output gain + double pre_i; // previous sample, I + double pre_q; // previous sample, Q + // for dc removal + double tau; + double mtau; + double onem_mtau; + double fmdc; + // de-emphasis, single-pole RC + int deemph_run; + double tau_de; // 75.0e-6 (Americas) or 50.0e-6 (elsewhere) + double mde; + double onem_mde; + double deemph_z; + // for audio filter + double* audio; + FIRCORE paud; + int nc_aud; + int mp_aud; + double afgain; + // detector limiter + WCPAGC plim; + int lim_run; + double lim_gain; + double lim_pre_gain; +} wfmd, *WFMD; + +extern WFMD create_wfmd ( int run, int size, double* in, double* out, int rate, double deviation, + double f_low, double f_high, double tau, int deemph_run, double tau_de, double afgain, + int nc_aud, int mp_aud); + +extern void destroy_wfmd (WFMD a); + +extern void flush_wfmd (WFMD a); + +extern void xwfmd (WFMD a); + +extern void setBuffers_wfmd (WFMD a, double* in, double* out); + +extern void setSamplerate_wfmd (WFMD a, int rate); + +extern void setSize_wfmd (WFMD a, int size); + +// RXA Properties + +extern __declspec (dllexport) void SetRXAWFMDeviation (int channel, double deviation); + +extern __declspec (dllexport) void SetRXAWFMNCaud (int channel, int nc); + +extern __declspec (dllexport) void SetRXAWFMMPaud (int channel, int mp); + +extern __declspec (dllexport) void SetRXAWFMAFFilter (int channel, double low, double high); + +extern __declspec (dllexport) void SetRXAWFMDeemphRun (int channel, int run); + +extern __declspec (dllexport) void SetRXAWFMDeemphTau (int channel, double tau); + +extern __declspec (dllexport) void SetRXAWFMLimRun (int channel, int run); + +extern __declspec (dllexport) void SetRXAWFMLimGain (int channel, double gaindB); + +#endif diff --git a/wfmmod.c b/wfmmod.c new file mode 100644 index 0000000..0ceb3bc --- /dev/null +++ b/wfmmod.c @@ -0,0 +1,250 @@ +/* wfmmod.c + +This file is part of a program that implements a Software-Defined Radio. + +Copyright (C) 2013, 2016, 2023 Warren Pratt, NR0V + +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.h" + +// the modulated spectrum spans +/-(deviation + f_high); at the sample rates used for +// narrowband modes that exceeds Nyquist, so the bandpass degenerates to a passthrough. +static double bpfc_wfmmod (double samplerate, double deviation, double f_high) +{ + double fc = deviation + f_high; + double max_fc = 0.45 * samplerate; + if (fc > max_fc) fc = max_fc; + return fc; +} + +void calc_wfmmod (WFMMOD a) +{ + // pre-emphasis + if (a->tau_pre > 0.0) a->pmult = exp(-1.0 / (a->samplerate * a->tau_pre)); + else a->pmult = 0.0; + a->pnorm = 1.0 / (1.0 - a->pmult); + a->pre_z = 0.0; + // mod + a->sphase = 0.0; + a->sdelta = TWOPI * a->deviation / a->samplerate; + // bandpass + a->bp_fc = bpfc_wfmmod (a->samplerate, a->deviation, a->f_high); +} + +WFMMOD create_wfmmod (int run, int size, double* in, double* out, int rate, double dev, double f_low, double f_high, + int pre_run, double tau_pre, int bp_run, int nc, int mp) +{ + WFMMOD a = (WFMMOD) malloc0 (sizeof (wfmmod)); + double* impulse; + a->run = run; + a->size = size; + a->in = in; + a->out = out; + a->samplerate = (double)rate; + a->deviation = dev; + a->f_low = f_low; + a->f_high = f_high; + a->pre_run = pre_run; + a->tau_pre = tau_pre; + a->bp_run = bp_run; + a->nc = nc; + a->mp = mp; + calc_wfmmod (a); + impulse = fir_bandpass(a->nc, -a->bp_fc, +a->bp_fc, a->samplerate, 0, 1, 1.0 / (2 * a->size)); + a->p = create_fircore (a->size, a->out, a->out, a->nc, a->mp, impulse); + _aligned_free (impulse); + return a; +} + +void destroy_wfmmod (WFMMOD a) +{ + destroy_fircore (a->p); + _aligned_free (a); +} + +void flush_wfmmod (WFMMOD a) +{ + a->pre_z = 0.0; + a->sphase = 0.0; + flush_fircore (a->p); +} + +void xwfmmod (WFMMOD a) +{ + int i; + double aud, dp; + if (a->run) + { + for (i = 0; i < a->size; i++) + { + aud = a->in[2 * i + 0]; + if (a->pre_run) + { + dp = a->pnorm * (aud - a->pmult * a->pre_z); + a->pre_z = aud; + aud = dp; + } + dp = aud * a->sdelta; + a->sphase += dp; + // at the wide deviation, |dp| exceeds TWOPI once samplerate < deviation, + // so one subtraction is not enough to bring sphase back into range + while (a->sphase >= TWOPI) a->sphase -= TWOPI; + while (a->sphase < 0.0 ) a->sphase += TWOPI; + a->out[2 * i + 0] = 0.7071 * cos (a->sphase); + a->out[2 * i + 1] = 0.7071 * sin (a->sphase); + } + if (a->bp_run) + xfircore (a->p); + } + else if (a->in != a->out) + memcpy (a->out, a->in, a->size * sizeof (complex)); +} + +void setBuffers_wfmmod (WFMMOD a, double* in, double* out) +{ + a->in = in; + a->out = out; + calc_wfmmod (a); + setBuffers_fircore (a->p, a->out, a->out); +} + +void setSamplerate_wfmmod (WFMMOD a, int rate) +{ + double* impulse; + a->samplerate = rate; + calc_wfmmod (a); + impulse = fir_bandpass(a->nc, -a->bp_fc, +a->bp_fc, a->samplerate, 0, 1, 1.0 / (2 * a->size)); + setImpulse_fircore (a->p, impulse, 1); + _aligned_free (impulse); +} + +void setSize_wfmmod (WFMMOD a, int size) +{ + double* impulse; + a->size = size; + calc_wfmmod (a); + setSize_fircore (a->p, a->size); + impulse = fir_bandpass(a->nc, -a->bp_fc, +a->bp_fc, a->samplerate, 0, 1, 1.0 / (2 * a->size)); + setImpulse_fircore (a->p, impulse, 1); + _aligned_free (impulse); +} + +/******************************************************************************************************** +* * +* TXA Properties * +* * +********************************************************************************************************/ + +PORT +void SetTXAWFMDeviation (int channel, double deviation) +{ + WFMMOD a = txa[channel].wfmmod.p; + double bp_fc = bpfc_wfmmod (a->samplerate, deviation, a->f_high); + double* impulse = fir_bandpass (a->nc, -bp_fc, +bp_fc, a->samplerate, 0, 1, 1.0 / (2 * a->size)); + setImpulse_fircore (a->p, impulse, 0); + _aligned_free (impulse); + EnterCriticalSection (&ch[channel].csDSP); + a->deviation = deviation; + // mod + a->sphase = 0.0; + a->sdelta = TWOPI * a->deviation / a->samplerate; + // bandpass + a->bp_fc = bp_fc; + setUpdate_fircore (a->p); + LeaveCriticalSection (&ch[channel].csDSP); +} + +PORT +void SetTXAWFMNC (int channel, int nc) +{ + WFMMOD a; + double* impulse; + EnterCriticalSection (&ch[channel].csDSP); + a = txa[channel].wfmmod.p; + if (a->nc != nc) + { + a->nc = nc; + impulse = fir_bandpass (a->nc, -a->bp_fc, +a->bp_fc, a->samplerate, 0, 1, 1.0 / (2 * a->size)); + setNc_fircore (a->p, a->nc, impulse); + _aligned_free (impulse); + } + LeaveCriticalSection (&ch[channel].csDSP); +} + +PORT +void SetTXAWFMMP (int channel, int mp) +{ + WFMMOD a; + a = txa[channel].wfmmod.p; + if (a->mp != mp) + { + a->mp = mp; + setMp_fircore (a->p, a->mp); + } +} + +PORT +void SetTXAWFMAFFreqs (int channel, double low, double high) +{ + WFMMOD a; + double* impulse; + EnterCriticalSection (&ch[channel].csDSP); + a = txa[channel].wfmmod.p; + if (a->f_low != low || a->f_high != high) + { + a->f_low = low; + a->f_high = high; + a->bp_fc = bpfc_wfmmod (a->samplerate, a->deviation, a->f_high); + impulse = fir_bandpass (a->nc, -a->bp_fc, +a->bp_fc, a->samplerate, 0, 1, 1.0 / (2 * a->size)); + setImpulse_fircore (a->p, impulse, 1); + _aligned_free (impulse); + } + LeaveCriticalSection (&ch[channel].csDSP); +} + +PORT +void SetTXAWFMPreEmphRun (int channel, int run) +{ + WFMMOD a = txa[channel].wfmmod.p; + EnterCriticalSection (&ch[channel].csDSP); + if (a->pre_run != run) + { + a->pre_run = run; + a->pre_z = 0.0; + } + LeaveCriticalSection (&ch[channel].csDSP); +} + +PORT +void SetTXAWFMPreEmphTau (int channel, double tau) +{ + WFMMOD a = txa[channel].wfmmod.p; + EnterCriticalSection (&ch[channel].csDSP); + if (a->tau_pre != tau && tau > 0.0) + { + a->tau_pre = tau; + a->pmult = exp(-1.0 / (a->samplerate * a->tau_pre)); + a->pnorm = 1.0 / (1.0 - a->pmult); + a->pre_z = 0.0; + } + LeaveCriticalSection (&ch[channel].csDSP); +} diff --git a/wfmmod.h b/wfmmod.h new file mode 100644 index 0000000..68a4b05 --- /dev/null +++ b/wfmmod.h @@ -0,0 +1,86 @@ +/* wfmmod.h + +This file is part of a program that implements a Software-Defined Radio. + +Copyright (C) 2013, 2016, 2023 Warren Pratt, NR0V + +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 + +*/ + +#ifndef _wfmmod_h +#define _wfmmod_h +#include "firmin.h" +typedef struct _wfmmod +{ + int run; + int size; + double* in; + double* out; + double samplerate; + double deviation; + double f_low; + double f_high; + // pre-emphasis, single-pole RC; inverse of the receiver's de-emphasis + int pre_run; + double tau_pre; // 75.0e-6 (Americas) or 50.0e-6 (elsewhere) + double pmult; + double pnorm; + double pre_z; + // mod + double sphase; + double sdelta; + // bandpass + int bp_run; + double bp_fc; + int nc; + int mp; + FIRCORE p; +}wfmmod, *WFMMOD; + +extern WFMMOD create_wfmmod (int run, int size, double* in, double* out, int rate, double dev, double f_low, double f_high, + int pre_run, double tau_pre, int bp_run, int nc, int mp); + +extern void destroy_wfmmod (WFMMOD a); + +extern void flush_wfmmod (WFMMOD a); + +extern void xwfmmod (WFMMOD a); + +extern void setBuffers_wfmmod (WFMMOD a, double* in, double* out); + +extern void setSamplerate_wfmmod (WFMMOD a, int rate); + +extern void setSize_wfmmod (WFMMOD a, int size); + +// TXA Properties + +extern __declspec (dllexport) void SetTXAWFMDeviation (int channel, double deviation); + +extern __declspec (dllexport) void SetTXAWFMNC (int channel, int nc); + +extern __declspec (dllexport) void SetTXAWFMMP (int channel, int mp); + +extern __declspec (dllexport) void SetTXAWFMAFFreqs (int channel, double low, double high); + +extern __declspec (dllexport) void SetTXAWFMPreEmphRun (int channel, int run); + +extern __declspec (dllexport) void SetTXAWFMPreEmphTau (int channel, double tau); + +#endif From c1964817ddd97398d4fe5df64d42306b547faf86 Mon Sep 17 00:00:00 2001 From: Uladzimir Karpenka Date: Wed, 15 Jul 2026 11:03:03 +0300 Subject: [PATCH 2/3] wfm: compile wfmd/wfmmod in the Windows build Makefile.windows arrived from main, which predates the wfm module, so its SOURCES list omitted wfmd.c/wfmmod.c even though RXA.c/TXA.c reference their symbols -- the MinGW link failed on the missing objects. Add both to SOURCES, matching Makefile and Makefile.android. Co-Authored-By: Claude Opus 4.8 --- Makefile.windows | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Makefile.windows b/Makefile.windows index 0730dec..b5ee1cb 100644 --- a/Makefile.windows +++ b/Makefile.windows @@ -86,7 +86,7 @@ SOURCES = \ matchedCW.c meter.c meterlog10.c nbp.c nob.c nobII.c osctrl.c \ patchpanel.c resample.c rmatch.c rnnr.c RXA.c sbnr.c sender.c \ shift.c siphon.c slew.c snb.c ssql.c syncbuffs.c TXA.c \ - utilities.c varsamp.c version.c wcpAGC.c wisdom.c zetaHat.c + utilities.c varsamp.c version.c wcpAGC.c wfmd.c wfmmod.c wisdom.c zetaHat.c OBJS := $(addprefix $(OBJDIR)/, $(SOURCES:.c=.o)) From 3ad40a00dab447e16229b25f46511f4386e6f496 Mon Sep 17 00:00:00 2001 From: Uladzimir Karpenka Date: Fri, 17 Jul 2026 23:00:41 +0300 Subject: [PATCH 3/3] calcc: port eu2av robustness patches for PureSignal calibration Three fixes from Thetis-Enhanced (Yurij eu2av), measured on Orion MK2 / Anvelina PRO3 hardware: - drop overrange samples (env_TX*hw_scale > 1.0) before the cubic xbuilder fit - they distort the fit and produce a wrong rx_scale - optional median-ratio + MAD outlier rejection before the fit, controlled by new SetPSOutlierSigma (0 = off) - fallback rx_scale estimate from the top amplitude intervals when the xbuilder fit fails or is rejected by rxscheck Co-Authored-By: Claude Fable 5 --- calcc.c | 173 +++++++++++++++++++++++++++++++++++++++++++++++++++++++- calcc.h | 1 + wdsp.h | 1 + 3 files changed, 173 insertions(+), 2 deletions(-) diff --git a/calcc.c b/calcc.c index 6e2bf8d..87afce7 100644 --- a/calcc.c +++ b/calcc.c @@ -138,6 +138,7 @@ CALCC create_calcc (int channel, int runcal, int size, int rate, int ints, int s a->stbl = stbl; a->npsamps = npsamps; a->alpha = alpha; + a->outlier_sigma = 0.0; a->info = (int *) malloc0 (16 * sizeof (int)); a->binfo = (int *) malloc0 (16 * sizeof (int)); @@ -321,6 +322,125 @@ void rxscheck (int rints, double* tvec, double* coef, int* info) if (out < 0.00) *info |= 0x0020; } +// Yurij_eu2av: fallback rx_scale estimator. It averages the top few +// amplitude intervals (ignoring overrange samples) and linearly extrapolates +// to full TX scale (env_TX = 1/hw_scale). Used only if the cubic xbuilder +// fit fails or is rejected by rxscheck. +static int estimate_rx_scale_from_top_intervals(CALCC a, double* rx_scale_out) +{ + const int n_top = 4; + double sx[4], sy[4], sw[4]; + int valid = 0; + int b, j; + for (b = a->ints - 1; b >= 0 && valid < n_top; b--) + { + int base = b * a->spi; + double sum_x = 0.0, sum_y = 0.0; + int n = 0; + for (j = 0; j < a->spi; j++) + { + int k = base + j; + double nx = a->env_TX[k] * a->hw_scale; + if (nx > 1.0 || nx < 0.0) continue; + if (a->env_TX[k] < 1.0e-30 || a->env_RX[k] < 1.0e-30) continue; + sum_x += a->env_TX[k]; + sum_y += a->env_RX[k]; + n++; + } + if (n == 0) continue; + sx[valid] = sum_x / (double)n; + sy[valid] = sum_y / (double)n; + sw[valid] = (double)n; + valid++; + } + if (valid < 2) return -1; + { + double s_w = 0.0, s_x = 0.0, s_y = 0.0, s_xx = 0.0, s_xy = 0.0; + double det, aa, bb, target_x, y_at_target; + int i; + for (i = 0; i < valid; i++) + { + double w = sw[i]; + s_w += w; + s_x += w * sx[i]; + s_y += w * sy[i]; + s_xx += w * sx[i] * sx[i]; + s_xy += w * sx[i] * sy[i]; + } + det = s_w * s_xx - s_x * s_x; + if (fabs(det) < 1e-30) return -1; + bb = (s_w * s_xy - s_x * s_y) / det; + aa = (s_y - bb * s_x) / s_w; + target_x = 1.0 / a->hw_scale; + y_at_target = aa + bb * target_x; + if (y_at_target <= 1e-15) return -1; + *rx_scale_out = 1.0 / y_at_target; + } + return 0; +} + +// Yurij_eu2av: robust outlier rejection for the cubic-spline xbuilder. +// Fits rx = k*tx through the origin via median ratio, then rejects points +// whose residual exceeds sigma * MAD. +static int cmp_double(const void* a, const void* b) +{ + double da = *(const double*)a; + double db = *(const double*)b; + if (da < db) return -1; + if (da > db) return 1; + return 0; +} + +static double median_double(double* v, int n) +{ + if (n <= 0) return 0.0; + if (n % 2 == 1) + return v[n / 2]; + else + return 0.5 * (v[n / 2 - 1] + v[n / 2]); +} + +static int reject_outliers(double* tx, double* rx, int n, double sigma) +{ + const int min_points = 32; + int i, keep = 0; + double* ratios; + double* absres; + double med_ratio, med_absres, thr; + + if (n < min_points || sigma <= 0.0) return n; + + ratios = (double*)malloc0(n * sizeof(double)); + for (i = 0; i < n; i++) + ratios[i] = (tx[i] > 1.0e-30) ? rx[i] / tx[i] : 0.0; + qsort(ratios, n, sizeof(double), cmp_double); + med_ratio = median_double(ratios, n); + _aligned_free(ratios); + + if (fabs(med_ratio) < 1.0e-30) return n; + + absres = (double*)malloc0(n * sizeof(double)); + for (i = 0; i < n; i++) + absres[i] = fabs(rx[i] - med_ratio * tx[i]); + qsort(absres, n, sizeof(double), cmp_double); + med_absres = median_double(absres, n); + _aligned_free(absres); + + if (med_absres < 1.0e-30) return n; + + thr = sigma * med_absres; + for (i = 0; i < n; i++) + { + if (fabs(rx[i] - med_ratio * tx[i]) <= thr) + { + tx[keep] = tx[i]; + rx[keep] = rx[i]; + keep++; + } + } + return (keep >= min_points) ? keep : n; +} + void calc (CALCC a) { int i; @@ -336,21 +456,60 @@ void calc (CALCC a) double tvec[3]; double txrxcoefs[4 * 2]; double rx_scale; + int xb_ok = 0; + double* tx_filt; + double* rx_filt; + int n_filt = 0; if (a->ints < 16) rints = 1; else rints = 2; ix = rints - 1; for (i = 0; i <= rints; i++) tvec[i] = (double)i / (double)rints / a->hw_scale; dx = tvec[rints] - tvec[rints - 1]; - xbuilder(a->ccbld, a->nsamps, a->env_TX, a->env_RX, rints, tvec, &(a->binfo[0]), txrxcoefs, a->ptol); + + // Yurij_eu2av: build a filtered dataset with overrange samples removed + // before running xbuilder. Overrange env_TX*hw_scale > 1.0 can distort + // the cubic fit and produce an incorrect rx_scale. + tx_filt = (double*)malloc0(a->nsamps * sizeof(double)); + rx_filt = (double*)malloc0(a->nsamps * sizeof(double)); + for (i = 0; i < a->nsamps; i++) + { + double nx = a->env_TX[i] * a->hw_scale; + if (nx > 1.0 || nx < 0.0) continue; + if (a->env_TX[i] < 1.0e-30 || a->env_RX[i] < 1.0e-30) continue; + tx_filt[n_filt] = a->env_TX[i]; + rx_filt[n_filt] = a->env_RX[i]; + n_filt++; + } + + // Yurij_eu2av: optional outlier rejection before cubic-spline fit. + if (a->outlier_sigma > 0.0) + n_filt = reject_outliers(tx_filt, rx_filt, n_filt, a->outlier_sigma); + + xbuilder(a->ccbld, n_filt, tx_filt, rx_filt, rints, tvec, &(a->binfo[0]), txrxcoefs, a->ptol); rxscheck (rints, tvec, txrxcoefs, &a->binfo[7]); if ((a->binfo[0] == 0) && (a->binfo[7] == 0)) + { rx_scale = 1.0 / (txrxcoefs[4 * ix + 0] + dx * (txrxcoefs[4 * ix + 1] + dx * (txrxcoefs[4 * ix + 2] + dx * txrxcoefs[4 * ix + 3]))); - else + xb_ok = 1; + } + else if (estimate_rx_scale_from_top_intervals(a, &rx_scale) == 0) + { + // Yurij_eu2av: xbuilder failed, but the bucket-average fallback + // gave a usable rx_scale. Keep binfo[0] bit 0 set for diagnostics. + a->binfo[0] |= 0x0001; + xb_ok = 1; + } + + _aligned_free(tx_filt); + _aligned_free(rx_filt); + + if (!xb_ok) { a->scOK = 0; goto cleanup; } + if (a->stbl && _InterlockedAnd (&a->ctrl.running, 1)) a->rx_scale = a->alpha * a->rx_scale + (1.0 - a->alpha) * rx_scale; else @@ -1046,6 +1205,16 @@ void SetPSPtol (int channel, double ptol) LeaveCriticalSection (&txa[channel].calcc.cs_update); } +PORT +void SetPSOutlierSigma (int channel, double sigma) +{ + // Yurij_eu2av: 0.0 disables the pre-xbuilder outlier filter. + if (sigma < 0.0) sigma = 0.0; + EnterCriticalSection (&txa[channel].calcc.cs_update); + txa[channel].calcc.p->outlier_sigma = sigma; + LeaveCriticalSection (&txa[channel].calcc.cs_update); +} + PORT void GetPSDisp (int channel, double* x, double* ym, double* yc, double* ys, double* cm, double* cc, double* cs) { diff --git a/calcc.h b/calcc.h index 7495430..06d0b49 100644 --- a/calcc.h +++ b/calcc.h @@ -48,6 +48,7 @@ typedef struct _calcc double hw_scale; double rx_scale; double alpha; + double outlier_sigma; int tsamps; double* env_TX; diff --git a/wdsp.h b/wdsp.h index 0826681..1d56c32 100644 --- a/wdsp.h +++ b/wdsp.h @@ -277,6 +277,7 @@ extern void SetPSHWPeak (int channel, double peak); extern void GetPSHWPeak (int channel, double* peak); extern void GetPSMaxTX (int channel, double* maxtx); extern void SetPSPtol (int channel, double ptol); +extern void SetPSOutlierSigma (int channel, double sigma); extern void GetPSDisp (int channel, double* x, double* ym, double* yc, double* ys, double* cm, double* cc, double* cs); extern void SetPSFeedbackRate (int channel, int rate); extern void SetPSPinMode (int channel, int pin);