Merge branch 'wfm' into cpu-optimization

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 <noreply@anthropic.com>
2026-07-18 10:51:50 +03:00 · 2026-07-18 10:47:44 +03:00 · 2026-07-15 11:03:03 +03:00 · 2026-07-15 10:50:34 +03:00 · 2026-07-15 10:50:15 +03:00 · 2026-07-10 07:43:07 +03:00
29 changed files with 1750 additions and 243 deletions
@@ -144,6 +144,8 @@ utilities.c \
 varsamp.c \
 version.c \
 wcpAGC.c \
+wfmd.c \
+wfmmod.c \
 wisdom.c \
 zetaHat.c

@@ -119,6 +119,8 @@ utilities.c \
 varsamp.c \
 version.c \
 wcpAGC.c \
+wfmd.c \
+wfmmod.c \
 wisdom.c \
 zetaHat.c

@@ -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))

@@ -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,6 +989,7 @@ 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;
@@ -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);
 }
@@ -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;
@@ -358,6 +358,21 @@ void create_txa (int channel)
 		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;
@@ -820,6 +846,7 @@ int TXAUslewCheck (int channel)
 {
 	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);
 }
@@ -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
@@ -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;
@@ -1087,7 +1087,15 @@ void __cdecl sendbuf(void *arg)
 					LeaveCriticalSection(&(a->BufferControlSection[a->ss][a->LO]));
 				}
 			}
-		Sleep(1);
+		//
+		// Block until a Spectrum*() call announces new samples. This used to be
+		// Sleep(1), i.e. 1000 wakeups per second spent re-reading the same flags
+		// -- 0.67% of a core even with no data arriving at all.
+		//
+		// Whoever sets end_dispatcher also signals the semaphore, so the wait
+		// below always has a way out.
+		//
+		WaitForSingleObject(a->Sem_BuffReady, INFINITE);
 	}
 	InterlockedBitTestAndReset(&a->dispatcher, 0);
 	_endthread();
@@ -1202,6 +1210,9 @@ void SetAnalyzer (	int disp,			// display identifier

 	EnterCriticalSection(&a->SetAnalyzerSection);
 	a->end_dispatcher = 1;
+	// wake the dispatcher out of its blocking wait so it can observe the flag;
+	// it does not take SetAnalyzerSection, so holding it here is safe
+	ReleaseSemaphore(a->Sem_BuffReady, 1, 0);
 	while (InterlockedAnd(&a->dispatcher, 1))
 		Sleep(1);
 	a->stop = 1;
@@ -1347,6 +1358,7 @@ void XCreateAnalyzer(	int disp,
 			a->hSnapEvent[i][j] = CreateEvent(NULL, FALSE, FALSE, TEXT("snap"));
 			a->snap[i][j] = 0;
 		}
+	a->Sem_BuffReady = CreateSemaphore(0, 0, 1000, 0);
 	InitializeCriticalSectionAndSpinCount(&a->ResampleSection, 0);
 	InitializeCriticalSectionAndSpinCount(&a->SetAnalyzerSection, 0);
 	InitializeCriticalSectionAndSpinCount(&a->StitchSection, 0);
@@ -1434,6 +1446,7 @@ void DestroyAnalyzer(int disp)
 	int i, j;

 	a->end_dispatcher = 1;
+	ReleaseSemaphore(a->Sem_BuffReady, 1, 0);
 	while (InterlockedAnd(&a->dispatcher, 1))
 		Sleep(1);

@@ -1497,6 +1510,7 @@ void DestroyAnalyzer(int disp)
 	for (i = 0; i < a->max_stitch; i++)
 		for (j = 0; j < a->max_num_fft; j++)
 			CloseHandle(a->hSnapEvent[i][j]);
+	CloseHandle(a->Sem_BuffReady);

 	_aligned_free ((void *) a->pnum_threads);

@@ -1633,6 +1647,7 @@ void CloseBuffer(int disp, int ss, int LO)
 	if((a->IQin_index[ss][LO] += a->buff_size) >= a->bsize) //REQUIRES buff_size IS A SUB-MULTIPLE OF SIZE OF INPUT SAMPLE BUFFS!
 		a->IQin_index[ss][LO] = 0;

+	ReleaseSemaphore(a->Sem_BuffReady, 1, 0);	// new samples: let the dispatcher run
 	if (!InterlockedAnd(&a->dispatcher, 1))
 	{
 		InterlockedBitTestAndSet (&a->dispatcher, 0);
@@ -1672,6 +1687,7 @@ void Spectrum(int disp, int ss, int LO, dINREAL* pI, dINREAL* pQ)
 	if((a->IQin_index[ss][LO] += a->buff_size) >= a->bsize) //REQUIRES buff_size IS A SUB-MULTIPLE OF SIZE OF INPUT SAMPLE BUFFS!
 		a->IQin_index[ss][LO] = 0;

+	ReleaseSemaphore(a->Sem_BuffReady, 1, 0);	// new samples: let the dispatcher run
 	if (!InterlockedAnd(&a->dispatcher, 1))
 	{
 		InterlockedBitTestAndSet(&a->dispatcher, 0);
@@ -1717,6 +1733,7 @@ void Spectrum2(int run, int disp, int ss, int LO, dINREAL* pbuff)
 		if((a->IQin_index[ss][LO] += a->buff_size) >= a->bsize) //REQUIRES buff_size IS A SUB-MULTIPLE OF SIZE OF INPUT SAMPLE BUFFS!
 			a->IQin_index[ss][LO] = 0;

+		ReleaseSemaphore(a->Sem_BuffReady, 1, 0);	// new samples: let the dispatcher run
 		if (!InterlockedAnd(&a->dispatcher, 1))
 		{
 			InterlockedBitTestAndSet(&a->dispatcher, 0);
@@ -1763,6 +1780,7 @@ void Spectrum0(int run, int disp, int ss, int LO, double* pbuff)
 		if((a->IQin_index[ss][LO] += a->buff_size) >= a->bsize) //REQUIRES buff_size IS A SUB-MULTIPLE OF SIZE OF INPUT SAMPLE BUFFS!
 			a->IQin_index[ss][LO] = 0;

+		ReleaseSemaphore(a->Sem_BuffReady, 1, 0);	// new samples: let the dispatcher run
 		if (!InterlockedAnd(&a->dispatcher, 1))
 		{
 			InterlockedBitTestAndSet(&a->dispatcher, 0);
@@ -121,6 +121,9 @@ typedef struct _dp
 	HANDLE hSnapEvent[dMAX_STITCH][dMAX_NUM_FFT];			// mutex handles; mutexes will be used to signal a snap is complete
 	double *snap_buff[dMAX_STITCH][dMAX_NUM_FFT];			// pointers to buffers for the snap

+	HANDLE Sem_BuffReady;									// signalled when input samples arrive, so the
+															//   dispatcher can block instead of polling
+
 	CRITICAL_SECTION PB_ControlsSection[dMAX_PIXOUTS];
 	CRITICAL_SECTION SetAnalyzerSection;
 	CRITICAL_SECTION BufferControlSection[dMAX_STITCH][dMAX_NUM_FFT];
@@ -26,6 +26,45 @@ warren@wpratt.com

 #include "comm.h"

+/* Filter output and tap-window energy over a unit-stride run of the delay line.
+
+   The delay line is indexed (in_idx + j + delay) & mask, which wraps at most
+   once across the tap window; xanf() splits the window at the wrap so both
+   halves are contiguous here.  Four independent accumulator pairs keep the FMAs
+   off a single dependency chain and let the vectorizer in -- a 'y += w[j]*x[j]'
+   reduction cannot be reassociated without -ffast-math, which this library must
+   not enable (it relies on IEEE semantics for 0/0 = NaN and x/0 = Inf). */
+static inline void anf_dot (const double* WDSP_RESTRICT w,
+	const double* WDSP_RESTRICT x, int n, double* py, double* psigma)
+{
+	double y0 = 0.0, y1 = 0.0, y2 = 0.0, y3 = 0.0;
+	double s0 = 0.0, s1 = 0.0, s2 = 0.0, s3 = 0.0;
+	int j = 0;
+	for (; j <= n - 4; j += 4)
+	{
+		y0 += w[j + 0] * x[j + 0];	s0 += x[j + 0] * x[j + 0];
+		y1 += w[j + 1] * x[j + 1];	s1 += x[j + 1] * x[j + 1];
+		y2 += w[j + 2] * x[j + 2];	s2 += x[j + 2] * x[j + 2];
+		y3 += w[j + 3] * x[j + 3];	s3 += x[j + 3] * x[j + 3];
+	}
+	for (; j < n; j++)
+	{
+		y0 += w[j] * x[j];
+		s0 += x[j] * x[j];
+	}
+	*py		+= (y0 + y1) + (y2 + y3);
+	*psigma += (s0 + s1) + (s2 + s3);
+}
+
+/* Leaky-LMS tap update over the same unit-stride run. */
+static inline void anf_update (double* WDSP_RESTRICT w,
+	const double* WDSP_RESTRICT x, int n, double c0, double c1)
+{
+	int j;
+	for (j = 0; j < n; j++)
+		w[j] = c0 * w[j] + c1 * x[j];
+}
+
 ANF create_anf	(
 				int run,
 				int position,
@@ -81,53 +120,81 @@ void destroy_anf (ANF a)

 void xanf(ANF a, int position)
 {
-	int i, j, idx;
+	int i;
 	double c0, c1;
 	double y, error, sigma, inv_sigp;
 	double nel, nev;
 	if (a->run && (a->position == position))
 	{
-		for (i = 0; i < a->buff_size; i++)
-		{
-			a->d[a->in_idx] = a->in_buff[2 * i + 0];
+		const int n_taps = a->n_taps;
+		const int dline_size = a->dline_size;
+		const int mask = a->mask;
+		const int delay = a->delay;
+		const int buff_size = a->buff_size;
+		const double two_mu = a->two_mu;
+		const double gamma = a->gamma;
+		const double den_mult = a->den_mult;
+		const double lincr = a->lincr;
+		const double ldecr = a->ldecr;
+		const double lidx_min = a->lidx_min;
+		const double lidx_max = a->lidx_max;
+		/* in_buff and out_buff are the same buffer in RXA, so neither may be
+		   marked restrict; d and w are private to the struct. */
+		const double* in_buff = a->in_buff;
+		double* out_buff = a->out_buff;
+		double* WDSP_RESTRICT d = a->d;
+		double* WDSP_RESTRICT w = a->w;
+		int in_idx = a->in_idx;
+		double lidx = a->lidx;
+		double ngamma = a->ngamma;

-			y = 0;
-			sigma = 0;
-
-			for (j = 0; j < a->n_taps; j++)
+		for (i = 0; i < buff_size; i++)
 		{
-				idx = (a->in_idx + j + a->delay) & a->mask;
-				y += a->w[j] * a->d[idx];
-				sigma += a->d[idx] * a->d[idx];
-			}
+			double dsamp;
+			int base, n1;
+
+			dsamp = in_buff[2 * i + 0];
+			d[in_idx] = dsamp;
+
+			base = (in_idx + delay) & mask;
+			if ((n1 = dline_size - base) > n_taps) n1 = n_taps;
+
+			y = 0.0;
+			sigma = 0.0;
+			anf_dot (w, d + base, n1, &y, &sigma);
+			if (n1 < n_taps)
+				anf_dot (w + n1, d, n_taps - n1, &y, &sigma);
+
 			inv_sigp = 1.0 / (sigma + 1e-10);
-			error = a->d[a->in_idx] - y;
+			error = dsamp - y;

-			a->out_buff[2 * i + 0] = error;
-			a->out_buff[2 * i + 1] = 0.0;
+			out_buff[2 * i + 0] = error;
+			out_buff[2 * i + 1] = 0.0;

-			if((nel = error * (1.0 - a->two_mu * sigma * inv_sigp)) < 0.0) nel = -nel;
-			if((nev = a->d[a->in_idx] - (1.0 - a->two_mu * a->ngamma) * y - a->two_mu * error * sigma * inv_sigp) < 0.0) nev = -nev;
+			if((nel = error * (1.0 - two_mu * sigma * inv_sigp)) < 0.0) nel = -nel;
+			if((nev = dsamp - (1.0 - two_mu * ngamma) * y - two_mu * error * sigma * inv_sigp) < 0.0) nev = -nev;
 			if (nev < nel)
 			{
-				if ((a->lidx += a->lincr) > a->lidx_max) a->lidx = a->lidx_max;
+				if ((lidx += lincr) > lidx_max) lidx = lidx_max;
 			}
 			else
 			{
-				if ((a->lidx -= a->ldecr) < a->lidx_min) a->lidx = a->lidx_min;
+				if ((lidx -= ldecr) < lidx_min) lidx = lidx_min;
 			}
-			a->ngamma = a->gamma * (a->lidx * a->lidx) * (a->lidx * a->lidx) * a->den_mult;
+			ngamma = gamma * (lidx * lidx) * (lidx * lidx) * den_mult;

-			c0 = 1.0 - a->two_mu * a->ngamma;
-			c1 = a->two_mu * error * inv_sigp;
+			c0 = 1.0 - two_mu * ngamma;
+			c1 = two_mu * error * inv_sigp;

-			for (j = 0; j < a->n_taps; j++)
-			{
-				idx = (a->in_idx + j + a->delay) & a->mask;
-				a->w[j] = c0 * a->w[j] + c1 * a->d[idx];
-			}
-			a->in_idx = (a->in_idx + a->mask) & a->mask;
+			anf_update (w, d + base, n1, c0, c1);
+			if (n1 < n_taps)
+				anf_update (w + n1, d, n_taps - n1, c0, c1);
+
+			in_idx = (in_idx + mask) & mask;
 		}
+		a->in_idx = in_idx;
+		a->lidx = lidx;
+		a->ngamma = ngamma;
 	}
 	else if (a->in_buff != a->out_buff)
 		memcpy (a->out_buff, a->in_buff, a->buff_size * sizeof (complex));
@@ -26,6 +26,45 @@ warren@wpratt.com

 #include "comm.h"

+/* Filter output and tap-window energy over a unit-stride run of the delay line.
+
+   The delay line is indexed (in_idx + j + delay) & mask, which wraps at most
+   once across the tap window; xanr() splits the window at the wrap so both
+   halves are contiguous here.  Four independent accumulator pairs keep the FMAs
+   off a single dependency chain and let the vectorizer in -- an 'y += w[j]*x[j]'
+   reduction cannot be reassociated without -ffast-math, which this library must
+   not enable (it relies on IEEE semantics for 0/0 = NaN and x/0 = Inf). */
+static inline void anr_dot (const double* WDSP_RESTRICT w,
+	const double* WDSP_RESTRICT x, int n, double* py, double* psigma)
+{
+	double y0 = 0.0, y1 = 0.0, y2 = 0.0, y3 = 0.0;
+	double s0 = 0.0, s1 = 0.0, s2 = 0.0, s3 = 0.0;
+	int j = 0;
+	for (; j <= n - 4; j += 4)
+	{
+		y0 += w[j + 0] * x[j + 0];	s0 += x[j + 0] * x[j + 0];
+		y1 += w[j + 1] * x[j + 1];	s1 += x[j + 1] * x[j + 1];
+		y2 += w[j + 2] * x[j + 2];	s2 += x[j + 2] * x[j + 2];
+		y3 += w[j + 3] * x[j + 3];	s3 += x[j + 3] * x[j + 3];
+	}
+	for (; j < n; j++)
+	{
+		y0 += w[j] * x[j];
+		s0 += x[j] * x[j];
+	}
+	*py		+= (y0 + y1) + (y2 + y3);
+	*psigma += (s0 + s1) + (s2 + s3);
+}
+
+/* Leaky-LMS tap update over the same unit-stride run. */
+static inline void anr_update (double* WDSP_RESTRICT w,
+	const double* WDSP_RESTRICT x, int n, double c0, double c1)
+{
+	int j;
+	for (j = 0; j < n; j++)
+		w[j] = c0 * w[j] + c1 * x[j];
+}
+
 ANR create_anr	(
 				int run,
 				int position,
@@ -81,53 +120,81 @@ void destroy_anr (ANR a)

 void xanr (ANR a, int position)
 {
-	int i, j, idx;
+	int i;
 	double c0, c1;
 	double y, error, sigma, inv_sigp;
 	double nel, nev;
 	if (a->run && (a->position == position))
 	{
-		for (i = 0; i < a->buff_size; i++)
-		{
-			a->d[a->in_idx] = a->in_buff[2 * i + 0];
+		const int n_taps = a->n_taps;
+		const int dline_size = a->dline_size;
+		const int mask = a->mask;
+		const int delay = a->delay;
+		const int buff_size = a->buff_size;
+		const double two_mu = a->two_mu;
+		const double gamma = a->gamma;
+		const double den_mult = a->den_mult;
+		const double lincr = a->lincr;
+		const double ldecr = a->ldecr;
+		const double lidx_min = a->lidx_min;
+		const double lidx_max = a->lidx_max;
+		/* in_buff and out_buff are the same buffer in RXA, so neither may be
+		   marked restrict; d and w are private to the struct. */
+		const double* in_buff = a->in_buff;
+		double* out_buff = a->out_buff;
+		double* WDSP_RESTRICT d = a->d;
+		double* WDSP_RESTRICT w = a->w;
+		int in_idx = a->in_idx;
+		double lidx = a->lidx;
+		double ngamma = a->ngamma;

-			y = 0;
-			sigma = 0;
-
-			for (j = 0; j < a->n_taps; j++)
+		for (i = 0; i < buff_size; i++)
 		{
-				idx = (a->in_idx + j + a->delay) & a->mask;
-				y += a->w[j] * a->d[idx];
-				sigma += a->d[idx] * a->d[idx];
-			}
+			double dsamp;
+			int base, n1;
+
+			dsamp = in_buff[2 * i + 0];
+			d[in_idx] = dsamp;
+
+			base = (in_idx + delay) & mask;
+			if ((n1 = dline_size - base) > n_taps) n1 = n_taps;
+
+			y = 0.0;
+			sigma = 0.0;
+			anr_dot (w, d + base, n1, &y, &sigma);
+			if (n1 < n_taps)
+				anr_dot (w + n1, d, n_taps - n1, &y, &sigma);
+
 			inv_sigp = 1.0 / (sigma + 1e-10);
-			error = a->d[a->in_idx] - y;
+			error = dsamp - y;

-			a->out_buff[2 * i + 0] = y;
-			a->out_buff[2 * i + 1] = 0.0;
+			out_buff[2 * i + 0] = y;
+			out_buff[2 * i + 1] = 0.0;

-			if((nel = error * (1.0 - a->two_mu * sigma * inv_sigp)) < 0.0) nel = -nel;
-			if((nev = a->d[a->in_idx] - (1.0 - a->two_mu * a->ngamma) * y - a->two_mu * error * sigma * inv_sigp) < 0.0) nev = -nev;
+			if((nel = error * (1.0 - two_mu * sigma * inv_sigp)) < 0.0) nel = -nel;
+			if((nev = dsamp - (1.0 - two_mu * ngamma) * y - two_mu * error * sigma * inv_sigp) < 0.0) nev = -nev;
 			if (nev < nel)
 			{
-				if ((a->lidx += a->lincr) > a->lidx_max) a->lidx = a->lidx_max;
+				if ((lidx += lincr) > lidx_max) lidx = lidx_max;
 			}
 			else
 			{
-				if ((a->lidx -= a->ldecr) < a->lidx_min) a->lidx = a->lidx_min;
+				if ((lidx -= ldecr) < lidx_min) lidx = lidx_min;
 			}
-			a->ngamma = a->gamma * (a->lidx * a->lidx) * (a->lidx * a->lidx) * a->den_mult;
+			ngamma = gamma * (lidx * lidx) * (lidx * lidx) * den_mult;

-			c0 = 1.0 - a->two_mu * a->ngamma;
-			c1 = a->two_mu * error * inv_sigp;
+			c0 = 1.0 - two_mu * ngamma;
+			c1 = two_mu * error * inv_sigp;

-			for (j = 0; j < a->n_taps; j++)
-			{
-				idx = (a->in_idx + j + a->delay) & a->mask;
-				a->w[j] = c0 * a->w[j] + c1 * a->d[idx];
-			}
-			a->in_idx = (a->in_idx + a->mask) & a->mask;
+			anr_update (w, d + base, n1, c0, c1);
+			if (n1 < n_taps)
+				anr_update (w + n1, d, n_taps - n1, c0, c1);
+
+			in_idx = (in_idx + mask) & mask;
 		}
+		a->in_idx = in_idx;
+		a->lidx = lidx;
+		a->ngamma = ngamma;
 	}
 	else if (a->in_buff != a->out_buff)
 		memcpy (a->out_buff, a->in_buff, a->buff_size * sizeof (complex));
@@ -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)
 {
@@ -48,6 +48,7 @@ typedef struct _calcc
 	double hw_scale;
 	double rx_scale;
 	double alpha;
+	double outlier_sigma;

 	int tsamps;
 	double* env_TX;
@@ -336,18 +336,36 @@ void xcfcomp (CFCOMP a, int pos)
 	if (a->run && pos == a->position)
 	{
 		int i, j, k, sbuff, sbegin;
-		for (i = 0; i < 2 * a->bsize; i += 2)
+		/* Each ring index below steps by one and, since iasize >= fsize and
+		   oasize >= incr always hold, wraps at most once per loop. The '% size'
+		   per step was therefore an integer division for nothing -- about 5100
+		   of them per call at fsize = 2048. Walk contiguous runs instead. */
+		const int iasize = a->iasize;
+		const int oasize = a->oasize;
+		const int fsize = a->fsize;
+		const int incr = a->incr;
+		const int bsize = a->bsize;
+		const int ovrlp = a->ovrlp;
+		const double pregain = a->pregain;
+		const double postgain = a->postgain;
+
+		for (i = 0, j = a->iainidx; i < 2 * bsize; i += 2)
 		{
-			a->inaccum[a->iainidx] = a->in[i];
-			a->iainidx = (a->iainidx + 1) % a->iasize;
+			a->inaccum[j] = a->in[i];
+			if (++j == iasize) j = 0;
 		}
-		a->nsamps += a->bsize;
-		while (a->nsamps >= a->fsize)
+		a->iainidx = j;
+		a->nsamps += bsize;
+		while (a->nsamps >= fsize)
 		{
-			for (i = 0, j = a->iaoutidx; i < a->fsize; i++, j = (j + 1) % a->iasize)
-				a->forfftin[i] = a->pregain * a->window[i] * a->inaccum[j];
-			a->iaoutidx = (a->iaoutidx + a->incr) % a->iasize;
-			a->nsamps -= a->incr;
+			int n1 = iasize - a->iaoutidx;
+			if (n1 > fsize) n1 = fsize;
+			for (i = 0; i < n1; i++)
+				a->forfftin[i] = pregain * a->window[i] * a->inaccum[a->iaoutidx + i];
+			for (; i < fsize; i++)
+				a->forfftin[i] = pregain * a->window[i] * a->inaccum[i - n1];
+			if ((a->iaoutidx += incr) >= iasize) a->iaoutidx -= iasize;
+			a->nsamps -= incr;
 			fftw_execute (a->Rfor);
 			calc_mask(a);
 			for (i = 0; i < a->msize; i++)
@@ -356,30 +374,41 @@ void xcfcomp (CFCOMP a, int pos)
 				a->revfftin[2 * i + 1] = a->mask[i] * a->forfftout[2 * i + 1];
 			}
 			fftw_execute (a->Rrev);
-			for (i = 0; i < a->fsize; i++)
-				a->save[a->saveidx][i] = a->postgain * a->window[i] * a->revfftout[i];
-			for (i = a->ovrlp; i > 0; i--)
+			for (i = 0; i < fsize; i++)
+				a->save[a->saveidx][i] = postgain * a->window[i] * a->revfftout[i];
+			for (i = ovrlp; i > 0; i--)
 			{
-				sbuff = (a->saveidx + i) % a->ovrlp;
-				sbegin = a->incr * (a->ovrlp - i);
-				for (j = sbegin, k = a->oainidx; j < a->incr + sbegin; j++, k = (k + 1) % a->oasize)
+				const double* WDSP_RESTRICT sv;
+				double* WDSP_RESTRICT oa = a->outaccum;
+				int m1;
+				sbuff = (a->saveidx + i) % ovrlp;
+				sbegin = incr * (ovrlp - i);
+				sv = a->save[sbuff] + sbegin;
+				m1 = oasize - a->oainidx;
+				if (m1 > incr) m1 = incr;
+				k = a->oainidx;
+				if (i == ovrlp)
 				{
-					if ( i == a->ovrlp)
-						a->outaccum[k]	= a->save[sbuff][j];
+					for (j = 0; j < m1; j++)		oa[k + j] = sv[j];
+					for (; j < incr; j++)			oa[j - m1] = sv[j];
+				}
 				else
-						a->outaccum[k] += a->save[sbuff][j];
-				}
-			}
-			a->saveidx = (a->saveidx + 1) % a->ovrlp;
-			a->oainidx = (a->oainidx + a->incr) % a->oasize;
-		}
-		for (i = 0; i < a->bsize; i++)
 				{
-			a->out[2 * i + 0] = a->outaccum[a->oaoutidx];
-			a->out[2 * i + 1] = 0.0;
-			a->oaoutidx = (a->oaoutidx + 1) % a->oasize;
+					for (j = 0; j < m1; j++)		oa[k + j] += sv[j];
+					for (; j < incr; j++)			oa[j - m1] += sv[j];
 				}
 			}
+			if (++a->saveidx == ovrlp) a->saveidx = 0;
+			if ((a->oainidx += incr) >= oasize) a->oainidx -= oasize;
+		}
+		for (i = 0, k = a->oaoutidx; i < bsize; i++)
+		{
+			a->out[2 * i + 0] = a->outaccum[k];
+			a->out[2 * i + 1] = 0.0;
+			if (++k == oasize) k = 0;
+		}
+		a->oaoutidx = k;
+	}
 	else if (a->out != a->in)
 		memcpy (a->out, a->in, a->bsize * sizeof (complex));
 }
@@ -75,12 +75,15 @@ warren@wpratt.com
 #include "sbnr.h" // NR3 + NR4 support
 #include "emph.h"
 #include "eq.h"
+#include "fastmath.h"
 #include "fcurve.h"
 #include "fir.h"
 #include "firmin.h"
 #include "fmd.h"
 #include "fmmod.h"
 #include "fmsq.h"
+#include "wfmd.h"
+#include "wfmmod.h"
 #include "gain.h"
 #include "gaussian.h"
 #include "gen.h"
@@ -150,6 +153,14 @@ warren@wpratt.com
 #define PI								3.1415926535897932
 #define TWOPI							6.2831853071795864

+// Non-aliasing qualifier for DSP buffers. Spelled __restrict rather than
+// restrict because the JNI translation unit is compiled as -std=gnu89.
+#if defined(__GNUC__) || defined(__clang__) || defined(_MSC_VER)
+#define WDSP_RESTRICT					__restrict
+#else
+#define WDSP_RESTRICT
+#endif
+
 // miscellaneous
 typedef double complex[2];
 #define PORT							__declspec( dllexport )
@@ -557,6 +557,7 @@ void calc_emnr(EMNR a)
 	a->ae.psi		 = 20.0;
 	a->ae.t2 = 0.20;
 	a->ae.nmask = (double *)malloc0(a->ae.msize * sizeof(double));
+	a->ae.csum	= (double *)malloc0((a->ae.msize + 1) * sizeof(double));
 	//
 	// post2
 	a->post2.run = 0;
@@ -580,6 +581,7 @@ void decalc_emnr(EMNR a)
 		_aligned_free(a->post2.noise_frame);
 		_aligned_free(a->post2.w);
 		// ae
+	_aligned_free(a->ae.csum);
 	_aligned_free(a->ae.nmask);
 	// npl
 	_aligned_free(a->npl.D);
@@ -868,26 +870,26 @@ void aepf(EMNR a)
 	else
 		N = 1 + 2 * (int)(0.5 + a->ae.psi * (1.0 - zetaT / a->ae.zetaThresh));
 	n = N / 2;
-	for (k = 0; k < n; k++)
+	/* Each of the three spans below averages mask[] over a window that is
+	   symmetric about k and clipped at the array ends.  Taking them straight
+	   from a prefix sum makes each output one subtraction rather than a walk of
+	   up to N = 2*psi + 1 taps, so the pass is O(msize) instead of O(msize*N). */
 	{
-		a->ae.nmask[k] = 0.0;
-		for (m = 0; m <= 2 * k; m++)
-			a->ae.nmask[k] += a->mask[m];
-		a->ae.nmask[k] /= (double)(2 * k + 1);
-	}
-	for (k = n; k < (a->ae.msize - n); k++)
-	{
-		a->ae.nmask[k] = 0.0;
-		for (m = k - n; m <= (k + n); m++)
-			a->ae.nmask[k] += a->mask[m];
-		a->ae.nmask[k] /= (double)N;
-	}
-	for (k = a->ae.msize - n; k < a->ae.msize; k++)
-	{
-		a->ae.nmask[k] = 0.0;
-		for (m = (a->ae.msize - 1); m >= (-a->ae.msize + 2 * k + 1); m--)
-			a->ae.nmask[k] += a->mask[m];
-		a->ae.nmask[k] /= (double)(2 * (a->ae.msize - k) - 1);
+		const int msize = a->ae.msize;
+		const double* WDSP_RESTRICT mask = a->mask;
+		double* WDSP_RESTRICT nmask = a->ae.nmask;
+		double* WDSP_RESTRICT csum = a->ae.csum;
+
+		csum[0] = 0.0;
+		for (k = 0; k < msize; k++)
+			csum[k + 1] = csum[k] + mask[k];
+
+		for (k = 0; k < n; k++)					// window [0, 2k]
+			nmask[k] = (csum[2 * k + 1] - csum[0]) / (double)(2 * k + 1);
+		for (k = n; k < (msize - n); k++)		// window [k-n, k+n]
+			nmask[k] = (csum[k + n + 1] - csum[k - n]) / (double)N;
+		for (k = msize - n; k < msize; k++)		// window [2k+1-msize, msize-1]
+			nmask[k] = (csum[msize] - csum[2 * k + 1 - msize]) / (double)(2 * (msize - k) - 1);
 	}
 	memcpy (a->mask, a->ae.nmask, a->ae.msize * sizeof (double));
 	if (a->g.gain_method == 3 && zetaT < a->ae.t2)
@@ -998,52 +1000,54 @@ void SetRXAEMNRpost2Rate(int channel, double tc)
 *								End Post-Processing Functions											*
 ********************************************************************************************************/

-double getKey(double* type, double gamma, double xi)
+/* Locate v on the table's 0.25 dB grid: n1/n2 bracket it, d is the fraction.
+
+   v is compared against the table's dB limits first, so on the interpolating
+   path the log argument is positive and normal and wdsp_log10 applies. */
+static inline void keyIndex (double v, int* n1, int* n2, double* d)
 {
-	int ngamma1, ngamma2, nxi1, nxi2;
-	double tg, tx, dg, dx;
 	const double dmin = 0.001;
 	const double dmax = 1000.0;
-	if (gamma <= dmin)
+	if (v <= dmin)
 	{
-		ngamma1 = ngamma2 = 0;
-		tg = 0.0;
+		*n1 = *n2 = 0;
+		*d = 0.0;
 	}
-	else if (gamma >= dmax)
+	else if (v >= dmax)
 	{
-		ngamma1 = ngamma2 = 240;
-		tg = 60.0;
+		*n1 = *n2 = 240;
+		*d = 0.0;
 	}
 	else
 	{
-		tg = 10.0 * log10(gamma / dmin);
-		ngamma1 = (int)(4.0 * tg);
-		ngamma2 = ngamma1 + 1;
+		double f = 40.0 * wdsp_log10 (v / dmin);		// 4 * (10 * log10)
+		int i = (int)f;
+		/* clamp so n2 cannot address the next row of the 241x241 table */
+		if (i > 239) i = 239;
+		*n1 = i;
+		*n2 = i + 1;
+		*d = f - (double)i;
 	}
-	if (xi <= dmin)
-	{
-		nxi1 = nxi2 = 0;
-		tx = 0.0;
-	}
-	else if (xi >= dmax)
-	{
-		nxi1 = nxi2 = 240;
-		tx = 60.0;
-	}
-	else
-	{
-		tx = 10.0 * log10(xi / dmin);
-		nxi1 = (int)(4.0 * tx);
-		nxi2 = nxi1 + 1;
-	}
-	dg = (tg - 0.25 * ngamma1) / 0.25;
-	dx = (tx - 0.25 * nxi1) / 0.25;
+}
+
+static inline double keyLerp (const double* type, int ngamma1, int ngamma2, double dg,
+	int nxi1, int nxi2, double dx)
+{
 	return (1.0 - dg)  * (1.0 - dx) * type[241 * nxi1 + ngamma1]
 		+  (1.0 - dg)  *		dx	* type[241 * nxi2 + ngamma1]
 		+		  dg   * (1.0 - dx) * type[241 * nxi1 + ngamma2]
 		+		  dg   *		dx	* type[241 * nxi2 + ngamma2];
 }

+double getKey(double* type, double gamma, double xi)
+{
+	int ngamma1, ngamma2, nxi1, nxi2;
+	double dg, dx;
+	keyIndex (gamma, &ngamma1, &ngamma2, &dg);
+	keyIndex (xi,	 &nxi1,	   &nxi2,	 &dx);
+	return keyLerp (type, ngamma1, ngamma2, dg, nxi1, nxi2, dx);
+}
+
 int getZeta( EMNR a, double gamma, double eps, double* zeta)
 {
 	int index, i_gamma, i_xi;
@@ -1133,13 +1137,20 @@ void calc_gain (EMNR a)
 	case 2:
 		{
 			double gamma, eps_hat, eps_p;
+			int ngamma1, ngamma2, nxi1, nxi2, npi1, npi2;
+			double dg, dx, dp;
 			for (k = 0; k < a->g.msize; k++)
 			{
 				gamma = min(a->g.lambda_y[k] / a->g.lambda_d[k], a->g.gamma_max);
 				eps_hat = a->g.alpha * a->g.prev_mask[k] * a->g.prev_mask[k] * a->g.prev_gamma[k]
 					+ (1.0 - a->g.alpha) * max(gamma - 1.0, a->g.eps_floor);
 				eps_p = eps_hat / (1.0 - a->g.q);
-				a->g.mask[k] = getKey(a->g.GG, gamma, eps_hat) * getKey(a->g.GGS, gamma, eps_p);
+				/* both lookups share gamma, so locate it once */
+				keyIndex (gamma,   &ngamma1, &ngamma2, &dg);
+				keyIndex (eps_hat, &nxi1,	 &nxi2,	   &dx);
+				keyIndex (eps_p,   &npi1,	 &npi2,	   &dp);
+				a->g.mask[k] = keyLerp (a->g.GG,  ngamma1, ngamma2, dg, nxi1, nxi2, dx)
+							 * keyLerp (a->g.GGS, ngamma1, ngamma2, dg, npi1, npi2, dp);
 				a->g.prev_gamma[k] = gamma;
 				a->g.prev_mask[k] = a->g.mask[k];
 			}
@@ -1203,18 +1214,34 @@ void xemnr (EMNR a, int pos)
 	{
 		int i, j, k, sbuff, sbegin;
 		double g1;
-		for (i = 0; i < 2 * a->bsize; i += 2)
+		/* The ring indices below advance by one per iteration and wrap at most
+		   once per loop, so a '% size' each step is an integer division for
+		   nothing (iasize = 3584 here, not a power of two). Walk contiguous runs
+		   and wrap between them instead. */
+		const int iasize = a->iasize;
+		const int oasize = a->oasize;
+		const int fsize = a->fsize;
+		const int incr = a->incr;
+		const int bsize = a->bsize;
+		const int ovrlp = a->ovrlp;
+
+		for (i = 0, j = a->iainidx; i < 2 * bsize; i += 2)
 		{
-			a->inaccum[a->iainidx] = a->in[i];
-			a->iainidx = (a->iainidx + 1) % a->iasize;
+			a->inaccum[j] = a->in[i];
+			if (++j == iasize) j = 0;
 		}
-		a->nsamps += a->bsize;
-		while (a->nsamps >= a->fsize)
+		a->iainidx = j;
+		a->nsamps += bsize;
+		while (a->nsamps >= fsize)
 		{
-			for (i = 0, j = a->iaoutidx; i < a->fsize; i++, j = (j + 1) % a->iasize)
-				a->forfftin[i] = a->window[i] * a->inaccum[j];
-			a->iaoutidx = (a->iaoutidx + a->incr) % a->iasize;
-			a->nsamps -= a->incr;
+			int n1 = iasize - a->iaoutidx;
+			if (n1 > fsize) n1 = fsize;
+			for (i = 0; i < n1; i++)
+				a->forfftin[i] = a->window[i] * a->inaccum[a->iaoutidx + i];
+			for (; i < fsize; i++)
+				a->forfftin[i] = a->window[i] * a->inaccum[i - n1];
+			if ((a->iaoutidx += incr) >= iasize) a->iaoutidx -= iasize;
+			a->nsamps -= incr;
 			fftw_execute (a->Rfor);
 			calc_gain(a);
 			for (i = 0; i < a->msize; i++)
@@ -1225,30 +1252,41 @@ void xemnr (EMNR a, int pos)
 			}
 						post2(a);
 			fftw_execute (a->Rrev);
-			for (i = 0; i < a->fsize; i++)
+			for (i = 0; i < fsize; i++)
 				a->save[a->saveidx][i] = a->window[i] * a->revfftout[i];
-			for (i = a->ovrlp; i > 0; i--)
+			for (i = ovrlp; i > 0; i--)
 			{
-				sbuff = (a->saveidx + i) % a->ovrlp;
-				sbegin = a->incr * (a->ovrlp - i);
-				for (j = sbegin, k = a->oainidx; j < a->incr + sbegin; j++, k = (k + 1) % a->oasize)
+				const double* WDSP_RESTRICT sv;
+				double* WDSP_RESTRICT oa = a->outaccum;
+				int m1;
+				sbuff = (a->saveidx + i) % ovrlp;
+				sbegin = incr * (ovrlp - i);
+				sv = a->save[sbuff] + sbegin;
+				m1 = oasize - a->oainidx;
+				if (m1 > incr) m1 = incr;
+				k = a->oainidx;
+				if (i == ovrlp)
 				{
-					if ( i == a->ovrlp)
-						a->outaccum[k]	= a->save[sbuff][j];
+					for (j = 0; j < m1; j++)		oa[k + j] = sv[j];
+					for (; j < incr; j++)			oa[j - m1] = sv[j];
+				}
 				else
-						a->outaccum[k] += a->save[sbuff][j];
-				}
-			}
-			a->saveidx = (a->saveidx + 1) % a->ovrlp;
-			a->oainidx = (a->oainidx + a->incr) % a->oasize;
-		}
-		for (i = 0; i < a->bsize; i++)
 				{
-			a->out[2 * i + 0] = a->outaccum[a->oaoutidx];
-			a->out[2 * i + 1] = 0.0;
-			a->oaoutidx = (a->oaoutidx + 1) % a->oasize;
+					for (j = 0; j < m1; j++)		oa[k + j] += sv[j];
+					for (; j < incr; j++)			oa[j - m1] += sv[j];
 				}
 			}
+			if (++a->saveidx == ovrlp) a->saveidx = 0;
+			if ((a->oainidx += incr) >= oasize) a->oainidx -= oasize;
+		}
+		for (i = 0, k = a->oaoutidx; i < bsize; i++)
+		{
+			a->out[2 * i + 0] = a->outaccum[k];
+			a->out[2 * i + 1] = 0.0;
+			if (++k == oasize) k = 0;
+		}
+		a->oaoutidx = k;
+	}
 	else if (a->out != a->in)
 		memcpy (a->out, a->in, a->bsize * sizeof (complex));
 }
@@ -185,6 +185,7 @@ typedef struct _emnr
 		double zetaThresh;
 		double psi;
 		double* nmask;
+		double* csum;		// prefix sums of mask[], msize + 1 entries
 		double t2;
 	} ae;
 		struct _post2
@@ -0,0 +1,74 @@
+/*	fastmath.h
+
+This file is part of a program that implements a Software-Defined Radio.
+
+Copyright (C) 2013, 2024, 2025 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 _fastmath_h
+#define _fastmath_h
+
+#include <stdint.h>
+#include <string.h>
+
+/* log10 for strictly positive, finite, normal arguments; about 2.4x the
+   throughput of libm's.
+
+   Accurate to within 2.0e-13 absolute over (0, 1], which is far tighter than
+   mlog10() in meterlog10.c (2.1e-4 absolute -- a 2048-entry table with no
+   interpolation).  That table is fine for driving a dB meter display, but too
+   coarse where the result feeds arithmetic rather than a readout.
+
+   Caller must guarantee x > 0 and normal.  Zero, negatives, denormals,
+   infinities and NaN are not handled.  */
+static inline double wdsp_log10 (double x)
+{
+	uint64_t bits;
+	double m, s, s2, p;
+	int e;
+
+	memcpy (&bits, &x, sizeof (bits));
+	e = (int)((bits >> 52) & 0x7FF) - 1023;
+	/* clear the exponent field, leaving the mantissa in [1, 2) */
+	bits = (bits & 0x000FFFFFFFFFFFFFULL) | 0x3FF0000000000000ULL;
+	memcpy (&m, &bits, sizeof (m));
+
+	/* Recentre onto [sqrt(1/2), sqrt(2)) so the series stays in its
+	   fast-converging range; |s| <= 0.1716 afterwards. */
+	if (m > 1.4142135623730951)
+	{
+		m *= 0.5;
+		e += 1;
+	}
+
+	/* log(m) = 2 * atanh(s), s = (m-1)/(m+1) */
+	s = (m - 1.0) / (m + 1.0);
+	s2 = s * s;
+	p = 2.0 * (s + s * s2 * (3.3333333333333331e-01 + s2 * (2.0000000000000001e-01
+		+ s2 * (1.4285714285714285e-01 + s2 * (1.1111111111111110e-01
+		+ s2 * (9.0909090909090912e-02 + s2 * 7.6923076923076927e-02))))));
+
+	/* log10(x) = (log(m) + e * ln2) / ln10 */
+	return (p + (double)e * 6.9314718055994531e-01) * 4.3429448190325182e-01;
+}
+
+#endif
@@ -616,19 +616,31 @@ void xphrot (PHROT a)
 	if (a->run)
 	{
 		int i, n;
-		for (i = 0; i < a->size; i++)
+		const int size = a->size;
+		const int nstages = a->nstages;
+		const double b0 = a->b0, b1 = a->b1, a1 = a->a1;
+		/* in and out are the same buffer in TXA, so neither may be restrict */
+		const double* in = a->in;
+		double* out = a->out;
+		/* x0[]/y0[] never carried state between samples: x0[n] was only ever the
+		   previous stage's output and y0[n] this stage's. Keep that single value
+		   in a register and cascade it, leaving x1/y1 as the actual filter state.
+		   Stores through out[] could alias the struct's doubles, so hoist the
+		   coefficients too. */
+		double* WDSP_RESTRICT x1 = a->x1;
+		double* WDSP_RESTRICT y1 = a->y1;
+
+		for (i = 0; i < size; i++)
 		{
-			a->x0[0] = a->in[2 * i + 0];
-			for (n = 0; n < a->nstages; n++)
+			double v = in[2 * i + 0];
+			for (n = 0; n < nstages; n++)
 			{
-				if (n > 0) a->x0[n] = a->y0[n - 1];
-				a->y0[n]	= a->b0 * a->x0[n]
-							+ a->b1 * a->x1[n]
-							- a->a1 * a->y1[n];
-				a->y1[n] = a->y0[n];
-				a->x1[n] = a->x0[n];
+				double y = b0 * v + b1 * x1[n] - a1 * y1[n];
+				x1[n] = v;
+				y1[n] = y;
+				v = y;
 			}
-			a->out[2 * i + 0] = a->y0[a->nstages - 1];
+			out[2 * i + 0] = v;
 		}
 	}
 	else if (a->out != a->in)
@@ -39,9 +39,17 @@ john.d.melton@googlemail.com
 #if defined(linux) || defined(__APPLE__)

 void QueueUserWorkItem(void *function,void *context,int flags) {
-	pthread_t t;
-	pthread_create(&t, NULL, function, context);
-	pthread_join(t, NULL);
+	//
+	// The Windows call queues the work item on a thread pool and returns at
+	// once, so callers get their items run in parallel. This shim spawned a
+	// thread and immediately joined it, which is a plain synchronous call that
+	// happens to cost a thread creation (~15 us) and delivers no parallelism.
+	// Call the function directly: same ordering, no thread.
+	//
+	// The cast matches the one pthread_create() performed here before.
+	//
+	(void)flags;
+	((void *(*)(void *))function)(context);
 }

 static inline void init_crit_section(pthread_mutex_t *mutex) {
@@ -32,6 +32,37 @@ warren@wpratt.com
 *																								*
 ************************************************************************************************/

+/* Accumulate n taps of a unit-stride complex dot product into *pI / *pQ.
+
+   Four independent accumulator pairs are carried so the FMAs are not serialized
+   on a single dependency chain, and so the compiler is free to vectorize: a
+   plain 'I += h[j]*x[j]' reduction cannot be reassociated without -ffast-math,
+   which this library must not enable (it relies on IEEE semantics for 0/0 = NaN
+   and x/0 = Inf).  Summation order therefore differs from a strict left-to-right
+   reduction, at the usual pairwise-summation accuracy gain. */
+static inline void resample_dot (const double* WDSP_RESTRICT hp,
+	const double* WDSP_RESTRICT xI, const double* WDSP_RESTRICT xQ,
+	int n, double* pI, double* pQ)
+{
+	double i0 = 0.0, i1 = 0.0, i2 = 0.0, i3 = 0.0;
+	double q0 = 0.0, q1 = 0.0, q2 = 0.0, q3 = 0.0;
+	int j = 0;
+	for (; j <= n - 4; j += 4)
+	{
+		i0 += hp[j + 0] * xI[j + 0];	q0 += hp[j + 0] * xQ[j + 0];
+		i1 += hp[j + 1] * xI[j + 1];	q1 += hp[j + 1] * xQ[j + 1];
+		i2 += hp[j + 2] * xI[j + 2];	q2 += hp[j + 2] * xQ[j + 2];
+		i3 += hp[j + 3] * xI[j + 3];	q3 += hp[j + 3] * xQ[j + 3];
+	}
+	for (; j < n; j++)
+	{
+		i0 += hp[j] * xI[j];
+		q0 += hp[j] * xQ[j];
+	}
+	*pI += (i0 + i1) + (i2 + i3);
+	*pQ += (q0 + q1) + (q2 + q3);
+}
+
 void calc_resample (RESAMPLE a)
 {
 	int x, y, z;
@@ -71,7 +102,8 @@ void calc_resample (RESAMPLE a)
 		for (k = 0; k < a->ncoef; k += a->L)
 			a->h[i++] = impulse[j + k];
 	a->ringsize = a->cpp;
-	a->ring = (double *)malloc0(a->ringsize * sizeof(complex));
+	a->ringI = (double *)malloc0(a->ringsize * sizeof(double));
+	a->ringQ = (double *)malloc0(a->ringsize * sizeof(double));
 	a->idx_in = a->ringsize - 1;
 	a->phnum = 0;
 	_aligned_free(impulse);
@@ -79,7 +111,8 @@ void calc_resample (RESAMPLE a)

 void decalc_resample (RESAMPLE a)
 {
-	_aligned_free(a->ring);
+	_aligned_free(a->ringQ);
+	_aligned_free(a->ringI);
 	_aligned_free(a->h);
 }

@@ -112,7 +145,8 @@ void destroy_resample (RESAMPLE a)
 PORT
 void flush_resample (RESAMPLE a)
 {
-	memset (a->ring, 0, a->ringsize * sizeof (complex));
+	memset (a->ringI, 0, a->ringsize * sizeof (double));
+	memset (a->ringQ, 0, a->ringsize * sizeof (double));
 	a->idx_in = a->ringsize - 1;
 	a->phnum = 0;
 }
@@ -123,40 +157,49 @@ int xresample (RESAMPLE a)
 	int outsamps = 0;
 	if (a->run)
 	{
-		int i, j, n;
-		int idx_out;
+		int i, n1;
 		double I, Q;

-		int cpp = a->cpp;
+		const int cpp = a->cpp;
+		const int ringsize = a->ringsize;
+		const int L = a->L;
+		const int M = a->M;
+		const int size = a->size;
+		const double* WDSP_RESTRICT h = a->h;
+		const double* WDSP_RESTRICT in = a->in;
+		double* WDSP_RESTRICT ringI = a->ringI;
+		double* WDSP_RESTRICT ringQ = a->ringQ;
+		double* WDSP_RESTRICT out = a->out;
 		int idx_in = a->idx_in;
-		int ringsize = a->ringsize;
-		double* h = a->h;
-		double* ring = a->ring;
+		int phnum = a->phnum;

-		for (i = 0; i < a->size; i++)
+		for (i = 0; i < size; i++)
 		{
-			ring[2 * idx_in + 0] = a->in[2 * i + 0];
-			ring[2 * idx_in + 1] = a->in[2 * i + 1];
-			while (a->phnum < a->L)
+			ringI[idx_in] = in[2 * i + 0];
+			ringQ[idx_in] = in[2 * i + 1];
+			while (phnum < L)
 			{
+				const double* WDSP_RESTRICT hp = h + cpp * phnum;
+				/* The tap loop walks the ring forward from idx_in and wraps at
+				   most once.  Split it at the wrap point so both halves are
+				   unit-stride: the wrap test that used to sit inside the loop
+				   made the address non-affine and blocked vectorization. */
+				if ((n1 = ringsize - idx_in) > cpp) n1 = cpp;
 				I = 0.0;
 				Q = 0.0;
-				n = cpp * a->phnum;
-				for (j = 0; j < cpp; j++)
-				{
-					if ((idx_out = idx_in + j) >= ringsize) idx_out -= ringsize;
-					I += h[n + j] * ring[2 * idx_out + 0];
-					Q += h[n + j] * ring[2 * idx_out + 1];
-				}
-				a->out[2 * outsamps + 0] = I;
-				a->out[2 * outsamps + 1] = Q;
+				resample_dot (hp, ringI + idx_in, ringQ + idx_in, n1, &I, &Q);
+				if (n1 < cpp)
+					resample_dot (hp + n1, ringI, ringQ, cpp - n1, &I, &Q);
+				out[2 * outsamps + 0] = I;
+				out[2 * outsamps + 1] = Q;
 				outsamps++;
-				a->phnum += a->M;
+				phnum += M;
 			}
-			a->phnum -= a->L;
-			if (--idx_in < 0) idx_in = a->ringsize - 1;
+			phnum -= L;
+			if (--idx_in < 0) idx_in = ringsize - 1;
 		}
 		a->idx_in = idx_in;
+		a->phnum = phnum;
 	}
 	else if (a->in != a->out)
 		memcpy (a->out, a->in, a->size * sizeof (complex));
@@ -240,6 +283,24 @@ void destroy_resampleV (void* ptr)
 *																								*
 ************************************************************************************************/

+/* Real-valued counterpart of resample_dot(). */
+static inline void resampleF_dot (const double* WDSP_RESTRICT hp,
+	const double* WDSP_RESTRICT x, int n, double* pI)
+{
+	double i0 = 0.0, i1 = 0.0, i2 = 0.0, i3 = 0.0;
+	int j = 0;
+	for (; j <= n - 4; j += 4)
+	{
+		i0 += hp[j + 0] * x[j + 0];
+		i1 += hp[j + 1] * x[j + 1];
+		i2 += hp[j + 2] * x[j + 2];
+		i3 += hp[j + 3] * x[j + 3];
+	}
+	for (; j < n; j++)
+		i0 += hp[j] * x[j];
+	*pI += (i0 + i1) + (i2 + i3);
+}
+
 RESAMPLEF create_resampleF ( int run, int size, float* in, float* out, int in_rate, int out_rate)
 {
 	RESAMPLEF a = (RESAMPLEF) malloc0 (sizeof (resampleF));
@@ -305,31 +366,46 @@ int xresampleF (RESAMPLEF a)
 	int outsamps = 0;
 	if (a->run)
 	{
-		int i, j, n;
-		int idx_out;
+		int i;
 		double I;

-		for (i = 0; i < a->size; i++)
-		{
-			a->ring[a->idx_in] = (double)a->in[i];
+		const int cpp = a->cpp;
+		const int ringsize = a->ringsize;
+		const int L = a->L;
+		const int M = a->M;
+		const int size = a->size;
+		const double* WDSP_RESTRICT h = a->h;
+		const float* WDSP_RESTRICT in = a->in;
+		double* WDSP_RESTRICT ring = a->ring;
+		float* WDSP_RESTRICT out = a->out;
+		int idx_in = a->idx_in;
+		int phnum = a->phnum;
+		int n1;

-			while (a->phnum < a->L)
+		for (i = 0; i < size; i++)
 		{
+			ring[idx_in] = (double)in[i];
+
+			while (phnum < L)
+			{
+				const double* WDSP_RESTRICT hp = h + cpp * phnum;
+				/* see resample_dot(): split at the ring wrap so both halves are
+				   unit-stride, and carry independent accumulators */
+				if ((n1 = ringsize - idx_in) > cpp) n1 = cpp;
 				I = 0.0;
-				n = a->cpp * a->phnum;
-				for (j = 0; j < a->cpp; j++)
-				{
-					if ((idx_out = a->idx_in + j) >= a->ringsize) idx_out -= a->ringsize;
-					I += a->h[n + j] * a->ring[idx_out];
-				}
-				a->out[outsamps] = (float)I;
+				resampleF_dot (hp, ring + idx_in, n1, &I);
+				if (n1 < cpp)
+					resampleF_dot (hp + n1, ring, cpp - n1, &I);
+				out[outsamps] = (float)I;

 				outsamps++;
-				a->phnum += a->M;
+				phnum += M;
 			}
-			a->phnum -= a->L;
-			if (--a->idx_in < 0) a->idx_in = a->ringsize - 1;
+			phnum -= L;
+			if (--idx_in < 0) idx_in = ringsize - 1;
 		}
+		a->idx_in = idx_in;
+		a->phnum = phnum;
 	}
 	else if (a->in != a->out)
 		memcpy (a->out, a->in, a->size * sizeof (float));
@@ -52,7 +52,9 @@ typedef struct _resample
 	int M;				// decimation factor
 	double* h;			// coefficients
 	int ringsize;		// number of complex pairs the ring buffer holds
-	double* ring;		// ring buffer
+	double* ringI;		// ring buffer, in-phase
+	double* ringQ;		// ring buffer, quadrature (split from I so the tap loop
+						//   reads unit-stride and vectorizes)
 	int cpp;			// coefficients of the phase
 	int phnum;			// phase number
 } resample, *RESAMPLE;
@@ -59,10 +59,27 @@ void calc_varsamp (VARSAMP a)
 		fc_norm_low = a->fc_low / norm_rate;
 	a->rsize = (int)(140.0 * norm_rate / min_rate);
 	a->ncoef = a->rsize + 1;
-	a->ncoef += (a->R - 1) * (a->ncoef - 1);
-	a->h = fir_bandpass(a->ncoef, fc_norm_low, fc_norm_high, (double)a->R, 1, 0, (double)a->R * a->gain);
-	// print_impulse ("imp.txt", a->ncoef, a->h, 0, 0);
-	a->ring = (double *)malloc0(a->rsize * sizeof(complex));
+	a->ncoef += (a->R - 1) * (a->ncoef - 1);		// = R * rsize + 1
+	{
+		/* Store the coefficients transposed into phases.  hshift() walks
+		   h[hidx + m*R] for m = 0..rsize-1, which strides by R doubles -- 8 KB
+		   at R = 1024 -- over a 1.1 MB table, so every tap is its own cache
+		   line.  Transposing makes each phase contiguous; hshift() interpolates
+		   between phases hidx and hidx+1, hence R+1 of them. */
+		int p, m;
+		const int R = a->R, rsize = a->rsize;
+		double* h = fir_bandpass(a->ncoef, fc_norm_low, fc_norm_high, (double)R, 1, 0, (double)R * a->gain);
+		// every element is written below, so skip malloc0()'s memset of ~1 MB
+		a->hp = (double *)_aligned_malloc ((size_t)(R + 1) * rsize * sizeof (double), 16);
+		// walk h forward (p is its fast axis) so the prefetcher sees a linear
+		// stream; h is cold here, straight from fir_bandpass()
+		for (m = 0; m < rsize; m++)
+			for (p = 0; p <= R; p++)
+				a->hp[(size_t)p * rsize + m] = h[p + (size_t)m * R];
+		_aligned_free (h);
+	}
+	a->ringI = (double *)malloc0(a->rsize * sizeof(double));
+	a->ringQ = (double *)malloc0(a->rsize * sizeof(double));
 	a->idx_in = a->rsize - 1;
 	a->h_offset = 0.0;
 	a->hs = (double *)malloc0 (a->rsize * sizeof (double));
@@ -72,8 +89,9 @@ void calc_varsamp (VARSAMP a)
 void decalc_varsamp (VARSAMP a)
 {
 	_aligned_free (a->hs);
-	_aligned_free (a->ring);
-	_aligned_free (a->h);
+	_aligned_free (a->ringQ);
+	_aligned_free (a->ringI);
+	_aligned_free (a->hp);
 }

 VARSAMP create_varsamp ( int run, int size, double* in, double* out,
@@ -105,22 +123,60 @@ void destroy_varsamp (VARSAMP a)

 void flush_varsamp (VARSAMP a)
 {
-	memset (a->ring, 0, a->rsize * sizeof (complex));
+	memset (a->ringI, 0, a->rsize * sizeof (double));
+	memset (a->ringQ, 0, a->rsize * sizeof (double));
 	a->idx_in = a->rsize - 1;
 	a->h_offset = 0.0;
 	a->isamps = 0.0;
 }

+/* Accumulate n taps of a unit-stride complex dot product into *pI / *pQ.
+
+   Four independent accumulator pairs keep the FMAs off a single dependency
+   chain and let the vectorizer in: an 'I += h[j]*x[j]' reduction cannot be
+   reassociated without -ffast-math, which this library must not enable (it
+   relies on IEEE semantics for 0/0 = NaN and x/0 = Inf). */
+static inline void varsamp_dot (const double* WDSP_RESTRICT hp,
+	const double* WDSP_RESTRICT xI, const double* WDSP_RESTRICT xQ,
+	int n, double* pI, double* pQ)
+{
+	double i0 = 0.0, i1 = 0.0, i2 = 0.0, i3 = 0.0;
+	double q0 = 0.0, q1 = 0.0, q2 = 0.0, q3 = 0.0;
+	int j = 0;
+	for (; j <= n - 4; j += 4)
+	{
+		i0 += hp[j + 0] * xI[j + 0];	q0 += hp[j + 0] * xQ[j + 0];
+		i1 += hp[j + 1] * xI[j + 1];	q1 += hp[j + 1] * xQ[j + 1];
+		i2 += hp[j + 2] * xI[j + 2];	q2 += hp[j + 2] * xQ[j + 2];
+		i3 += hp[j + 3] * xI[j + 3];	q3 += hp[j + 3] * xQ[j + 3];
+	}
+	for (; j < n; j++)
+	{
+		i0 += hp[j] * xI[j];
+		q0 += hp[j] * xQ[j];
+	}
+	*pI += (i0 + i1) + (i2 + i3);
+	*pQ += (q0 + q1) + (q2 + q3);
+}
+
 void hshift (VARSAMP a)
 {
-	int i, j, k;
+	int m;
 	int hidx;
 	double frac, pos;
+	const int rsize = a->rsize;
+	const double* WDSP_RESTRICT h0;
+	const double* WDSP_RESTRICT h1;
+	double* WDSP_RESTRICT hs = a->hs;
+	/* h_offset is normalized to [0,1) by the caller, so hidx is in [0, R-1]
+	   and phase hidx+1 <= R exists. */
 	pos = (double)a->R * a->h_offset;
 	hidx = (int)(pos);
 	frac = pos - (double)hidx;
-	for (i = a->rsize - 1, j = hidx, k = hidx + 1; i >= 0; i--, j += a->R, k += a->R)
-		a->hs[i] = a->h[j] + frac * (a->h[k] - a->h[j]);
+	h0 = a->hp + (size_t)hidx * rsize;
+	h1 = h0 + rsize;
+	for (m = 0; m < rsize; m++)
+		hs[rsize - 1 - m] = h0[m] + frac * (h1[m] - h0[m]);
 }

 int xvarsamp (VARSAMP a, double var)
@@ -140,13 +196,21 @@ int xvarsamp (VARSAMP a, double var)
 	else			a->dicvar = 0.0;
 	if (a->run)
 	{
-		int i, j;
-		int idx_out;
+		int i, n1;
 		double I, Q;
+		const int rsize = a->rsize;
+		/* a->hs is rewritten by hshift() on every output sample, so it must not
+		   be hoisted behind a restrict pointer here; varsamp_dot() re-reads it. */
+		const double* in = a->in;
+		double* out = a->out;
+		double* WDSP_RESTRICT ringI = a->ringI;
+		double* WDSP_RESTRICT ringQ = a->ringQ;
+		int idx_in = a->idx_in;
+
 		for (i = 0; i < a->size; i++)
 		{
-			a->ring[2 * a->idx_in + 0] = a->in[2 * i + 0];
-			a->ring[2 * a->idx_in + 1] = a->in[2 * i + 1];
+			ringI[idx_in] = in[2 * i + 0];
+			ringQ[idx_in] = in[2 * i + 1];
 			a->inv_cvar += a->dicvar;
 			picvar = (uint64_t*)(&a->inv_cvar);
 			N = *picvar & 0xffffffffffff0000;
@@ -160,20 +224,21 @@ int xvarsamp (VARSAMP a, double var)
 				a->h_offset += a->delta;
 				while (a->h_offset >= 1.0) a->h_offset -= 1.0;
 				while (a->h_offset <  0.0) a->h_offset += 1.0;
-				for (j = 0; j < a->rsize; j++)
-				{
-					if ((idx_out = a->idx_in + j) >= a->rsize) idx_out -= a->rsize;
-					I += a->hs[j] * a->ring[2 * idx_out + 0];
-					Q += a->hs[j] * a->ring[2 * idx_out + 1];
-				}
-				a->out[2 * outsamps + 0] = I;
-				a->out[2 * outsamps + 1] = Q;
+				/* the ring wraps at most once over rsize taps; split it so both
+				   halves are unit-stride */
+				n1 = rsize - idx_in;
+				varsamp_dot (a->hs, ringI + idx_in, ringQ + idx_in, n1, &I, &Q);
+				if (n1 < rsize)
+					varsamp_dot (a->hs + n1, ringI, ringQ, rsize - n1, &I, &Q);
+				out[2 * outsamps + 0] = I;
+				out[2 * outsamps + 1] = Q;
 				outsamps++;
 				a->isamps += a->inv_cvar;
 			}
 			a->isamps -= 1.0;
-			if (--a->idx_in < 0) a->idx_in = a->rsize - 1;
+			if (--idx_in < 0) idx_in = rsize - 1;
 		}
+		a->idx_in = idx_in;
 	}
 	else if (a->in != a->out)
 		memcpy (a->out, a->in, a->size * sizeof (complex));
@@ -41,9 +41,13 @@ typedef struct _varsamp
 	double gain;
 	int idx_in;
 	int ncoef;
-	double* h;
+	double* hp;			// coefficients, polyphase: hp[p * rsize + m] = h[p + m * R],
+						//   p = 0..R.  hshift() then reads two adjacent phases
+						//   contiguously instead of striding by R.
 	int rsize;
-	double* ring;
+	double* ringI;		// ring buffer, in-phase
+	double* ringQ;		// ring buffer, quadrature (split from I so the tap loop
+						//   reads unit-stride and vectorizes)
 	double var;
 	int varmode;
 	double cvar;
@@ -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);
@@ -517,6 +518,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
 //
@@ -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);
+}
@@ -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
@@ -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);
+}
@@ -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
Author	SHA1	Message	Date
vladimir	544e75798c	Merge branch 'wfm' into cpu-optimization	2026-07-18 10:51:50 +03:00
vladimir	3ad40a00da	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 <noreply@anthropic.com>	2026-07-18 10:47:44 +03:00
Uladzimir Karpenka	c1964817dd	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 <noreply@anthropic.com>	2026-07-15 11:03:03 +03:00
Uladzimir Karpenka	81f7750def	Merge remote-tracking branch 'origin/main' into cpu-optimization	2026-07-15 10:50:34 +03:00
Uladzimir Karpenka	d24b5ae316	Merge remote-tracking branch 'origin/main' into wfm	2026-07-15 10:50:15 +03:00
Uladzimir Karpenka	65cb3c386e	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 <noreply@anthropic.com>	2026-07-10 07:43:07 +03:00
Uladzimir Karpenka	f39eea6b01	phrot, cfcomp: cascade in a register; drop idiv from cfcomp's ring walks xphrot() kept x0[] and y0[] in the struct, but neither carried state between samples: x0[n] was always the previous stage's output and y0[n] this stage's. Cascade that single value in a register and leave x1/y1 as the real filter state. Stores through out[] may alias the struct's doubles, so hoist the coefficients as well. Eight first-order sections then run without touching memory except for the state. xcfcomp() is structurally the same overlap-add loop as xemnr(), and had the same defect: four ring indices advanced with '% size' per step, ~5100 integer divisions per call at fsize = 2048, which a profile showed dominating the block (1820 samples in xcfcomp against 362 in calc_mask and ~520 in the FFTs). The indices step by one, and iasize >= fsize and oasize >= incr always hold, so they wrap at most once per loop: walk contiguous runs and wrap between them. Measured in situ on an Apple M1 Pro, 512-sample buffers, cost of turning the block on, best of 5: phrot 17219 ns -> 7047 ns 2.44x cfcomp 28918 ns -> 13637 ns 2.12x Both are bit-identical: phrot preserves the operation order, and cfcomp only changes integer index arithmetic. The RX chain is unchanged bit-for-bit. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>	2026-07-10 00:31:32 +03:00
Uladzimir Karpenka	91df5b1f2d	varsamp: build the phase table without malloc0's memset, read h forward calc_varsamp() gained a transpose of the coefficient table, which showed up as ~0.3 ms on create_varsamp()/setInRate_varsamp() (1.0 ms -> 1.4 ms). It is a one-off setup cost, not on the sample path, but it is easy to trim: every element of hp is written, so malloc0()'s memset of ~1 MB is dead, and h is cold straight out of fir_bandpass(), so walk it along its fast axis and let the prefetcher work. Output is bit-identical; xvarsamp() is unchanged at ~34 us/buffer. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>	2026-07-10 00:13:27 +03:00
Uladzimir Karpenka	12dc701604	analyzer: block the dispatcher on a semaphore, stop spawning a thread per frame Two pieces of pure overhead, neither of them doing any DSP. The POSIX QueueUserWorkItem() shim spawned a thread and immediately joined it. That is a synchronous call -- the parallelism the Windows thread pool provides is absent here either way -- so it only bought a pthread_create() and its stack mmap, ~15 us, for every FFT frame. Call the function directly. (The Windows build is untouched and still gets its thread pool; a real pool for POSIX would be a separate change, and only pays off for num_stitch > 1.) The dispatcher thread ran `for (each ss, LO) {...} Sleep(1);`, and Sleep(1) is usleep(1000), so it woke 1000 times a second to re-read the same flags. That burned 0.63% of a core even with no samples arriving. Give it a semaphore instead, signalled by the four Spectrum*() entry points when new samples land, and by SetAnalyzer/DestroyAnalyzer after they raise end_dispatcher so the blocking wait always has a way out. SetAnalyzer holds SetAnalyzerSection while it waits for the dispatcher to quit, and the dispatcher never takes that section, so signalling from under it is safe. Measured on an Apple M1 Pro; 16384-point complex FFT, 1024-sample buffers, ~40 pixel frames/s, CPU of all threads via getrusage: idle, no samples at all 0.63% of a core -> 0.00% under load 3.54% -> 2.65% Output pixels are bit-identical. Stressed with 3 create/destroy cycles and 24 on-the-fly SetAnalyzer reconfigurations while a second thread fed samples: no deadlock. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>	2026-07-09 23:59:20 +03:00
Uladzimir Karpenka	645cbbb2d1	varsamp: transpose the coefficients into phases, vectorize the tap loop hshift() rebuilds the whole interpolated tap set on every output sample, reading h[hidx + mR] for m = 0..rsize-1. With R = 1024 (what rmatch asks for) that strides 8 KB at a time through a 1.1 MB table, so every one of the 2rsize reads is its own cache line. It cost more than the filter it was feeding: 258 ns per output sample against 146 ns for the tap loop. Store the coefficients transposed instead, hp[prsize + m] = h[p + mR], so the two phases hshift() interpolates between are each contiguous. R+1 phases are needed since it reads hidx and hidx+1, and h_offset is kept in [0,1) by the caller so hidx <= R-1. The untransposed h is freed; the impulse cache hands back a copy, so varsamp owns it. Net memory is unchanged. The tap loop had the same wrap test per tap as resample.c did, so split it at the wrap and carry four independent accumulator pairs; the ring is split into I/Q so the taps load unit-stride. Note a->hs is rewritten by hshift() inside the sample loop, so it must not be hoisted behind a restrict pointer in xvarsamp(). Measured on an Apple M1 Pro, 512-sample buffers, best of 5: 48k -> 48k varmode=0 138580 ns -> 37370 ns 3.71x 48k -> 48k varmode=1 138390 ns -> 32333 ns 4.28x 48k -> 44.1k varmode=1 141473 ns -> 38220 ns 3.70x 44.1k -> 48k varmode=1 152307 ns -> 45263 ns 3.36x hshift() is numerically identical -- same coefficients, different layout. Only the reassociated tap sum rounds differently: worst deviation 7.2e-16 over four rate configurations, an SNR of 344 dB. Driven end to end through rmatch's public API, output SNR is 306 dB and total energy matches. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>	2026-07-09 23:19:36 +03:00
Uladzimir Karpenka	fd2ba84e7d	emnr: linear-time aepf, drop idiv from ring walks, halve getKey's logs Three separate costs, found with a sampling profile of the RX chain with emnr forced on: aepf() averaged mask[] over a window of N = 2psi + 1 = 41 bins by walking the window for every one of the 2049 bins, i.e. O(msizeN). Its three spans are all symmetric windows clipped at the array ends, so take each from a prefix sum instead: one subtraction per output, O(msize). xemnr() advanced four ring indices with a '% size' per step. iasize is 4096 and oasize 1024 here, and neither is known to the compiler, so each step was a real integer division -- ~8700 of them per frame. The indices step by one and, since iasize >= fsize and oasize >= incr always hold, wrap at most once per loop, so walk contiguous runs and wrap between them. calc_gain() called getKey() twice per bin with the same gamma, so the gamma row index and its log10 were computed twice. Split getKey into keyIndex() + keyLerp() and locate gamma once. The remaining logs go through wdsp_log10() (new fastmath.h), accurate to 2e-13 against libm and ~2.4x its throughput; gamma and xi are bracketed against the table limits first, so the argument is always positive and normal. Also clamp the row index so the second bilinear corner cannot address the next row of the 241x241 table. Measured in situ on an Apple M1 Pro, 512-sample buffers, cost of turning emnr on, best of 5: baseline 73661 ns + aepf, ring walks 48927 ns 1.51x + getKey 37852 ns 1.95x Output is not bit-identical, as the prefix sum and the reassociated logs round differently. Over 300 buffers with emnr alone the worst deviation is 4.0e-09, an SNR of 196 dB; perturbing a single input sample of the unmodified code by one ulp diverges it from itself by 1.4e-08 (186 dB), so this change disturbs the chain less than the last bit of the input does. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>	2026-07-09 23:08:20 +03:00
Uladzimir Karpenka	3faad14fc1	anr, anf: vectorize the LMS tap and update loops Both loops indexed the delay line as (in_idx + j + delay) & mask, one masked index per tap, which made the address non-affine and stopped the vectorizer. The window wraps at most once, so split it at the wrap and walk two contiguous runs instead. As in resample.c, the y/sigma reduction cannot be reassociated without -ffast-math (which this library must not enable, see linux_port.h), so carry four independent accumulator pairs to break the FMA dependency chain and let the vectorizer in. in_buff and out_buff alias in RXA -- both are midbuff -- so only the private d/w arrays are marked restrict. Measured in situ on an Apple M1 Pro, 512-sample buffers, cost of turning the block on, best of 5: anr 59303 ns -> 20253 ns 2.93x anf 55073 ns -> 20511 ns 2.68x Summation order changes, so output is not bit-identical: over 300 buffers of the full RX chain the worst deviation is 2.4e-07, an SNR of 153 dB. An LMS filter is an adaptive feedback loop, so its trajectory is chaotic. As a control, perturbing a single input sample of the unmodified code by one ulp diverges it from itself by 6.6e-07, an SNR of 144.6 dB -- i.e. this change disturbs the filter less than the last bit of the input does. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>	2026-07-09 22:52:33 +03:00
Uladzimir Karpenka	9960a1d0e3	resample: vectorize the polyphase tap loop The tap loop wrapped the ring with a test on every tap: if ((idx_out = idx_in + j) >= ringsize) idx_out -= ringsize; which made the address non-affine and stopped the vectorizer. Split the walk at the wrap point instead, so both halves are unit-stride, and split the complex ring into separate I/Q arrays so the taps load contiguously rather than through a de-interleaving ld2. The dot product also could not be vectorized as written: reassociating an fp reduction needs -ffast-math, which this library must not enable (it relies on IEEE semantics for 0/0 = NaN and x/0 = Inf, see linux_port.h). Carry four independent accumulator pairs instead, which both breaks the FMA dependency chain and lets the vectorizer in on any compiler. Measured on an Apple M1 Pro, 512-sample DSP buffers, best of 5: xresample, decimation to 48 kHz before after speedup 192k -> 48k (561 taps) 342.0 us 84.2 us 4.06x 384k -> 48k (1121 taps) 708.9 us 175.1 us 4.05x 576k -> 48k (1681 taps) 1074.4 us 266.9 us 4.03x 768k -> 48k (2241 taps) 1438.9 us 360.0 us 4.00x full xrxa() chain, 576k input 1148.2 us 337.6 us 3.40x 10.76% 3.16% of one core xresampleF (float, host audio) 2.6x - 3.4x Summation order changes, so the double path is not bit-identical: over 400 buffers of the full RX chain the worst deviation is 1.1e-12, an SNR of 251 dB. The float path is bit-identical, as the cast to float absorbs the difference. With the resampler bypassed (48k in, 48k out) the chain is unchanged bit-for-bit. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>	2026-07-09 22:44:08 +03:00