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
71 changed files with 2202 additions and 277 deletions
@@ -7,6 +7,10 @@ java/build/

 # third_party — only fftw is downloaded manually, rnnoise and libspecbleach are tracked
 third_party/fftw/
+third_party/**/*.o
+third_party/**/*.a
+third_party/fftw-3.3.11.tar.gz
+third_party/fftw-3.3.11/

 # macOS
 .DS_Store
@@ -144,6 +144,8 @@ utilities.c \
 varsamp.c \
 version.c \
 wcpAGC.c \
+wfmd.c \
+wfmmod.c \
 wisdom.c \
 zetaHat.c

@@ -5,17 +5,39 @@
 # Optional overrides:
 #   ANDROID_NDK=/path/to/ndk ANDROID_API=24 ANDROID_ABIS="arm64-v8a x86_64"
 #
-# Requires FFTW source in third_party/fftw/
-#   Download from https://www.fftw.org/download.html and extract so that
-#   third_party/fftw/configure exists.
+# FFTW source is downloaded and extracted automatically.
+# Override FFTW_SRC to use an existing directory.

-ANDROID_NDK ?= /home/vladimir/Android/Sdk/ndk/29.0.14206865
+# Auto-detect NDK path (can always be overridden by setting ANDROID_NDK explicitly):
+#   1. ANDROID_NDK_HOME environment variable
+#   2. Latest NDK under ANDROID_SDK_ROOT or ANDROID_HOME
+#   3. Default SDK locations: ~/Android/Sdk (Linux) or ~/Library/Android/sdk (macOS)
+ifeq ($(origin ANDROID_NDK),undefined)
+  ifdef ANDROID_NDK_HOME
+    ANDROID_NDK := $(ANDROID_NDK_HOME)
+  else
+    _SDK_ROOT := $(or $(ANDROID_SDK_ROOT),$(ANDROID_HOME),\
+                   $(wildcard $(HOME)/Android/Sdk),\
+                   $(wildcard $(HOME)/Library/Android/sdk))
+    ifdef _SDK_ROOT
+      ANDROID_NDK := $(lastword $(sort $(wildcard $(_SDK_ROOT)/ndk/*)))
+    endif
+  endif
+endif
+ifndef ANDROID_NDK
+  $(error Cannot find Android NDK. Set ANDROID_NDK, ANDROID_NDK_HOME, ANDROID_SDK_ROOT, or ANDROID_HOME)
+endif
 ANDROID_API ?= 24
 ANDROID_ABIS ?= arm64-v8a armeabi-v7a x86_64
 ANDROID_HOST_TAG ?= $(shell uname -s | tr '[:upper:]' '[:lower:]' | sed 's/darwin/darwin/;s/linux/linux/')-$(shell uname -m | sed 's/aarch64/arm64/;s/x86_64/x86_64/')

-JBR_BIN ?= /opt/android-studio/jbr/bin
-JAVAC ?= $(JBR_BIN)/javac
+# Auto-detect javac: JAVA_HOME > system PATH > Android Studio JBR (Linux/macOS)
+_JAVAC_CANDIDATES := \
+  $(if $(JAVA_HOME),$(wildcard $(JAVA_HOME)/bin/javac)) \
+  $(shell command -v javac 2>/dev/null) \
+  $(wildcard /opt/android-studio/jbr/bin/javac) \
+  $(wildcard /Applications/Android\ Studio.app/Contents/jbr/Contents/Home/bin/javac)
+JAVAC ?= $(firstword $(_JAVAC_CANDIDATES))

 TOOLCHAIN := $(ANDROID_NDK)/toolchains/llvm/prebuilt/$(ANDROID_HOST_TAG)

@@ -24,7 +46,10 @@ COMMON_CPPFLAGS ?= -I. -I third_party/rnnoise/include -I third_party/libspecblea
 ANDROID_JNI_CFLAGS ?= -std=gnu89 -Wno-implicit-function-declaration -Wno-int-conversion \
  -Wno-incompatible-pointer-types -Wno-incompatible-pointer-types-discards-qualifiers

-FFTW_SRC ?= third_party/fftw
+FFTW_VERSION  := 3.3.11
+FFTW_TAR      := fftw-$(FFTW_VERSION).tar.gz
+FFTW_URL      := https://fftw.org/pub/fftw/$(FFTW_TAR)
+FFTW_SRC      ?= third_party/fftw-$(FFTW_VERSION)
 FFTW_MAKEJOBS ?= $(shell nproc 2>/dev/null || sysctl -n hw.ncpu 2>/dev/null || echo 4)

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

@@ -148,9 +175,9 @@ ANDROID_OBJ_ROOT  := obj/android
 ANDROID_LIB_ROOT  := lib/android
 ANDROID_JAVA_OUT  := lib/android/java

-.PHONY: all clean java-classes check-fftw-src
+.PHONY: all clean distclean java-classes

-all: check-fftw-src \
+all: $(FFTW_SRC)/configure \
  $(foreach abi,$(ANDROID_ABIS),\
    $(ANDROID_LIB_ROOT)/$(abi)/libfftw3.so \
    $(ANDROID_LIB_ROOT)/$(abi)/libfftw3f.so \
@@ -158,15 +185,21 @@ all: check-fftw-src \
    $(ANDROID_LIB_ROOT)/$(abi)/libwdspj.so) \
  java-classes

-check-fftw-src:
-	@test -f $(FFTW_SRC)/configure || { \
-	  echo ""; \
-	  echo "ERROR: FFTW source not found at $(FFTW_SRC)/configure"; \
-	  echo "Download and extract FFTW from https://www.fftw.org/download.html"; \
-	  echo "so that $(FFTW_SRC)/configure exists."; \
-	  echo ""; \
-	  exit 1; \
-	}
+# ── Автоматическое получение FFTW ────────────────────────────────────────────
+
+third_party/$(FFTW_TAR):
+	mkdir -p third_party
+	@echo ">>> Скачиваем FFTW $(FFTW_VERSION)..."
+	@if command -v wget >/dev/null 2>&1; then \
+	  wget -q --show-progress -O $@ "$(FFTW_URL)"; \
+	else \
+	  curl -L --progress-bar -o $@ "$(FFTW_URL)"; \
+	fi
+
+$(FFTW_SRC)/configure: third_party/$(FFTW_TAR)
+	@echo ">>> Распаковываем FFTW..."
+	tar xf $< -C third_party/
+	@touch $@

 # abi_target: compiler triple for NDK clang wrapper
 define abi_target
@@ -226,11 +259,12 @@ ABI_$(1)_WDSP_OBJS := $$(patsubst %.c,$$(ABI_$(1)_OBJDIR)/%.o,$$(ALL_C_SOURCES))
 ABI_$(1)_JNI_OBJ   := $$(ABI_$(1)_OBJDIR)/$(JNI_SOURCE:.c=.o)

 # --- Build fftw3 (double) ---
-$$(ABI_$(1)_FFTW_STAMP):
+$$(ABI_$(1)_FFTW_STAMP): $(FFTW_SRC)/configure
 	@mkdir -p $$(ABI_$(1)_OBJDIR)/fftw-build $$(ABI_$(1)_FFTW_PREFIX)
 	cd $$(ABI_$(1)_OBJDIR)/fftw-build && \
 	  $(abspath $(FFTW_SRC))/configure \
-	    --host=$(call fftw_host,$(1)) \
+	    --build=$(shell uname -m)-linux-gnu \
+	    --host=$(strip $(call fftw_host,$(1))) \
 	    CC="$$(ABI_$(1)_CC)" \
 	    CFLAGS="-fPIC" \
 	    --prefix=$(abspath $$(ABI_$(1)_FFTW_PREFIX)) \
@@ -240,11 +274,12 @@ $$(ABI_$(1)_FFTW_STAMP):
 	@touch $$@

 # --- Build fftw3f (float/single) ---
-$$(ABI_$(1)_FFTWF_STAMP):
+$$(ABI_$(1)_FFTWF_STAMP): $(FFTW_SRC)/configure
 	@mkdir -p $$(ABI_$(1)_OBJDIR)/fftwf-build $$(ABI_$(1)_FFTWF_PREFIX)
 	cd $$(ABI_$(1)_OBJDIR)/fftwf-build && \
 	  $(abspath $(FFTW_SRC))/configure \
-	    --host=$(call fftw_host,$(1)) \
+	    --build=$(shell uname -m)-linux-gnu \
+	    --host=$(strip $(call fftw_host,$(1))) \
 	    CC="$$(ABI_$(1)_CC)" \
 	    CFLAGS="-fPIC" \
 	    --prefix=$(abspath $$(ABI_$(1)_FFTWF_PREFIX)) \
@@ -311,3 +346,6 @@ java-classes: $(JAVA_SOURCE)

 clean:
 	rm -rf $(ANDROID_OBJ_ROOT) $(ANDROID_LIB_ROOT)
+
+distclean: clean
+	-rm -rf $(FFTW_SRC) third_party/$(FFTW_TAR)
@@ -0,0 +1,294 @@
+#
+# Makefile.windows — кросс-компиляция WDSP для 64-bit Windows на Linux (MinGW-w64)
+#
+# Использование:
+#   make -f Makefile.windows            # DLL + статическая библиотека
+#   make -f Makefile.windows dll        # только DLL
+#   make -f Makefile.windows static     # только статическая библиотека
+#   make -f Makefile.windows clean      # удалить артефакты сборки
+#   make -f Makefile.windows distclean  # + удалить собранный и скачанный FFTW
+#
+# Зависимости (установить до запуска):
+#   Arch:    sudo pacman -S mingw-w64-gcc         (тянет mingw-w64-headers автоматически)
+#   Debian:  sudo apt install mingw-w64            (мета-пакет: компилятор + заголовки + CRT)
+#   Fedora:  sudo dnf install mingw64-gcc          (тянет mingw64-headers автоматически)
+#   + wget или curl, tar
+#
+
+MINGW_PREFIX ?= x86_64-w64-mingw32
+
+CC     := $(MINGW_PREFIX)-gcc
+AR     := $(MINGW_PREFIX)-ar
+RANLIB := $(MINGW_PREFIX)-ranlib
+
+CFLAGS  ?= -O3 -Wno-parentheses -Wno-incompatible-pointer-types
+LDFLAGS ?=
+
+# ── FFTW: сборка из исходников ────────────────────────────────────────────────
+
+FFTW_VERSION := 3.3.11
+FFTW_TAR     := fftw-$(FFTW_VERSION).tar.gz
+FFTW_URL     := https://fftw.org/pub/fftw/$(FFTW_TAR)
+
+# Абсолютные пути — нужны для configure --prefix и sub-make -C
+FFTW_SRC     := $(CURDIR)/third_party/fftw-$(FFTW_VERSION)
+FFTW_INST    := $(CURDIR)/third_party/fftw-win64
+
+# Раздельные build-директории для double и float precision
+FFTW_BUILD_D := $(FFTW_INST)/build-double
+FFTW_BUILD_F := $(FFTW_INST)/build-float
+
+# FFTW_MARCH задаёт минимальный набор инструкций для Windows-бинарника.
+# haswell (2013+) гарантирует AVX2+FMA — оптимальный выбор для SDR на современном ПК.
+# Для совместимости со старыми машинами замените на -march=sandybridge (AVX, 2011+)
+# или -march=core2 (SSSE3, 2007+).
+FFTW_MARCH ?= -march=haswell
+
+FFTW_CONF_COMMON = \
+  --host=$(MINGW_PREFIX) \
+  CC=$(CC) \
+  CFLAGS="$(FFTW_MARCH)" \
+  --prefix=$(FFTW_INST) \
+  --enable-shared \
+  --disable-static \
+  --with-our-malloc \
+  --enable-sse2 \
+  --enable-avx \
+  --enable-avx2 \
+  --disable-fortran \
+  --quiet
+
+# Sentinel-файлы для отслеживания готовности FFTW
+FFTW_HEADER := $(FFTW_INST)/include/fftw3.h
+FFTW_LIB_D  := $(FFTW_INST)/lib/libfftw3.dll.a
+FFTW_LIB_F  := $(FFTW_INST)/lib/libfftw3f.dll.a
+FFTW_DLL_D  := $(FFTW_INST)/bin/libfftw3-3.dll
+FFTW_DLL_F  := $(FFTW_INST)/bin/libfftw3f-3.dll
+
+# ── Выходные директории и файлы ───────────────────────────────────────────────
+
+OBJDIR := obj_win
+OUTDIR := lib_win
+
+DLL        := $(OUTDIR)/libwdsp.dll
+IMPLIB     := $(OUTDIR)/libwdsp.dll.a
+STATIC_LIB := $(OUTDIR)/libwdsp.a
+
+# ── Исходники (копия из основного Makefile) ───────────────────────────────────
+
+SOURCES = \
+  amd.c ammod.c amsq.c analyzer.c anf.c anr.c apfshadow.c \
+  bandpass.c calcc.c calculus.c cblock.c cfcomp.c cfir.c \
+  channel.c cmath.c compress.c delay.c dexp.c div.c doublepole.c \
+  eer.c emnr.c emph.c eq.c fcurve.c FDnoiseIQ.c fir.c firmin.c \
+  fmd.c fmmod.c fmsq.c gain.c gaussian.c gen.c icfir.c iir.c \
+  impulse_cache.c iobuffs.c iqc.c linux_port.c lmath.c main.c \
+  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 wfmd.c wfmmod.c wisdom.c zetaHat.c
+
+OBJS := $(addprefix $(OBJDIR)/, $(SOURCES:.c=.o))
+
+# ── NR3/NR4: исходники third-party ────────────────────────────────────────────
+#
+# Компилируем rnnoise и libspecbleach напрямую в obj_win/, минуя их sub-make.
+# Это исключает конфликт с Linux-сборкой: third_party/*.a содержат ELF-объекты,
+# и если они "свежее" исходников, sub-make считает их up-to-date и не пересобирает.
+
+RNNOISE_DIR    := third_party/rnnoise
+SPECBLEACH_DIR := third_party/libspecbleach
+
+RNNOISE_SRCS := \
+  src/denoise.c src/celt_lpc.c src/kiss_fft.c src/nnet.c \
+  src/nnet_default.c src/parse_lpcnet_weights.c src/pitch.c \
+  src/rnn.c src/rnnoise_data.c \
+  src/rnnoise_data_1.c src/rnnoise_data_2.c src/rnnoise_data_3.c \
+  src/rnnoise_data_4.c src/rnnoise_data_5.c src/rnnoise_data_6.c \
+  src/rnnoise_tables.c
+
+SPECBLEACH_SRCS := \
+  src/processors/specbleach_adenoiser.c \
+  src/processors/specbleach_denoiser.c \
+  src/processors/adaptivedenoiser/adaptive_denoiser.c \
+  src/shared/gain_estimation/gain_estimators.c \
+  src/shared/noise_estimation/adaptive_noise_estimator.c \
+  src/shared/pre_estimation/absolute_hearing_thresholds.c \
+  src/shared/pre_estimation/critical_bands.c \
+  src/shared/pre_estimation/masking_estimator.c \
+  src/shared/pre_estimation/noise_scaling_criterias.c \
+  src/shared/pre_estimation/spectral_smoother.c \
+  src/shared/pre_estimation/transient_detector.c \
+  src/shared/post_estimation/noise_floor_manager.c \
+  src/shared/post_estimation/postfilter.c \
+  src/shared/post_estimation/spectral_whitening.c \
+  src/shared/utils/denoise_mixer.c \
+  src/shared/utils/general_utils.c \
+  src/shared/utils/spectral_features.c \
+  src/shared/utils/spectral_utils.c \
+  src/shared/stft/stft_processor.c \
+  src/shared/stft/fft_transform.c \
+  src/shared/stft/stft_buffer.c \
+  src/shared/stft/stft_windows.c
+
+RNNOISE_OBJS    := $(addprefix $(OBJDIR)/rnnoise/, $(RNNOISE_SRCS:.c=.o))
+SPECBLEACH_OBJS := $(addprefix $(OBJDIR)/specbleach/, $(SPECBLEACH_SRCS:.c=.o))
+
+RNNOISE_LIB    := $(OUTDIR)/librnnoise.a
+SPECBLEACH_LIB := $(OUTDIR)/libspecbleach.a
+
+# ── Флаги компиляции ──────────────────────────────────────────────────────────
+
+NR34_INC  := -I$(RNNOISE_DIR)/include -I$(SPECBLEACH_DIR)/include
+FFTW_INC  := -I$(FFTW_INST)/include
+
+COMPILE   := $(CC) $(CFLAGS) $(FFTW_INC) $(NR34_INC)
+
+NR34_DEPS := $(RNNOISE_LIB) $(SPECBLEACH_LIB)
+
+# avrt нужен т.к. comm.h включает <avrt.h> под _WIN32
+LINK_LIBS := $(NR34_DEPS) \
+             -L$(FFTW_INST)/lib -lfftw3 -lfftw3f \
+             -lavrt -lm
+
+# ── Цели ──────────────────────────────────────────────────────────────────────
+
+.PHONY: all dll static fftw fftw-double fftw-float third-party check-tools clean distclean
+
+all: dll static
+
+dll:    check-tools fftw third-party $(DLL)
+static: check-tools fftw third-party $(STATIC_LIB)
+
+# ── Проверка инструментов ─────────────────────────────────────────────────────
+
+check-tools:
+	@command -v $(CC) >/dev/null 2>&1 || { \
+	  echo ""; \
+	  echo "ОШИБКА: $(CC) не найден. Установите mingw-w64:"; \
+	  echo "  Arch:    sudo pacman -S mingw-w64-gcc"; \
+	  echo "  Debian:  sudo apt install mingw-w64"; \
+	  echo "  Fedora:  sudo dnf install mingw64-gcc"; \
+	  echo ""; exit 1; }
+	@echo '#include <windows.h>' | $(CC) -x c - -fsyntax-only -Wno-pragma-once-outside-header 2>/dev/null || { \
+	  echo ""; \
+	  echo "ОШИБКА: Windows API заголовки не найдены."; \
+	  echo "Установите пакет с заголовками MinGW-w64:"; \
+	  echo "  Arch:    sudo pacman -S mingw-w64-headers"; \
+	  echo "  Debian:  sudo apt install mingw-w64          # не gcc-mingw-w64-x86-64!"; \
+	  echo "  Fedora:  sudo dnf install mingw64-headers"; \
+	  echo ""; exit 1; }
+	@(command -v wget >/dev/null 2>&1 || command -v curl >/dev/null 2>&1) || { \
+	  echo "ОШИБКА: требуется wget или curl"; exit 1; }
+
+# ── FFTW: скачать исходники → собрать double → собрать float ──────────────────
+
+fftw: fftw-double fftw-float
+
+# Шаг 1: скачать tarball
+third_party/$(FFTW_TAR):
+	mkdir -p third_party
+	@echo ">>> Скачиваем FFTW $(FFTW_VERSION)..."
+	@if command -v wget >/dev/null 2>&1; then \
+	  wget -q --show-progress -O $@ "$(FFTW_URL)"; \
+	else \
+	  curl -L --progress-bar -o $@ "$(FFTW_URL)"; \
+	fi
+
+# Шаг 2: распаковать
+$(FFTW_SRC)/configure: third_party/$(FFTW_TAR)
+	@echo ">>> Распаковываем FFTW..."
+	tar xf $< -C third_party/
+	@touch $@
+
+# Шаг 3a: configure для double precision
+$(FFTW_BUILD_D)/Makefile: $(FFTW_SRC)/configure
+	mkdir -p $(FFTW_BUILD_D)
+	@echo ">>> Конфигурируем FFTW (double)..."
+	cd $(FFTW_BUILD_D) && $(FFTW_SRC)/configure $(FFTW_CONF_COMMON)
+
+# Шаг 4a: сборка и установка double precision
+fftw-double: $(FFTW_LIB_D)
+$(FFTW_LIB_D): $(FFTW_BUILD_D)/Makefile
+	@echo ">>> Собираем FFTW (double)..."
+	$(MAKE) -C $(FFTW_BUILD_D) install
+
+# Шаг 3b: configure для float precision (--enable-float)
+$(FFTW_BUILD_F)/Makefile: $(FFTW_SRC)/configure
+	mkdir -p $(FFTW_BUILD_F)
+	@echo ">>> Конфигурируем FFTW (float)..."
+	cd $(FFTW_BUILD_F) && $(FFTW_SRC)/configure $(FFTW_CONF_COMMON) --enable-float
+
+# Шаг 4b: сборка и установка float precision
+# Зависит от fftw-double чтобы install-шаги не конкурировали за prefix
+fftw-float: $(FFTW_LIB_F)
+$(FFTW_LIB_F): $(FFTW_BUILD_F)/Makefile $(FFTW_LIB_D)
+	@echo ">>> Собираем FFTW (float)..."
+	$(MAKE) -C $(FFTW_BUILD_F) install
+
+# ── Third-party библиотеки (кросс-компиляция) ─────────────────────────────────
+
+third-party: $(NR34_DEPS)
+
+$(RNNOISE_LIB): $(RNNOISE_OBJS) | $(OUTDIR)
+	$(AR) rv $@ $^
+	$(RANLIB) $@
+
+$(SPECBLEACH_LIB): $(SPECBLEACH_OBJS) | $(OUTDIR)
+	$(AR) rv $@ $^
+	$(RANLIB) $@
+
+# % в GNU make совпадает через '/', поэтому паттерн покрывает вложенные пути
+$(OBJDIR)/rnnoise/%.o: $(RNNOISE_DIR)/%.c
+	@mkdir -p $(@D)
+	$(CC) $(CFLAGS) -I$(RNNOISE_DIR)/include -I$(RNNOISE_DIR)/src -c -o $@ $<
+
+$(OBJDIR)/specbleach/%.o: $(SPECBLEACH_DIR)/%.c | $(FFTW_LIB_D)
+	@mkdir -p $(@D)
+	$(CC) $(CFLAGS) -I$(FFTW_INST)/include \
+	  -I$(SPECBLEACH_DIR)/include -I$(SPECBLEACH_DIR)/src -I$(SPECBLEACH_DIR)/src/shared \
+	  -c -o $@ $<
+
+# ── Основная библиотека ───────────────────────────────────────────────────────
+
+$(OUTDIR) $(OBJDIR):
+	mkdir -p $@
+
+$(DLL): $(OBJS) $(NR34_DEPS) $(FFTW_LIB_D) $(FFTW_LIB_F) | $(OUTDIR)
+	@echo ">>> Линкуем $@..."
+	$(CC) -shared \
+	  -Wl,--out-implib,$(IMPLIB) \
+	  $(LDFLAGS) \
+	  -o $@ \
+	  $(OBJS) $(LINK_LIBS)
+	@cp $(FFTW_DLL_D) $(FFTW_DLL_F) $(OUTDIR)/
+	@echo ""
+	@echo "=== Готово ==="
+	@echo "  DLL:        $(DLL)"
+	@echo "  Import lib: $(IMPLIB)"
+	@echo "  Заголовок:  wdsp.h"
+	@echo ""
+	@echo "  FFTW DLL скопированы в $(OUTDIR)/ — включите их в дистрибутив"
+	@echo "  рядом с libwdsp.dll."
+
+$(STATIC_LIB): $(OBJS) | $(OUTDIR)
+	$(AR) rv $@ $(OBJS)
+	$(RANLIB) $@
+	@echo ""
+	@echo "=== Готово ==="
+	@echo "  Статическая: $(STATIC_LIB)"
+	@echo "  Линкуйте приложение с:"
+	@echo "    $(STATIC_LIB) $(NR34_DEPS)"
+	@echo "    -L$(FFTW_INST)/lib -lfftw3 -lfftw3f -lavrt -lm"
+
+$(OBJDIR)/%.o: %.c | $(OBJDIR) $(FFTW_LIB_D)
+	$(COMPILE) -c -o $@ $<
+
+# ── Очистка ───────────────────────────────────────────────────────────────────
+
+clean:
+	-rm -rf $(OBJDIR) $(OUTDIR)
+
+# distclean удаляет собранный и скачанный FFTW — при следующем запуске пересоберётся
+distclean: clean
+	-rm -rf $(FFTW_INST) $(FFTW_SRC) third_party/$(FFTW_TAR)
@@ -13,14 +13,18 @@ wdsp/
 ├── org_openhpsdr_dsp_Wdsp.c/.h  — JNI bridge
 ├── Makefile                     — host build (Linux / macOS)
 ├── Makefile.android             — Android cross-build
+├── Makefile.windows             — Windows cross-build (MinGW-w64, runs on Linux)
 ├── java/
 │   └── org/openhpsdr/dsp/Wdsp.java
 ├── third_party/
-│   ├── fftw/              — FFTW source (download manually, see below)
+│   ├── fftw-3.3.11/       — FFTW source (downloaded automatically)
+│   ├── fftw-win64/        — FFTW Windows binaries (downloaded automatically)
 │   ├── rnnoise/           — RNNoise noise suppression
 │   └── libspecbleach/     — spectral noise reduction
 ├── obj/                   — object files (generated)
-└── lib/                   — built libraries (generated)
+├── obj_win/               — Windows object files (generated)
+├── lib/                   — built libraries (generated)
+└── lib_win/               — Windows libraries (generated)
 ```

 ---
@@ -81,20 +85,95 @@ make clean

 ---

+## Windows build (кросс-компиляция с Linux)
+
+Сборка выполняется на Linux с помощью MinGW-w64. FFTW скачивается и собирается из исходников автоматически.
+
+### Зависимости
+
+```bash
+# Arch
+sudo pacman -S mingw-w64-gcc        # тянет mingw-w64-headers автоматически
+
+# Debian / Ubuntu — нужен мета-пакет mingw-w64, а не только компилятор
+sudo apt install mingw-w64
+
+# Fedora
+sudo dnf install mingw64-gcc        # тянет mingw64-headers автоматически
+```
+
+Также нужен `wget` или `curl`.
+
+> **Debian/Ubuntu:** пакет `gcc-mingw-w64-x86-64` содержит только компилятор без Windows API заголовков — `Windows.h` не будет найден. Используйте `mingw-w64`.
+
+### Сборка
+
+```bash
+# DLL + статическая библиотека (рекомендуется)
+make -f Makefile.windows
+
+# Только DLL
+make -f Makefile.windows dll
+
+# Только статическая библиотека
+make -f Makefile.windows static
+```
+
+Результат в `lib_win/`:
+
+| Файл | Назначение |
+|---|---|
+| `libwdsp.dll` | DLL для Windows |
+| `libwdsp.dll.a` | Import library для MinGW (`-lwdsp`) |
+| `libwdsp.a` | Статическая библиотека |
+| `libfftw3-3.dll` | Runtime dependency — распространять вместе с DLL |
+| `libfftw3f-3.dll` | Runtime dependency — распространять вместе с DLL |
+
+### Линковка приложения
+
+**MinGW:**
+```bash
+gcc myapp.c -Llib_win -lwdsp -o myapp.exe
+```
+
+**MSVC** — переименовать `libwdsp.dll.a` в `libwdsp.lib` и подключить обычным образом.
+
+**Статически** (MinGW):
+```bash
+gcc myapp.c libwdsp.a third_party/rnnoise/librnnoise.a \
+    third_party/libspecbleach/libspecbleach.a \
+    -Lthird_party/fftw-win64 -lfftw3-3 -lfftw3f-3 -lavrt -lm -o myapp.exe
+```
+
+### Параметры
+
+| Переменная | По умолчанию | Описание |
+|---|---|---|
+| `MINGW_PREFIX` | `x86_64-w64-mingw32` | Префикс MinGW toolchain |
+| `FFTW_VERSION` | `3.3.11` | Версия FFTW |
+| `CFLAGS` | `-O3 -Wno-parentheses` | Флаги компилятора |
+
+```bash
+make -f Makefile.windows MINGW_PREFIX=i686-w64-mingw32   # 32-bit Windows
+```
+
+### Очистка
+
+```bash
+make -f Makefile.windows clean       # артефакты сборки
+make -f Makefile.windows distclean   # + удалить скачанный FFTW
+```
+
+---
+
 ## Android build

 ### Prerequisites

 1. **Android NDK** r23 or newer
-2. **FFTW source** — download and extract into `third_party/fftw/`:
+2. **Java compiler** — for building the `.class` file (Android Studio's JBR or any JDK)

-```bash
-wget https://www.fftw.org/fftw-3.3.10.tar.gz
-tar xf fftw-3.3.10.tar.gz
-mv fftw-3.3.10 third_party/fftw
-```
-
-3. **Java compiler** — for building the `.class` file (Android Studio's JBR or any JDK)
+FFTW 3.3.11 скачивается и собирается автоматически. Чтобы использовать уже скачанный исходник, передайте `FFTW_SRC=/path/to/fftw-3.3.11`.

 ### Build

@@ -144,8 +223,8 @@ make -f Makefile.android \
 ### Clean

 ```bash
-make -f Makefile.android clean
-# removes obj/android/ and lib/android/
+make -f Makefile.android clean       # удаляет obj/android/ и lib/android/
+make -f Makefile.android distclean   # + удаляет скачанный FFTW
 ```

 ---
@@ -233,6 +233,23 @@ void create_rxa (int channel)
 		max(2048, ch[channel].dsp_size),				// number of coefficients for noise filter
 		0);												// minimum phase flag

+	// WFM demod
+	rxa[channel].wfmd.p = create_wfmd (
+		0,												// run
+		ch[channel].dsp_size,							// buffer size
+		rxa[channel].midbuff,							// pointer to input buffer
+		rxa[channel].midbuff,							// pointer to output buffer
+		ch[channel].dsp_rate,							// sample rate
+		75000.0,										// deviation
+		20.0,											// f_low
+		15000.0,										// f_high
+		0.02,											// tau - for dc removal
+		1,												// run de-emphasis
+		75.0e-6,										// de-emphasis time constant
+		0.5,											// audio gain
+		max(2048, ch[channel].dsp_size),				// # coefs for audio cutoff filter
+		0);												// min phase flag for audio cutoff filter
+
 	// snba
 	rxa[channel].snba.p = create_snba (
 		0,												// run
@@ -579,6 +596,7 @@ void destroy_rxa (int channel)
 	destroy_anf (rxa[channel].anf.p);
 	destroy_eqp (rxa[channel].eqp.p);
 	destroy_snba (rxa[channel].snba.p);
+	destroy_wfmd (rxa[channel].wfmd.p);
 	destroy_fmsq (rxa[channel].fmsq.p);
 	destroy_fmd (rxa[channel].fmd.p);
 	destroy_amd (rxa[channel].amd.p);
@@ -614,6 +632,7 @@ void flush_rxa (int channel)
 	flush_amd (rxa[channel].amd.p);
 	flush_fmd (rxa[channel].fmd.p);
 	flush_fmsq (rxa[channel].fmsq.p);
+	flush_wfmd (rxa[channel].wfmd.p);
 	flush_snba (rxa[channel].snba.p);
 	flush_eqp (rxa[channel].eqp.p);
 	flush_anf (rxa[channel].anf.p);
@@ -649,6 +668,7 @@ void xrxa (int channel)
 	xamd (rxa[channel].amd.p);
 	xfmd (rxa[channel].fmd.p);
 	xfmsq (rxa[channel].fmsq.p);
+	xwfmd (rxa[channel].wfmd.p);
 	xbpsnbain (rxa[channel].bpsnba.p, 1);
 	xbpsnbaout (rxa[channel].bpsnba.p, 1);
 	xsnba (rxa[channel].snba.p);
@@ -733,6 +753,7 @@ void setDSPSamplerate_rxa (int channel)
 	setSamplerate_fmd (rxa[channel].fmd.p, ch[channel].dsp_rate);
 	setBuffers_fmsq (rxa[channel].fmsq.p, rxa[channel].midbuff, rxa[channel].midbuff, rxa[channel].fmd.p->audio);
 	setSamplerate_fmsq (rxa[channel].fmsq.p, ch[channel].dsp_rate);
+	setSamplerate_wfmd (rxa[channel].wfmd.p, ch[channel].dsp_rate);
 	setSamplerate_snba (rxa[channel].snba.p, ch[channel].dsp_rate);
 	setSamplerate_eqp (rxa[channel].eqp.p, ch[channel].dsp_rate);
 	setSamplerate_anf (rxa[channel].anf.p, ch[channel].dsp_rate);
@@ -794,6 +815,8 @@ void setDSPBuffsize_rxa (int channel)
 	setSize_fmd (rxa[channel].fmd.p, ch[channel].dsp_size);
 	setBuffers_fmsq (rxa[channel].fmsq.p, rxa[channel].midbuff, rxa[channel].midbuff, rxa[channel].fmd.p->audio);
 	setSize_fmsq (rxa[channel].fmsq.p, ch[channel].dsp_size);
+	setBuffers_wfmd (rxa[channel].wfmd.p, rxa[channel].midbuff, rxa[channel].midbuff);
+	setSize_wfmd (rxa[channel].wfmd.p, ch[channel].dsp_size);
 	setBuffers_snba (rxa[channel].snba.p, rxa[channel].midbuff, rxa[channel].midbuff);
 	setSize_snba (rxa[channel].snba.p, ch[channel].dsp_size);
 	setBuffers_eqp (rxa[channel].eqp.p, rxa[channel].midbuff, rxa[channel].midbuff);
@@ -857,6 +880,7 @@ void SetRXAMode (int channel, int mode)
 		rxa[channel].mode = mode;
 		rxa[channel].amd.p->run	 = 0;
 		rxa[channel].fmd.p->run	 = 0;
+		rxa[channel].wfmd.p->run = 0;
 		rxa[channel].agc.p->run	 = 1;
 		switch (mode)
 		{
@@ -875,6 +899,10 @@ void SetRXAMode (int channel, int mode)
 			rxa[channel].fmd.p->run	 = 1;
 			rxa[channel].agc.p->run	 = 0;
 			break;
+		case RXA_WFM:
+			rxa[channel].wfmd.p->run = 1;
+			rxa[channel].agc.p->run	 = 0;
+			break;
 		default:

 			break;
@@ -961,13 +989,14 @@ void RXAbpsnbaCheck (int channel, int mode, int notch_run)
 			run_notches = 0;
 			break;
 		case RXA_FM:
+		case RXA_WFM:
 			f_low  = +a->abs_low_freq;
 			f_high = +a->abs_high_freq;
 			run_notches = 0;
 			break;
 		case RXA_DRM:
 		case RXA_SPEC:
-		
+
 			break;
 	}
 	// 'run' and 'position' are examined at run time; no filter changes required.
@@ -1010,6 +1039,7 @@ void RXAbpsnbaSet (int channel)
 			a->position = 1;
 			break;
 		case RXA_FM:
+		case RXA_WFM:
 			a->run = rxa[channel].snba.p->run;
 			a->position = 1;
 			break;
@@ -1046,6 +1076,7 @@ void RXASetNC (int channel, int nc)
 	SetRXAFMSQNC				(channel, nc);
 	SetRXAFMNCde				(channel, nc);
 	SetRXAFMNCaud				(channel, nc);
+	SetRXAWFMNCaud				(channel, nc);
 	SetChannelState (channel, oldstate, 0);
 }

@@ -1059,4 +1090,5 @@ void RXASetMP (int channel, int mp)
 	SetRXAFMSQMP				(channel, mp);
 	SetRXAFMMPde				(channel, mp);
 	SetRXAFMMPaud				(channel, mp);
+	SetRXAWFMMPaud				(channel, mp);
 }
@@ -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;
@@ -357,7 +357,22 @@ void create_txa (int channel)
 		1,											// run bandpass filter
 		max(2048, ch[channel].dsp_size),			// number coefficients for bandpass filter
 		0);											// minimum phase flag
-	
+
+	txa[channel].wfmmod.p = create_wfmmod (
+		0,											// run - OFF by default
+		ch[channel].dsp_size,						// size
+		txa[channel].midbuff,						// pointer to input buffer
+		txa[channel].midbuff,						// pointer to output buffer
+		ch[channel].dsp_rate,						// samplerate
+		75000.0,									// deviation
+		20.0,										// low cutoff frequency
+		15000.0,									// high cutoff frequency
+		1,											// run pre-emphasis
+		75.0e-6,									// pre-emphasis time constant
+		1,											// run bandpass filter
+		max(2048, ch[channel].dsp_size),			// number coefficients for bandpass filter
+		0);											// minimum phase flag
+
 	txa[channel].gen1.p = create_gen (
 		0,											// run
 		ch[channel].dsp_size,						// buffer size
@@ -490,6 +505,7 @@ void destroy_txa (int channel)
 	destroy_meter (txa[channel].alcmeter.p);
 	destroy_uslew (txa[channel].uslew.p);
 	destroy_gen (txa[channel].gen1.p);
+	destroy_wfmmod (txa[channel].wfmmod.p);
 	destroy_fmmod (txa[channel].fmmod.p);
 	destroy_ammod (txa[channel].ammod.p);
 	destroy_wcpagc (txa[channel].alc.p);
@@ -544,6 +560,7 @@ void flush_txa (int channel)
 	flush_wcpagc (txa[channel].alc.p);
 	flush_ammod (txa[channel].ammod.p);
 	flush_fmmod (txa[channel].fmmod.p);
+	flush_wfmmod (txa[channel].wfmmod.p);
 	flush_gen (txa[channel].gen1.p);
 	flush_uslew (txa[channel].uslew.p);
 	flush_meter (txa[channel].alcmeter.p);
@@ -580,6 +597,7 @@ void xtxa (int channel)
 	xammod (txa[channel].ammod.p);					// AM Modulator
 	xemphp (txa[channel].preemph.p, 1);				// FM pre-emphasis (second option)
 	xfmmod (txa[channel].fmmod.p);					// FM Modulator
+	xwfmmod (txa[channel].wfmmod.p);				// WFM Modulator (pre-emphasis is internal)
 	xgen (txa[channel].gen1.p);						// output signal generator (TUN and Two-tone)
 	xuslew (txa[channel].uslew.p);					// up-slew for AM, FM, and gens
 	xmeter (txa[channel].alcmeter.p);				// ALC Meter
@@ -653,6 +671,7 @@ void setDSPSamplerate_txa (int channel)
 	setSamplerate_wcpagc (txa[channel].alc.p, ch[channel].dsp_rate);
 	setSamplerate_ammod (txa[channel].ammod.p, ch[channel].dsp_rate);
 	setSamplerate_fmmod (txa[channel].fmmod.p, ch[channel].dsp_rate);
+	setSamplerate_wfmmod (txa[channel].wfmmod.p, ch[channel].dsp_rate);
 	setSamplerate_gen (txa[channel].gen1.p, ch[channel].dsp_rate);
 	setSamplerate_uslew (txa[channel].uslew.p, ch[channel].dsp_rate);
 	setSamplerate_meter (txa[channel].alcmeter.p, ch[channel].dsp_rate);
@@ -723,6 +742,8 @@ void setDSPBuffsize_txa (int channel)
 	setSize_ammod (txa[channel].ammod.p, ch[channel].dsp_size);
 	setBuffers_fmmod (txa[channel].fmmod.p, txa[channel].midbuff, txa[channel].midbuff);
 	setSize_fmmod (txa[channel].fmmod.p, ch[channel].dsp_size);
+	setBuffers_wfmmod (txa[channel].wfmmod.p, txa[channel].midbuff, txa[channel].midbuff);
+	setSize_wfmmod (txa[channel].wfmmod.p, ch[channel].dsp_size);
 	setBuffers_gen (txa[channel].gen1.p, txa[channel].midbuff, txa[channel].midbuff);
 	setSize_gen (txa[channel].gen1.p, ch[channel].dsp_size);
 	setBuffers_uslew (txa[channel].uslew.p, txa[channel].midbuff, txa[channel].midbuff);
@@ -758,6 +779,7 @@ void SetTXAMode (int channel, int mode)
 		txa[channel].mode = mode;
 		txa[channel].ammod.p->run	= 0;
 		txa[channel].fmmod.p->run	= 0;
+		txa[channel].wfmmod.p->run	= 0;
 		txa[channel].preemph.p->run = 0;
 		switch (mode)
 		{
@@ -779,6 +801,10 @@ void SetTXAMode (int channel, int mode)
 			txa[channel].fmmod.p->run	= 1;
 			txa[channel].preemph.p->run = 1;
 			break;
+		case TXA_WFM:
+			// wfmmod carries its own RC pre-emphasis; the shared emphp stays off
+			txa[channel].wfmmod.p->run	= 1;
+			break;
 		default:

 			break;
@@ -818,10 +844,11 @@ void TXAResCheck (int channel)

 int TXAUslewCheck (int channel)
 {
-	return	(txa[channel].ammod.p->run == 1) ||
-			(txa[channel].fmmod.p->run == 1) ||
-			(txa[channel].gen0.p->run  == 1) ||
-			(txa[channel].gen1.p->run  == 1);
+	return	(txa[channel].ammod.p->run  == 1) ||
+			(txa[channel].fmmod.p->run  == 1) ||
+			(txa[channel].wfmmod.p->run == 1) ||
+			(txa[channel].gen0.p->run   == 1) ||
+			(txa[channel].gen1.p->run   == 1);
 }

 void TXASetupBPFilters (int channel)
@@ -855,6 +882,7 @@ void TXASetupBPFilters (int channel)
 	case TXA_AM:
 	case TXA_SAM:
 	case TXA_FM:
+	case TXA_WFM:
 		if (txa[channel].compressor.p->run)
 		{
 			CalcBandpassFilter (txa[channel].bp0.p, 0.0, txa[channel].f_high, 2.0);
@@ -914,6 +942,7 @@ void TXASetNC (int channel, int nc)
 	SetTXAFMEmphNC				(channel, nc);
 	SetTXAEQNC					(channel, nc);
 	SetTXAFMNC					(channel, nc);
+	SetTXAWFMNC					(channel, nc);
 	SetTXACFIRNC				(channel, nc);
 	SetChannelState (channel, oldstate, 0);
 }
@@ -925,6 +954,7 @@ void TXASetMP (int channel, int mp)
 	SetTXAFMEmphMP				(channel, mp);
 	SetTXAEQMP					(channel, mp);
 	SetTXAFMMP					(channel, mp);
+	SetTXAWFMMP					(channel, mp);
 }

 PORT
@@ -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++)
+		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;
+
+		for (i = 0; i < buff_size; i++)
 		{
-			a->d[a->in_idx] = a->in_buff[2 * i + 0];
+			double dsamp;
+			int base, n1;

-			y = 0;
-			sigma = 0;
+			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);

-			for (j = 0; j < a->n_taps; j++)
-			{
-				idx = (a->in_idx + j + a->delay) & a->mask;
-				y += a->w[j] * a->d[idx];
-				sigma += a->d[idx] * a->d[idx];
-			}
 			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++)
+		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;
+
+		for (i = 0; i < buff_size; i++)
 		{
-			a->d[a->in_idx] = a->in_buff[2 * i + 0];
+			double dsamp;
+			int base, n1;

-			y = 0;
-			sigma = 0;
+			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);

-			for (j = 0; j < a->n_taps; j++)
-			{
-				idx = (a->in_idx + j + a->delay) & a->mask;
-				y += a->w[j] * a->d[idx];
-				sigma += a->d[idx] * a->d[idx];
-			}
 			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,29 +374,40 @@ 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];
-					else
-						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
+				{
+					for (j = 0; j < m1; j++)		oa[k + j] += sv[j];
+					for (; j < incr; j++)			oa[j - m1] += sv[j];
 				}
 			}
-			a->saveidx = (a->saveidx + 1) % a->ovrlp;
-			a->oainidx = (a->oainidx + a->incr) % a->oasize;
+			if (++a->saveidx == ovrlp) a->saveidx = 0;
+			if ((a->oainidx += incr) >= oasize) a->oainidx -= oasize;
 		}
-		for (i = 0; i < a->bsize; i++)
+		for (i = 0, k = a->oaoutidx; i < bsize; i++)
 		{
-			a->out[2 * i + 0] = a->outaccum[a->oaoutidx];
+			a->out[2 * i + 0] = a->outaccum[k];
 			a->out[2 * i + 1] = 0.0;
-			a->oaoutidx = (a->oaoutidx + 1) % a->oasize;
+			if (++k == oasize) k = 0;
 		}
+		a->oaoutidx = k;
 	}
 	else if (a->out != a->in)
 		memcpy (a->out, a->in, a->bsize * sizeof (complex));
@@ -33,7 +33,7 @@ warren@wpratt.com
 #endif

 #ifdef _WIN32
-#include <Windows.h>
+#include <windows.h>
 #include <process.h>
 #include <intrin.h>
 #endif
@@ -43,6 +43,11 @@ warren@wpratt.com
 #ifdef _WIN32
 #include <avrt.h>
 #endif
+
+#ifndef WDSP_FPE_GUARD
+#define WDSP_FPE_GUARD   ((void)0)
+#define WDSP_FPE_RESTORE ((void)0)
+#endif
 #include "fftw3.h"

 #include "amd.h"
@@ -70,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"
@@ -145,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,29 +1252,40 @@ 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];
-					else
-						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
+				{
+					for (j = 0; j < m1; j++)		oa[k + j] += sv[j];
+					for (; j < incr; j++)			oa[j - m1] += sv[j];
 				}
 			}
-			a->saveidx = (a->saveidx + 1) % a->ovrlp;
-			a->oainidx = (a->oainidx + a->incr) % a->oasize;
+			if (++a->saveidx == ovrlp) a->saveidx = 0;
+			if ((a->oainidx += incr) >= oasize) a->oainidx -= oasize;
 		}
-		for (i = 0; i < a->bsize; i++)
+		for (i = 0, k = a->oaoutidx; i < bsize; i++)
 		{
-			a->out[2 * i + 0] = a->outaccum[a->oaoutidx];
+			a->out[2 * i + 0] = a->outaccum[k];
 			a->out[2 * i + 1] = 0.0;
-			a->oaoutidx = (a->oaoutidx + 1) % a->oasize;
+			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