vladimir/wdsp
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Compare commits

base: vladimir/wdsp:ec93fbfcbb5446406ca1ccec84458e8cbaf194e9
Branches Tags
vladimir/wdsp:main vladimir/wdsp:cpu-optimization vladimir/wdsp:wfm
...
compare: vladimir/wdsp:main
Branches Tags
vladimir/wdsp:cpu-optimization vladimir/wdsp:main vladimir/wdsp:wfm

26 Commits
ec93fbfcbb ... main

Author SHA1 Message Date
vladimir 544e75798c Merge branch 'wfm' into cpu-optimization 2026-07-18 10:51:50 +03:00
vladimir 3ad40a00da calcc: port eu2av robustness patches for PureSignal calibration 
Three fixes from Thetis-Enhanced (Yurij eu2av), measured on Orion MK2 /
Anvelina PRO3 hardware:

- drop overrange samples (env_TX*hw_scale > 1.0) before the cubic
  xbuilder fit - they distort the fit and produce a wrong rx_scale
- optional median-ratio + MAD outlier rejection before the fit,
  controlled by new SetPSOutlierSigma (0 = off)
- fallback rx_scale estimate from the top amplitude intervals when the
  xbuilder fit fails or is rejected by rxscheck

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
 2026-07-18 10:47:44 +03:00
Uladzimir Karpenka c1964817dd wfm: compile wfmd/wfmmod in the Windows build 
Makefile.windows arrived from main, which predates the wfm module, so its
SOURCES list omitted wfmd.c/wfmmod.c even though RXA.c/TXA.c reference their
symbols -- the MinGW link failed on the missing objects. Add both to SOURCES,
matching Makefile and Makefile.android.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
 2026-07-15 11:03:03 +03:00
Uladzimir Karpenka 81f7750def Merge remote-tracking branch 'origin/main' into cpu-optimization 2026-07-15 10:50:34 +03:00
Uladzimir Karpenka d24b5ae316 Merge remote-tracking branch 'origin/main' into wfm 2026-07-15 10:50:15 +03:00
Uladzimir Karpenka 65cb3c386e wfm: add a wideband FM modulator and demodulator 
Model the pair on fmd.c / fmmod.c, but with the parts that a 75 kHz
deviation forces:

  - the demodulator discriminates with arg(x[n] * conj(x[n-1])) rather
    than a PLL; an omegaN in the tens of kHz cannot track 75 kHz.
  - emphasis is a one-pole RC (tau 75 us) on both ends, not fmd's 1/f
    fc_impulse FIR.  A 1/f FIR from f_low = 20 Hz would sit ~+57 dB at
    20 Hz, where broadcast FM specifies flat below the corner.
  - TXA_WFM leaves the shared preemph block off; wfmmod carries its own.
  - wfmmod clamps bp_fc = deviation + f_high to 0.45 * samplerate, since
    +/-90 kHz exceeds Nyquist at the rates the narrowband modes use.

Scope is mono: no 19 kHz pilot, no 38 kHz stereo subcarrier, no RDS.

Both mode enums are appended to so the ABI stays stable for clients.
The JNI bindings are deliberately left alone.

Verified against a wfmmod -> wfmd loopback: the discriminator is exact,
the dc-removal one-pole tracks |H_lp(f)| * ain * sdelta to ratio 1.000,
and a 700 + 1900 Hz two-tone comes back with THD+N ~ 3e-5 % at +/-37.5
kHz deviation.  Note the demodulator aliases unless samplerate exceeds
2 * deviation * |aud|max, so 192 kHz is the practical floor.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
 2026-07-10 07:43:07 +03:00
Uladzimir Karpenka f39eea6b01 phrot, cfcomp: cascade in a register; drop idiv from cfcomp's ring walks 
xphrot() kept x0[] and y0[] in the struct, but neither carried state between
samples: x0[n] was always the previous stage's output and y0[n] this stage's.
Cascade that single value in a register and leave x1/y1 as the real filter
state. Stores through out[] may alias the struct's doubles, so hoist the
coefficients as well. Eight first-order sections then run without touching
memory except for the state.

xcfcomp() is structurally the same overlap-add loop as xemnr(), and had the
same defect: four ring indices advanced with '% size' per step, ~5100 integer
divisions per call at fsize = 2048, which a profile showed dominating the
block (1820 samples in xcfcomp against 362 in calc_mask and ~520 in the FFTs).
The indices step by one, and iasize >= fsize and oasize >= incr always hold,
so they wrap at most once per loop: walk contiguous runs and wrap between them.

Measured in situ on an Apple M1 Pro, 512-sample buffers, cost of turning the
block on, best of 5:

    phrot    17219 ns -> 7047 ns   2.44x
    cfcomp   28918 ns -> 13637 ns  2.12x

Both are bit-identical: phrot preserves the operation order, and cfcomp only
changes integer index arithmetic. The RX chain is unchanged bit-for-bit.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
 2026-07-10 00:31:32 +03:00
Uladzimir Karpenka 91df5b1f2d varsamp: build the phase table without malloc0's memset, read h forward 
calc_varsamp() gained a transpose of the coefficient table, which showed up
as ~0.3 ms on create_varsamp()/setInRate_varsamp() (1.0 ms -> 1.4 ms). It is
a one-off setup cost, not on the sample path, but it is easy to trim: every
element of hp is written, so malloc0()'s memset of ~1 MB is dead, and h is
cold straight out of fir_bandpass(), so walk it along its fast axis and let
the prefetcher work.

Output is bit-identical; xvarsamp() is unchanged at ~34 us/buffer.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
 2026-07-10 00:13:27 +03:00
Uladzimir Karpenka 12dc701604 analyzer: block the dispatcher on a semaphore, stop spawning a thread per frame 
Two pieces of pure overhead, neither of them doing any DSP.

The POSIX QueueUserWorkItem() shim spawned a thread and immediately joined
it. That is a synchronous call -- the parallelism the Windows thread pool
provides is absent here either way -- so it only bought a pthread_create()
and its stack mmap, ~15 us, for every FFT frame. Call the function directly.
(The Windows build is untouched and still gets its thread pool; a real pool
for POSIX would be a separate change, and only pays off for num_stitch > 1.)

The dispatcher thread ran `for (each ss, LO) {...} Sleep(1);`, and Sleep(1)
is usleep(1000), so it woke 1000 times a second to re-read the same flags.
That burned 0.63% of a core even with no samples arriving. Give it a
semaphore instead, signalled by the four Spectrum*() entry points when new
samples land, and by SetAnalyzer/DestroyAnalyzer after they raise
end_dispatcher so the blocking wait always has a way out. SetAnalyzer holds
SetAnalyzerSection while it waits for the dispatcher to quit, and the
dispatcher never takes that section, so signalling from under it is safe.

Measured on an Apple M1 Pro; 16384-point complex FFT, 1024-sample buffers,
~40 pixel frames/s, CPU of all threads via getrusage:

    idle, no samples at all      0.63% of a core  ->  0.00%
    under load                   3.54%            ->  2.65%

Output pixels are bit-identical. Stressed with 3 create/destroy cycles and
24 on-the-fly SetAnalyzer reconfigurations while a second thread fed samples:
no deadlock.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
 2026-07-09 23:59:20 +03:00
Uladzimir Karpenka 645cbbb2d1 varsamp: transpose the coefficients into phases, vectorize the tap loop 
hshift() rebuilds the whole interpolated tap set on every output sample,
reading h[hidx + m*R] for m = 0..rsize-1. With R = 1024 (what rmatch asks
for) that strides 8 KB at a time through a 1.1 MB table, so every one of
the 2*rsize reads is its own cache line. It cost more than the filter it
was feeding: 258 ns per output sample against 146 ns for the tap loop.

Store the coefficients transposed instead, hp[p*rsize + m] = h[p + m*R],
so the two phases hshift() interpolates between are each contiguous. R+1
phases are needed since it reads hidx and hidx+1, and h_offset is kept in
[0,1) by the caller so hidx <= R-1. The untransposed h is freed; the
impulse cache hands back a copy, so varsamp owns it. Net memory is
unchanged.

The tap loop had the same wrap test per tap as resample.c did, so split it
at the wrap and carry four independent accumulator pairs; the ring is split
into I/Q so the taps load unit-stride.

Note a->hs is rewritten by hshift() inside the sample loop, so it must not
be hoisted behind a restrict pointer in xvarsamp().

Measured on an Apple M1 Pro, 512-sample buffers, best of 5:

    48k -> 48k   varmode=0    138580 ns -> 37370 ns   3.71x
    48k -> 48k   varmode=1    138390 ns -> 32333 ns   4.28x
    48k -> 44.1k varmode=1    141473 ns -> 38220 ns   3.70x
    44.1k -> 48k varmode=1    152307 ns -> 45263 ns   3.36x

hshift() is numerically identical -- same coefficients, different layout.
Only the reassociated tap sum rounds differently: worst deviation 7.2e-16
over four rate configurations, an SNR of 344 dB. Driven end to end through
rmatch's public API, output SNR is 306 dB and total energy matches.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
 2026-07-09 23:19:36 +03:00
Uladzimir Karpenka fd2ba84e7d emnr: linear-time aepf, drop idiv from ring walks, halve getKey's logs 
Three separate costs, found with a sampling profile of the RX chain with
emnr forced on:

aepf() averaged mask[] over a window of N = 2*psi + 1 = 41 bins by walking
the window for every one of the 2049 bins, i.e. O(msize*N). Its three spans
are all symmetric windows clipped at the array ends, so take each from a
prefix sum instead: one subtraction per output, O(msize).

xemnr() advanced four ring indices with a '% size' per step. iasize is 4096
and oasize 1024 here, and neither is known to the compiler, so each step was
a real integer division -- ~8700 of them per frame. The indices step by one
and, since iasize >= fsize and oasize >= incr always hold, wrap at most once
per loop, so walk contiguous runs and wrap between them.

calc_gain() called getKey() twice per bin with the same gamma, so the gamma
row index and its log10 were computed twice. Split getKey into keyIndex() +
keyLerp() and locate gamma once. The remaining logs go through wdsp_log10()
(new fastmath.h), accurate to 2e-13 against libm and ~2.4x its throughput;
gamma and xi are bracketed against the table limits first, so the argument
is always positive and normal. Also clamp the row index so the second
bilinear corner cannot address the next row of the 241x241 table.

Measured in situ on an Apple M1 Pro, 512-sample buffers, cost of turning
emnr on, best of 5:

    baseline                       73661 ns
    + aepf, ring walks             48927 ns   1.51x
    + getKey                       37852 ns   1.95x

Output is not bit-identical, as the prefix sum and the reassociated logs
round differently. Over 300 buffers with emnr alone the worst deviation is
4.0e-09, an SNR of 196 dB; perturbing a single input sample of the unmodified
code by one ulp diverges it from itself by 1.4e-08 (186 dB), so this change
disturbs the chain less than the last bit of the input does.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
 2026-07-09 23:08:20 +03:00
Uladzimir Karpenka 3faad14fc1 anr, anf: vectorize the LMS tap and update loops 
Both loops indexed the delay line as (in_idx + j + delay) & mask, one
masked index per tap, which made the address non-affine and stopped the
vectorizer. The window wraps at most once, so split it at the wrap and
walk two contiguous runs instead.

As in resample.c, the y/sigma reduction cannot be reassociated without
-ffast-math (which this library must not enable, see linux_port.h), so
carry four independent accumulator pairs to break the FMA dependency
chain and let the vectorizer in.

in_buff and out_buff alias in RXA -- both are midbuff -- so only the
private d/w arrays are marked restrict.

Measured in situ on an Apple M1 Pro, 512-sample buffers, cost of turning
the block on, best of 5:

    anr    59303 ns -> 20253 ns   2.93x
    anf    55073 ns -> 20511 ns   2.68x

Summation order changes, so output is not bit-identical: over 300 buffers
of the full RX chain the worst deviation is 2.4e-07, an SNR of 153 dB.
An LMS filter is an adaptive feedback loop, so its trajectory is
chaotic. As a control, perturbing a single input sample of the unmodified
code by one ulp diverges it from itself by 6.6e-07, an SNR of 144.6 dB --
i.e. this change disturbs the filter less than the last bit of the input
does.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
 2026-07-09 22:52:33 +03:00
Uladzimir Karpenka 9960a1d0e3 resample: vectorize the polyphase tap loop 
The tap loop wrapped the ring with a test on every tap:

    if ((idx_out = idx_in + j) >= ringsize) idx_out -= ringsize;

which made the address non-affine and stopped the vectorizer. Split the
walk at the wrap point instead, so both halves are unit-stride, and split
the complex ring into separate I/Q arrays so the taps load contiguously
rather than through a de-interleaving ld2.

The dot product also could not be vectorized as written: reassociating an
fp reduction needs -ffast-math, which this library must not enable (it
relies on IEEE semantics for 0/0 = NaN and x/0 = Inf, see linux_port.h).
Carry four independent accumulator pairs instead, which both breaks the
FMA dependency chain and lets the vectorizer in on any compiler.

Measured on an Apple M1 Pro, 512-sample DSP buffers, best of 5:

  xresample, decimation to 48 kHz     before      after    speedup
    192k -> 48k  (561 taps)          342.0 us    84.2 us     4.06x
    384k -> 48k  (1121 taps)         708.9 us   175.1 us     4.05x
    576k -> 48k  (1681 taps)        1074.4 us   266.9 us     4.03x
    768k -> 48k  (2241 taps)        1438.9 us   360.0 us     4.00x

  full xrxa() chain, 576k input      1148.2 us   337.6 us     3.40x
                                     10.76%       3.16%   of one core

  xresampleF (float, host audio)                          2.6x - 3.4x

Summation order changes, so the double path is not bit-identical: over
400 buffers of the full RX chain the worst deviation is 1.1e-12, an SNR
of 251 dB. The float path is bit-identical, as the cast to float absorbs
the difference. With the resampler bypassed (48k in, 48k out) the chain
is unchanged bit-for-bit.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
 2026-07-09 22:44:08 +03:00
vladimir bb7e0b3df6 add gitignore third_party 2026-06-12 21:44:25 +03:00
vladimir 4313006fa6 Auto-detect Android NDK and javac paths in Makefile.android 
Remove hardcoded /home/vladimir and /opt paths. NDK is now resolved via
ANDROID_NDK_HOME, ANDROID_SDK_ROOT/ANDROID_HOME, or the default SDK
location for Linux/macOS (latest version picked automatically). javac is
resolved via JAVA_HOME, system PATH, or Android Studio JBR on Linux/macOS.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
 2026-06-12 21:31:15 +03:00
vladimir 4f44118299 Remove committed build artifacts from third_party and fix Android FFTW cross-compile 
- Untrack all .o and .a files under third_party/ (rnnoise, libspecbleach)
- Add third_party/**/*.o and third_party/**/*.a to .gitignore
- Fix FFTW configure cross-compilation: add $(strip ...) around fftw_host
  macro call and explicit --build flag so configure detects cross-compile
  correctly (GNU make backslash-continuation in define blocks inserts a
  leading space, causing --host= arm-... to be parsed as empty --host)

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
 2026-06-12 21:27:13 +03:00
vladimir edd8991a3c Enable AVX2 for FFTW Windows build to suppress SSE2-only warning 
FFTW emits '#warning Only SSE and SSE2 are available' when cross-compiling
with a generic x86_64 target: the configure test for AVX intrinsics fails
without an explicit -mavx flag, so FFTW silently falls back to SSE2.

Fix: pass CFLAGS="-march=haswell" to configure, which enables SSE2/AVX/AVX2/FMA.
Add --enable-avx2 alongside existing --enable-avx.
FFTW_MARCH is overridable for older CPU targets:
  make -f Makefile.windows FFTW_MARCH=-march=sandybridge  # AVX only, 2011+
  make -f Makefile.windows FFTW_MARCH=-march=core2        # SSSE3, 2007+

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
 2026-06-02 13:39:24 +03:00
vladimir a3385bd1f9 Fix parallel build race: DLL link must wait for fftw-float install 
With -jN, the DLL link step could start before fftw-float finished
installing libfftw3f.dll.a, since $(DLL) only had an implicit ordering
through .o compilation (which needs FFTW_LIB_D but not FFTW_LIB_F).

Add both FFTW_LIB_D and FFTW_LIB_F as explicit prerequisites of $(DLL).

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
 2026-06-02 13:35:04 +03:00
vladimir b7d156ea40 Compile rnnoise and libspecbleach directly into obj_win/ for Windows 
Sub-make of third_party/ reported 'up to date' when Linux ELF objects
were present from a prior host build, causing the linker to receive
ELF archives instead of PE/COFF and producing undefined reference errors.

Fix: compile rnnoise and libspecbleach sources directly into obj_win/
with explicit pattern rules, producing lib_win/librnnoise.a and
lib_win/libspecbleach.a. These are fully isolated from third_party/
and never conflict with the Linux build artifacts.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
 2026-06-02 13:17:32 +03:00
vladimir d997990ceb Fix two Windows cross-compilation errors in analyzer.c 
WDSP_FPE_GUARD: the macro is defined only in linux_port.h under
#if defined(linux)||defined(__APPLE__). On the _WIN32 path it was
completely undefined. Added a fallback no-op definition in comm.h
guarded by #ifndef so it applies to Windows (and any future platform
that doesn't include linux_port.h).

volatile int* vs volatile LONG*: Win32 Interlocked functions expect
volatile LONG* (= volatile long*). The dispatcher field is volatile int.
On Windows LLP64 both are 32-bit so the operation is correct, but GCC 14
promotes this mismatch from warning to error. Suppressed with
-Wno-incompatible-pointer-types in Makefile.windows, consistent with
how MSVC handles it silently.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
 2026-06-02 12:54:56 +03:00
vladimir a77dbbb6e5 Fix Windows.h case on Linux: use lowercase windows.h 
On Linux (case-sensitive filesystem) MinGW-w64 installs the header as
windows.h (lowercase) while comm.h included <Windows.h> (capital W),
causing a fatal compile error when cross-compiling for Windows.

Changed to <windows.h> which works on both Linux/MinGW-w64 and native
Windows/MSVC (Windows header includes are case-insensitive on NTFS).

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
 2026-06-02 12:43:45 +03:00
vladimir 90c906eb3a Fix missing Windows.h: require full mingw-w64 package, not just compiler 
gcc-mingw-w64-x86-64 on Debian/Ubuntu ships only the compiler binary,
without Windows API headers — causing fatal error: Windows.h not found.
The fix is to install the mingw-w64 meta-package which includes headers,
CRT and runtime libraries.

Added a check-tools probe that compiles #include <windows.h> and prints
a clear error with the correct package name if headers are missing.
Updated README with a warning about the incomplete package.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
 2026-06-02 12:33:13 +03:00
vladimir f29970542b Pass CC explicitly to FFTW configure to suppress cross-tools warning 
Without CC=, autoconf issues "using cross tools not prefixed with host
triplet" because it cannot match the detected compiler against --host.
Passing CC=$(MINGW_PREFIX)-gcc directly resolves the ambiguity.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
 2026-06-02 12:26:39 +03:00
vladimir 69ea631daa Fix parallel build race condition in Makefile.windows 
With -jN, make started third-party and fftw targets in parallel.
libspecbleach depended on FFTW_HEADER which has no explicit rule
(it is a side effect of make install), causing 'No rule to make target'.
Object compilation also could start before FFTW headers were ready.

Fix: depend on FFTW_LIB_D (has an explicit rule) instead of FFTW_HEADER,
and add FFTW_LIB_D as order-only prerequisite for .o compilation.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
 2026-06-02 12:24:39 +03:00
vladimir 02926b36ce Auto-download FFTW source for Android build 
Makefile.android now downloads fftw-3.3.11.tar.gz automatically instead
of requiring the user to place sources in third_party/fftw/ manually.
The FFTW_SRC variable can still be overridden to use an existing copy.
Added distclean target. Updated README accordingly.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
 2026-06-02 12:15:59 +03:00
vladimir 5bbf9698a8 Build FFTW from source instead of downloading pre-built DLLs 
FFTW 3.3.11 has no pre-built Windows binaries, so Makefile.windows now
downloads the source tarball and cross-compiles it twice with MinGW-w64:
once for double precision and once for float (--enable-float), installing
both into third_party/fftw-win64/. Links with -lfftw3/-lfftw3f via
libtool-generated import libs instead of dlltool-generated ones.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
 2026-06-02 12:13:35 +03:00
71 changed files with 1998 additions and 353 deletions
Expand all files Collapse all files
.gitignore
+4
View File
@@ -7,6 +7,10 @@ java/build/
# third_party — only fftw is downloaded manually, rnnoise and libspecbleach are tracked # third_party — only fftw is downloaded manually, rnnoise and libspecbleach are tracked
third_party/fftw/ third_party/fftw/
third_party/**/*.o
third_party/**/*.a
third_party/fftw-3.3.11.tar.gz
third_party/fftw-3.3.11/
# macOS # macOS
.DS_Store .DS_Store
Makefile
+2
View File
@@ -144,6 +144,8 @@ utilities.c \
varsamp.c \ varsamp.c \
version.c \ version.c \
wcpAGC.c \ wcpAGC.c \
wfmd.c \
wfmmod.c \
wisdom.c \ wisdom.c \
zetaHat.c zetaHat.c
Makefile.android
+60 -22
View File
@@ -5,17 +5,39 @@
# Optional overrides: # Optional overrides:
# ANDROID_NDK=/path/to/ndk ANDROID_API=24 ANDROID_ABIS="arm64-v8a x86_64" # ANDROID_NDK=/path/to/ndk ANDROID_API=24 ANDROID_ABIS="arm64-v8a x86_64"
# #
# Requires FFTW source in third_party/fftw/ # FFTW source is downloaded and extracted automatically.
# Download from https://www.fftw.org/download.html and extract so that # Override FFTW_SRC to use an existing directory.
# third_party/fftw/configure exists.
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_API ?= 24
ANDROID_ABIS ?= arm64-v8a armeabi-v7a x86_64 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/') 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 # Auto-detect javac: JAVA_HOME > system PATH > Android Studio JBR (Linux/macOS)
JAVAC ?= $(JBR_BIN)/javac _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) 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 \ ANDROID_JNI_CFLAGS ?= -std=gnu89 -Wno-implicit-function-declaration -Wno-int-conversion \
-Wno-incompatible-pointer-types -Wno-incompatible-pointer-types-discards-qualifiers -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) FFTW_MAKEJOBS ?= $(shell nproc 2>/dev/null || sysctl -n hw.ncpu 2>/dev/null || echo 4)
WDSP_SOURCES = amd.c \ WDSP_SOURCES = amd.c \
@@ -94,6 +119,8 @@ utilities.c \
varsamp.c \ varsamp.c \
version.c \ version.c \
wcpAGC.c \ wcpAGC.c \
wfmd.c \
wfmmod.c \
wisdom.c \ wisdom.c \
zetaHat.c zetaHat.c
@@ -148,9 +175,9 @@ ANDROID_OBJ_ROOT := obj/android
ANDROID_LIB_ROOT := lib/android ANDROID_LIB_ROOT := lib/android
ANDROID_JAVA_OUT := lib/android/java 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),\ $(foreach abi,$(ANDROID_ABIS),\
$(ANDROID_LIB_ROOT)/$(abi)/libfftw3.so \ $(ANDROID_LIB_ROOT)/$(abi)/libfftw3.so \
$(ANDROID_LIB_ROOT)/$(abi)/libfftw3f.so \ $(ANDROID_LIB_ROOT)/$(abi)/libfftw3f.so \
@@ -158,15 +185,21 @@ all: check-fftw-src \
$(ANDROID_LIB_ROOT)/$(abi)/libwdspj.so) \ $(ANDROID_LIB_ROOT)/$(abi)/libwdspj.so) \
java-classes java-classes
check-fftw-src: # ── Автоматическое получение FFTW ────────────────────────────────────────────
@test -f $(FFTW_SRC)/configure || { \
echo ""; \ third_party/$(FFTW_TAR):
echo "ERROR: FFTW source not found at $(FFTW_SRC)/configure"; \ mkdir -p third_party
echo "Download and extract FFTW from https://www.fftw.org/download.html"; \ @echo ">>> Скачиваем FFTW $(FFTW_VERSION)..."
echo "so that $(FFTW_SRC)/configure exists."; \ @if command -v wget >/dev/null 2>&1; then \
echo ""; \ wget -q --show-progress -O $@ "$(FFTW_URL)"; \
exit 1; \ 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 # abi_target: compiler triple for NDK clang wrapper
define abi_target 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) ABI_$(1)_JNI_OBJ := $$(ABI_$(1)_OBJDIR)/$(JNI_SOURCE:.c=.o)
# --- Build fftw3 (double) --- # --- 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) @mkdir -p $$(ABI_$(1)_OBJDIR)/fftw-build $$(ABI_$(1)_FFTW_PREFIX)
cd $$(ABI_$(1)_OBJDIR)/fftw-build && \ cd $$(ABI_$(1)_OBJDIR)/fftw-build && \
$(abspath $(FFTW_SRC))/configure \ $(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)" \ CC="$$(ABI_$(1)_CC)" \
CFLAGS="-fPIC" \ CFLAGS="-fPIC" \
--prefix=$(abspath $$(ABI_$(1)_FFTW_PREFIX)) \ --prefix=$(abspath $$(ABI_$(1)_FFTW_PREFIX)) \
@@ -240,11 +274,12 @@ $$(ABI_$(1)_FFTW_STAMP):
@touch $$@ @touch $$@
# --- Build fftw3f (float/single) --- # --- 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) @mkdir -p $$(ABI_$(1)_OBJDIR)/fftwf-build $$(ABI_$(1)_FFTWF_PREFIX)
cd $$(ABI_$(1)_OBJDIR)/fftwf-build && \ cd $$(ABI_$(1)_OBJDIR)/fftwf-build && \
$(abspath $(FFTW_SRC))/configure \ $(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)" \ CC="$$(ABI_$(1)_CC)" \
CFLAGS="-fPIC" \ CFLAGS="-fPIC" \
--prefix=$(abspath $$(ABI_$(1)_FFTWF_PREFIX)) \ --prefix=$(abspath $$(ABI_$(1)_FFTWF_PREFIX)) \
@@ -311,3 +346,6 @@ java-classes: $(JAVA_SOURCE)
clean: clean:
rm -rf $(ANDROID_OBJ_ROOT) $(ANDROID_LIB_ROOT) rm -rf $(ANDROID_OBJ_ROOT) $(ANDROID_LIB_ROOT)
distclean: clean
-rm -rf $(FFTW_SRC) third_party/$(FFTW_TAR)
Makefile.windows
+169 -72
View File
@@ -6,31 +6,64 @@
# make -f Makefile.windows dll # только DLL # make -f Makefile.windows dll # только DLL
# make -f Makefile.windows static # только статическая библиотека # make -f Makefile.windows static # только статическая библиотека
# make -f Makefile.windows clean # удалить артефакты сборки # make -f Makefile.windows clean # удалить артефакты сборки
# make -f Makefile.windows distclean # + удалить скачанный FFTW # make -f Makefile.windows distclean # + удалить собранный и скачанный FFTW
# #
# Зависимости (установить до запуска): # Зависимости (установить до запуска):
# Arch: sudo pacman -S mingw-w64-gcc unzip # Arch: sudo pacman -S mingw-w64-gcc (тянет mingw-w64-headers автоматически)
# Debian: sudo apt install gcc-mingw-w64-x86-64 binutils-mingw-w64-x86-64 unzip # Debian: sudo apt install mingw-w64 (мета-пакет: компилятор + заголовки + CRT)
# Fedora: sudo dnf install mingw64-gcc mingw64-binutils unzip # Fedora: sudo dnf install mingw64-gcc (тянет mingw64-headers автоматически)
# + wget или curl # + wget или curl, tar
# #
MINGW_PREFIX ?= x86_64-w64-mingw32 MINGW_PREFIX ?= x86_64-w64-mingw32
CC := $(MINGW_PREFIX)-gcc CC := $(MINGW_PREFIX)-gcc
AR := $(MINGW_PREFIX)-ar AR := $(MINGW_PREFIX)-ar
RANLIB := $(MINGW_PREFIX)-ranlib RANLIB := $(MINGW_PREFIX)-ranlib
DLLTOOL := $(MINGW_PREFIX)-dlltool
CFLAGS ?= -O3 -Wno-parentheses CFLAGS ?= -O3 -Wno-parentheses -Wno-incompatible-pointer-types
LDFLAGS ?= LDFLAGS ?=
# ── FFTW Windows binaries ───────────────────────────────────────────────────── # ── FFTW: сборка из исходников ────────────────────────────────────────────────
FFTW_VERSION := 3.3.11 FFTW_VERSION := 3.3.11
FFTW_ZIP := fftw-$(FFTW_VERSION)-dll64.zip FFTW_TAR := fftw-$(FFTW_VERSION).tar.gz
FFTW_URL := https://fftw.org/pub/fftw/$(FFTW_ZIP) FFTW_URL := https://fftw.org/pub/fftw/$(FFTW_TAR)
FFTW_DIR := third_party/fftw-win64
# Абсолютные пути — нужны для 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
# ── Выходные директории и файлы ─────────────────────────────────────────────── # ── Выходные директории и файлы ───────────────────────────────────────────────
@@ -53,30 +86,74 @@ SOURCES = \
matchedCW.c meter.c meterlog10.c nbp.c nob.c nobII.c osctrl.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 \ 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 \ shift.c siphon.c slew.c snb.c ssql.c syncbuffs.c TXA.c \
utilities.c varsamp.c version.c wcpAGC.c wisdom.c zetaHat.c utilities.c varsamp.c version.c wcpAGC.c wfmd.c wfmmod.c wisdom.c zetaHat.c
OBJS := $(addprefix $(OBJDIR)/, $(SOURCES:.c=.o)) 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
# ── Флаги компиляции ────────────────────────────────────────────────────────── # ── Флаги компиляции ──────────────────────────────────────────────────────────
# Абсолютные пути, т.к. используются в sub-make из других директорий NR34_INC := -I$(RNNOISE_DIR)/include -I$(SPECBLEACH_DIR)/include
FFTW_ABS := $(CURDIR)/$(FFTW_DIR) FFTW_INC := -I$(FFTW_INST)/include
NR34_INC := -I third_party/rnnoise/include -I third_party/libspecbleach/include
FFTW_INC := -I$(FFTW_ABS)
COMPILE := $(CC) $(CFLAGS) $(FFTW_INC) $(NR34_INC) COMPILE := $(CC) $(CFLAGS) $(FFTW_INC) $(NR34_INC)
NR34_DEPS := third_party/rnnoise/librnnoise.a \ NR34_DEPS := $(RNNOISE_LIB) $(SPECBLEACH_LIB)
third_party/libspecbleach/libspecbleach.a
# avrt нужен т.к. comm.h включает <avrt.h> под _WIN32; fftw3f нужен libspecbleach # avrt нужен т.к. comm.h включает <avrt.h> под _WIN32
LINK_LIBS := $(NR34_DEPS) \ LINK_LIBS := $(NR34_DEPS) \
-L$(FFTW_ABS) -lfftw3-3 -lfftw3f-3 \ -L$(FFTW_INST)/lib -lfftw3 -lfftw3f \
-lavrt -lm -lavrt -lm
# ── Цели ────────────────────────────────────────────────────────────────────── # ── Цели ──────────────────────────────────────────────────────────────────────
.PHONY: all dll static fftw third-party check-tools clean distclean .PHONY: all dll static fftw fftw-double fftw-float third-party check-tools clean distclean
all: dll static all: dll static
@@ -90,79 +167,101 @@ check-tools:
echo ""; \ echo ""; \
echo "ОШИБКА: $(CC) не найден. Установите mingw-w64:"; \ echo "ОШИБКА: $(CC) не найден. Установите mingw-w64:"; \
echo " Arch: sudo pacman -S mingw-w64-gcc"; \ echo " Arch: sudo pacman -S mingw-w64-gcc"; \
echo " Debian: sudo apt install gcc-mingw-w64-x86-64"; \ echo " Debian: sudo apt install mingw-w64"; \
echo " Fedora: sudo dnf install mingw64-gcc"; \ echo " Fedora: sudo dnf install mingw64-gcc"; \
echo ""; exit 1; } echo ""; exit 1; }
@command -v $(DLLTOOL) >/dev/null 2>&1 || { \ @echo '#include <windows.h>' | $(CC) -x c - -fsyntax-only -Wno-pragma-once-outside-header 2>/dev/null || { \
echo "ОШИБКА: $(DLLTOOL) не найден"; exit 1; } 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) || { \ @(command -v wget >/dev/null 2>&1 || command -v curl >/dev/null 2>&1) || { \
echo "ОШИБКА: требуется wget или curl"; exit 1; } echo "ОШИБКА: требуется wget или curl"; exit 1; }
@command -v unzip >/dev/null 2>&1 || { \
echo "ОШИБКА: требуется unzip"; exit 1; }
# ── FFTW: скачать → распаковать → сгенерировать import-библиотеки ───────────── # ── FFTW: скачать исходники → собрать double → собрать float ──────────────────
fftw: $(FFTW_DIR)/libfftw3-3.dll.a $(FFTW_DIR)/libfftw3f-3.dll.a fftw: fftw-double fftw-float
$(FFTW_DIR)/$(FFTW_ZIP): # Шаг 1: скачать tarball
mkdir -p $(FFTW_DIR) third_party/$(FFTW_TAR):
@echo ">>> Скачиваем FFTW $(FFTW_VERSION) (Windows 64-bit)..." mkdir -p third_party
@echo ">>> Скачиваем FFTW $(FFTW_VERSION)..."
@if command -v wget >/dev/null 2>&1; then \ @if command -v wget >/dev/null 2>&1; then \
wget -q --show-progress -O $@ "$(FFTW_URL)"; \ wget -q --show-progress -O $@ "$(FFTW_URL)"; \
else \ else \
curl -L --progress-bar -o $@ "$(FFTW_URL)"; \ curl -L --progress-bar -o $@ "$(FFTW_URL)"; \
fi fi
$(FFTW_DIR)/fftw3.h: $(FFTW_DIR)/$(FFTW_ZIP) # Шаг 2: распаковать
@echo ">>> Распаковываем $(FFTW_ZIP)..." $(FFTW_SRC)/configure: third_party/$(FFTW_TAR)
unzip -o -q -d $(FFTW_DIR) $< @echo ">>> Распаковываем FFTW..."
tar xf $< -C third_party/
@touch $@ @touch $@
# MinGW ld ищет libfoo.dll.a при -lfoo, поэтому такое имя # Шаг 3a: configure для double precision
$(FFTW_DIR)/libfftw3-3.dll.a: $(FFTW_DIR)/fftw3.h $(FFTW_BUILD_D)/Makefile: $(FFTW_SRC)/configure
@echo ">>> Генерируем import lib: libfftw3-3.dll.a" mkdir -p $(FFTW_BUILD_D)
$(DLLTOOL) --input-def $(FFTW_DIR)/libfftw3-3.def \ @echo ">>> Конфигурируем FFTW (double)..."
--dllname libfftw3-3.dll \ cd $(FFTW_BUILD_D) && $(FFTW_SRC)/configure $(FFTW_CONF_COMMON)
--output-lib $@
$(FFTW_DIR)/libfftw3f-3.dll.a: $(FFTW_DIR)/fftw3.h # Шаг 4a: сборка и установка double precision
@echo ">>> Генерируем import lib: libfftw3f-3.dll.a" fftw-double: $(FFTW_LIB_D)
$(DLLTOOL) --input-def $(FFTW_DIR)/libfftw3f-3.def \ $(FFTW_LIB_D): $(FFTW_BUILD_D)/Makefile
--dllname libfftw3f-3.dll \ @echo ">>> Собираем FFTW (double)..."
--output-lib $@ $(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 библиотеки (кросс-компиляция) ─────────────────────────────────
third-party: $(NR34_DEPS) third-party: $(NR34_DEPS)
third_party/rnnoise/librnnoise.a: $(RNNOISE_LIB): $(RNNOISE_OBJS) | $(OUTDIR)
@echo ">>> Собираем rnnoise для Windows..." $(AR) rv $@ $^
$(MAKE) -C third_party/rnnoise \ $(RANLIB) $@
CC="$(CC)" AR="$(AR)" RANLIB="$(RANLIB)" \
CFLAGS="$(CFLAGS) -Iinclude -Isrc"
# libspecbleach/Makefile жёстко вызывает pkg-config через FFTWINCLUDE=`...` $(SPECBLEACH_LIB): $(SPECBLEACH_OBJS) | $(OUTDIR)
# и использует результат только в CFLAGS?= — переопределяем CFLAGS целиком, $(AR) rv $@ $^
# чтобы обойти pkg-config и прописать путь к нашему FFTW вручную. $(RANLIB) $@
third_party/libspecbleach/libspecbleach.a: $(FFTW_DIR)/fftw3.h
@echo ">>> Собираем libspecbleach для Windows..." # % в GNU make совпадает через '/', поэтому паттерн покрывает вложенные пути
$(MAKE) -C third_party/libspecbleach \ $(OBJDIR)/rnnoise/%.o: $(RNNOISE_DIR)/%.c
CC="$(CC)" AR="$(AR)" RANLIB="$(RANLIB)" \ @mkdir -p $(@D)
CFLAGS="$(CFLAGS) -I$(FFTW_ABS) -Iinclude -Isrc -Isrc/shared" $(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): $(OUTDIR) $(OBJDIR):
mkdir -p $@ mkdir -p $@
$(DLL): $(OBJS) $(NR34_DEPS) | $(OUTDIR) $(DLL): $(OBJS) $(NR34_DEPS) $(FFTW_LIB_D) $(FFTW_LIB_F) | $(OUTDIR)
@echo ">>> Линкуем $@..." @echo ">>> Линкуем $@..."
$(CC) -shared \ $(CC) -shared \
-Wl,--out-implib,$(IMPLIB) \ -Wl,--out-implib,$(IMPLIB) \
$(LDFLAGS) \ $(LDFLAGS) \
-o $@ \ -o $@ \
$(OBJS) $(LINK_LIBS) $(OBJS) $(LINK_LIBS)
@cp $(FFTW_DIR)/libfftw3-3.dll $(FFTW_DIR)/libfftw3f-3.dll $(OUTDIR)/ @cp $(FFTW_DLL_D) $(FFTW_DLL_F) $(OUTDIR)/
@echo "" @echo ""
@echo "=== Готово ===" @echo "=== Готово ==="
@echo " DLL: $(DLL)" @echo " DLL: $(DLL)"
@@ -180,18 +279,16 @@ $(STATIC_LIB): $(OBJS) | $(OUTDIR)
@echo " Статическая: $(STATIC_LIB)" @echo " Статическая: $(STATIC_LIB)"
@echo " Линкуйте приложение с:" @echo " Линкуйте приложение с:"
@echo " $(STATIC_LIB) $(NR34_DEPS)" @echo " $(STATIC_LIB) $(NR34_DEPS)"
@echo " -L$(FFTW_ABS) -lfftw3-3 -lfftw3f-3 -lavrt -lm" @echo " -L$(FFTW_INST)/lib -lfftw3 -lfftw3f -lavrt -lm"
$(OBJDIR)/%.o: %.c | $(OBJDIR) $(OBJDIR)/%.o: %.c | $(OBJDIR) $(FFTW_LIB_D)
$(COMPILE) -c -o $@ $< $(COMPILE) -c -o $@ $<
# ── Очистка ─────────────────────────────────────────────────────────────────── # ── Очистка ───────────────────────────────────────────────────────────────────
clean: clean:
-rm -rf $(OBJDIR) $(OUTDIR) -rm -rf $(OBJDIR) $(OUTDIR)
-$(MAKE) -C third_party/rnnoise clean
-$(MAKE) -C third_party/libspecbleach clean
# distclean удаляет и скачанный FFTW — при следующем запуске скачается заново # distclean удаляет собранный и скачанный FFTW — при следующем запуске пересоберётся
distclean: clean distclean: clean
-rm -rf $(FFTW_DIR) -rm -rf $(FFTW_INST) $(FFTW_SRC) third_party/$(FFTW_TAR)
README.md
+12 -16
View File
@@ -17,7 +17,7 @@ wdsp/
├── java/ ├── java/
│ └── org/openhpsdr/dsp/Wdsp.java │ └── org/openhpsdr/dsp/Wdsp.java
├── third_party/ ├── third_party/
│ ├── fftw/ — FFTW source (download manually, Android only) │ ├── fftw-3.3.11/ — FFTW source (downloaded automatically)
│ ├── fftw-win64/ — FFTW Windows binaries (downloaded automatically) │ ├── fftw-win64/ — FFTW Windows binaries (downloaded automatically)
│ ├── rnnoise/ — RNNoise noise suppression │ ├── rnnoise/ — RNNoise noise suppression
│ └── libspecbleach/ — spectral noise reduction │ └── libspecbleach/ — spectral noise reduction
@@ -87,23 +87,25 @@ make clean
## Windows build (кросс-компиляция с Linux) ## Windows build (кросс-компиляция с Linux)
Сборка выполняется на Linux с помощью MinGW-w64. FFTW скачивается автоматически. Сборка выполняется на Linux с помощью MinGW-w64. FFTW скачивается и собирается из исходников автоматически.
### Зависимости ### Зависимости
```bash ```bash
# Arch # Arch
sudo pacman -S mingw-w64-gcc unzip sudo pacman -S mingw-w64-gcc # тянет mingw-w64-headers автоматически
# Debian / Ubuntu # Debian / Ubuntu — нужен мета-пакет mingw-w64, а не только компилятор
sudo apt install gcc-mingw-w64-x86-64 binutils-mingw-w64-x86-64 unzip sudo apt install mingw-w64
# Fedora # Fedora
sudo dnf install mingw64-gcc mingw64-binutils unzip sudo dnf install mingw64-gcc # тянет mingw64-headers автоматически
``` ```
Также нужен `wget` или `curl`. Также нужен `wget` или `curl`.
> **Debian/Ubuntu:** пакет `gcc-mingw-w64-x86-64` содержит только компилятор без Windows API заголовков — `Windows.h` не будет найден. Используйте `mingw-w64`.
### Сборка ### Сборка
```bash ```bash
@@ -169,15 +171,9 @@ make -f Makefile.windows distclean # + удалить скачанный FFTW
### Prerequisites ### Prerequisites
1. **Android NDK** r23 or newer 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 FFTW 3.3.11 скачивается и собирается автоматически. Чтобы использовать уже скачанный исходник, передайте `FFTW_SRC=/path/to/fftw-3.3.11`.
wget https://www.fftw.org/fftw-3.3.11.tar.gz
tar xf fftw-3.3.11.tar.gz
mv fftw-3.3.11 third_party/fftw
```
3. **Java compiler** — for building the `.class` file (Android Studio's JBR or any JDK)
### Build ### Build
@@ -227,8 +223,8 @@ make -f Makefile.android \
### Clean ### Clean
```bash ```bash
make -f Makefile.android clean make -f Makefile.android clean # удаляет obj/android/ и lib/android/
# removes obj/android/ and lib/android/ make -f Makefile.android distclean # + удаляет скачанный FFTW
``` ```
--- ---
RXA.c
+32
View File
@@ -233,6 +233,23 @@ void create_rxa (int channel)
max(2048, ch[channel].dsp_size), // number of coefficients for noise filter max(2048, ch[channel].dsp_size), // number of coefficients for noise filter
0); // minimum phase flag 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 // snba
rxa[channel].snba.p = create_snba ( rxa[channel].snba.p = create_snba (
0, // run 0, // run
@@ -579,6 +596,7 @@ void destroy_rxa (int channel)
destroy_anf (rxa[channel].anf.p); destroy_anf (rxa[channel].anf.p);
destroy_eqp (rxa[channel].eqp.p); destroy_eqp (rxa[channel].eqp.p);
destroy_snba (rxa[channel].snba.p); destroy_snba (rxa[channel].snba.p);
destroy_wfmd (rxa[channel].wfmd.p);
destroy_fmsq (rxa[channel].fmsq.p); destroy_fmsq (rxa[channel].fmsq.p);
destroy_fmd (rxa[channel].fmd.p); destroy_fmd (rxa[channel].fmd.p);
destroy_amd (rxa[channel].amd.p); destroy_amd (rxa[channel].amd.p);
@@ -614,6 +632,7 @@ void flush_rxa (int channel)
flush_amd (rxa[channel].amd.p); flush_amd (rxa[channel].amd.p);
flush_fmd (rxa[channel].fmd.p); flush_fmd (rxa[channel].fmd.p);
flush_fmsq (rxa[channel].fmsq.p); flush_fmsq (rxa[channel].fmsq.p);
flush_wfmd (rxa[channel].wfmd.p);
flush_snba (rxa[channel].snba.p); flush_snba (rxa[channel].snba.p);
flush_eqp (rxa[channel].eqp.p); flush_eqp (rxa[channel].eqp.p);
flush_anf (rxa[channel].anf.p); flush_anf (rxa[channel].anf.p);
@@ -649,6 +668,7 @@ void xrxa (int channel)
xamd (rxa[channel].amd.p); xamd (rxa[channel].amd.p);
xfmd (rxa[channel].fmd.p); xfmd (rxa[channel].fmd.p);
xfmsq (rxa[channel].fmsq.p); xfmsq (rxa[channel].fmsq.p);
xwfmd (rxa[channel].wfmd.p);
xbpsnbain (rxa[channel].bpsnba.p, 1); xbpsnbain (rxa[channel].bpsnba.p, 1);
xbpsnbaout (rxa[channel].bpsnba.p, 1); xbpsnbaout (rxa[channel].bpsnba.p, 1);
xsnba (rxa[channel].snba.p); xsnba (rxa[channel].snba.p);
@@ -733,6 +753,7 @@ void setDSPSamplerate_rxa (int channel)
setSamplerate_fmd (rxa[channel].fmd.p, ch[channel].dsp_rate); 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); 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_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_snba (rxa[channel].snba.p, ch[channel].dsp_rate);
setSamplerate_eqp (rxa[channel].eqp.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); 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); 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); 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); 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); setBuffers_snba (rxa[channel].snba.p, rxa[channel].midbuff, rxa[channel].midbuff);
setSize_snba (rxa[channel].snba.p, ch[channel].dsp_size); setSize_snba (rxa[channel].snba.p, ch[channel].dsp_size);
setBuffers_eqp (rxa[channel].eqp.p, rxa[channel].midbuff, rxa[channel].midbuff); 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].mode = mode;
rxa[channel].amd.p->run = 0; rxa[channel].amd.p->run = 0;
rxa[channel].fmd.p->run = 0; rxa[channel].fmd.p->run = 0;
rxa[channel].wfmd.p->run = 0;
rxa[channel].agc.p->run = 1; rxa[channel].agc.p->run = 1;
switch (mode) switch (mode)
{ {
@@ -875,6 +899,10 @@ void SetRXAMode (int channel, int mode)
rxa[channel].fmd.p->run = 1; rxa[channel].fmd.p->run = 1;
rxa[channel].agc.p->run = 0; rxa[channel].agc.p->run = 0;
break; break;
case RXA_WFM:
rxa[channel].wfmd.p->run = 1;
rxa[channel].agc.p->run = 0;
break;
default: default:
break; break;
@@ -961,6 +989,7 @@ void RXAbpsnbaCheck (int channel, int mode, int notch_run)
run_notches = 0; run_notches = 0;
break; break;
case RXA_FM: case RXA_FM:
case RXA_WFM:
f_low = +a->abs_low_freq; f_low = +a->abs_low_freq;
f_high = +a->abs_high_freq; f_high = +a->abs_high_freq;
run_notches = 0; run_notches = 0;
@@ -1010,6 +1039,7 @@ void RXAbpsnbaSet (int channel)
a->position = 1; a->position = 1;
break; break;
case RXA_FM: case RXA_FM:
case RXA_WFM:
a->run = rxa[channel].snba.p->run; a->run = rxa[channel].snba.p->run;
a->position = 1; a->position = 1;
break; break;
@@ -1046,6 +1076,7 @@ void RXASetNC (int channel, int nc)
SetRXAFMSQNC (channel, nc); SetRXAFMSQNC (channel, nc);
SetRXAFMNCde (channel, nc); SetRXAFMNCde (channel, nc);
SetRXAFMNCaud (channel, nc); SetRXAFMNCaud (channel, nc);
SetRXAWFMNCaud (channel, nc);
SetChannelState (channel, oldstate, 0); SetChannelState (channel, oldstate, 0);
} }
@@ -1059,4 +1090,5 @@ void RXASetMP (int channel, int mp)
SetRXAFMSQMP (channel, mp); SetRXAFMSQMP (channel, mp);
SetRXAFMMPde (channel, mp); SetRXAFMMPde (channel, mp);
SetRXAFMMPaud (channel, mp); SetRXAFMMPaud (channel, mp);
SetRXAWFMMPaud (channel, mp);
} }
RXA.h
+6 -1
View File
@@ -41,7 +41,8 @@ enum rxaMode
RXA_SPEC, RXA_SPEC,
RXA_DIGL, RXA_DIGL,
RXA_SAM, RXA_SAM,
RXA_DRM RXA_DRM,
RXA_WFM
}; };
enum rxaMeterType enum rxaMeterType
@@ -121,6 +122,10 @@ struct _rxa
FMSQ p; FMSQ p;
} fmsq; } fmsq;
struct struct
{
WFMD p;
} wfmd;
struct
{ {
EQP p; EQP p;
} eqp; } eqp;
TXA.c
+34 -4
View File
@@ -358,6 +358,21 @@ void create_txa (int channel)
max(2048, ch[channel].dsp_size), // number coefficients for bandpass filter max(2048, ch[channel].dsp_size), // number coefficients for bandpass filter
0); // minimum phase flag 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 ( txa[channel].gen1.p = create_gen (
0, // run 0, // run
ch[channel].dsp_size, // buffer size ch[channel].dsp_size, // buffer size
@@ -490,6 +505,7 @@ void destroy_txa (int channel)
destroy_meter (txa[channel].alcmeter.p); destroy_meter (txa[channel].alcmeter.p);
destroy_uslew (txa[channel].uslew.p); destroy_uslew (txa[channel].uslew.p);
destroy_gen (txa[channel].gen1.p); destroy_gen (txa[channel].gen1.p);
destroy_wfmmod (txa[channel].wfmmod.p);
destroy_fmmod (txa[channel].fmmod.p); destroy_fmmod (txa[channel].fmmod.p);
destroy_ammod (txa[channel].ammod.p); destroy_ammod (txa[channel].ammod.p);
destroy_wcpagc (txa[channel].alc.p); destroy_wcpagc (txa[channel].alc.p);
@@ -544,6 +560,7 @@ void flush_txa (int channel)
flush_wcpagc (txa[channel].alc.p); flush_wcpagc (txa[channel].alc.p);
flush_ammod (txa[channel].ammod.p); flush_ammod (txa[channel].ammod.p);
flush_fmmod (txa[channel].fmmod.p); flush_fmmod (txa[channel].fmmod.p);
flush_wfmmod (txa[channel].wfmmod.p);
flush_gen (txa[channel].gen1.p); flush_gen (txa[channel].gen1.p);
flush_uslew (txa[channel].uslew.p); flush_uslew (txa[channel].uslew.p);
flush_meter (txa[channel].alcmeter.p); flush_meter (txa[channel].alcmeter.p);
@@ -580,6 +597,7 @@ void xtxa (int channel)
xammod (txa[channel].ammod.p); // AM Modulator xammod (txa[channel].ammod.p); // AM Modulator
xemphp (txa[channel].preemph.p, 1); // FM pre-emphasis (second option) xemphp (txa[channel].preemph.p, 1); // FM pre-emphasis (second option)
xfmmod (txa[channel].fmmod.p); // FM Modulator 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) xgen (txa[channel].gen1.p); // output signal generator (TUN and Two-tone)
xuslew (txa[channel].uslew.p); // up-slew for AM, FM, and gens xuslew (txa[channel].uslew.p); // up-slew for AM, FM, and gens
xmeter (txa[channel].alcmeter.p); // ALC Meter 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_wcpagc (txa[channel].alc.p, ch[channel].dsp_rate);
setSamplerate_ammod (txa[channel].ammod.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_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_gen (txa[channel].gen1.p, ch[channel].dsp_rate);
setSamplerate_uslew (txa[channel].uslew.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); 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); setSize_ammod (txa[channel].ammod.p, ch[channel].dsp_size);
setBuffers_fmmod (txa[channel].fmmod.p, txa[channel].midbuff, txa[channel].midbuff); setBuffers_fmmod (txa[channel].fmmod.p, txa[channel].midbuff, txa[channel].midbuff);
setSize_fmmod (txa[channel].fmmod.p, ch[channel].dsp_size); 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); setBuffers_gen (txa[channel].gen1.p, txa[channel].midbuff, txa[channel].midbuff);
setSize_gen (txa[channel].gen1.p, ch[channel].dsp_size); setSize_gen (txa[channel].gen1.p, ch[channel].dsp_size);
setBuffers_uslew (txa[channel].uslew.p, txa[channel].midbuff, txa[channel].midbuff); 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].mode = mode;
txa[channel].ammod.p->run = 0; txa[channel].ammod.p->run = 0;
txa[channel].fmmod.p->run = 0; txa[channel].fmmod.p->run = 0;
txa[channel].wfmmod.p->run = 0;
txa[channel].preemph.p->run = 0; txa[channel].preemph.p->run = 0;
switch (mode) switch (mode)
{ {
@@ -779,6 +801,10 @@ void SetTXAMode (int channel, int mode)
txa[channel].fmmod.p->run = 1; txa[channel].fmmod.p->run = 1;
txa[channel].preemph.p->run = 1; txa[channel].preemph.p->run = 1;
break; break;
case TXA_WFM:
// wfmmod carries its own RC pre-emphasis; the shared emphp stays off
txa[channel].wfmmod.p->run = 1;
break;
default: default:
break; break;
@@ -818,10 +844,11 @@ void TXAResCheck (int channel)
int TXAUslewCheck (int channel) int TXAUslewCheck (int channel)
{ {
return (txa[channel].ammod.p->run == 1) || return (txa[channel].ammod.p->run == 1) ||
(txa[channel].fmmod.p->run == 1) || (txa[channel].fmmod.p->run == 1) ||
(txa[channel].gen0.p->run == 1) || (txa[channel].wfmmod.p->run == 1) ||
(txa[channel].gen1.p->run == 1); (txa[channel].gen0.p->run == 1) ||
(txa[channel].gen1.p->run == 1);
} }
void TXASetupBPFilters (int channel) void TXASetupBPFilters (int channel)
@@ -855,6 +882,7 @@ void TXASetupBPFilters (int channel)
case TXA_AM: case TXA_AM:
case TXA_SAM: case TXA_SAM:
case TXA_FM: case TXA_FM:
case TXA_WFM:
if (txa[channel].compressor.p->run) if (txa[channel].compressor.p->run)
{ {
CalcBandpassFilter (txa[channel].bp0.p, 0.0, txa[channel].f_high, 2.0); 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); SetTXAFMEmphNC (channel, nc);
SetTXAEQNC (channel, nc); SetTXAEQNC (channel, nc);
SetTXAFMNC (channel, nc); SetTXAFMNC (channel, nc);
SetTXAWFMNC (channel, nc);
SetTXACFIRNC (channel, nc); SetTXACFIRNC (channel, nc);
SetChannelState (channel, oldstate, 0); SetChannelState (channel, oldstate, 0);
} }
@@ -925,6 +954,7 @@ void TXASetMP (int channel, int mp)
SetTXAFMEmphMP (channel, mp); SetTXAFMEmphMP (channel, mp);
SetTXAEQMP (channel, mp); SetTXAEQMP (channel, mp);
SetTXAFMMP (channel, mp); SetTXAFMMP (channel, mp);
SetTXAWFMMP (channel, mp);
} }
PORT PORT
TXA.h
+6 -1
View File
@@ -43,7 +43,8 @@ enum txaMode
TXA_SAM, TXA_SAM,
TXA_DRM, TXA_DRM,
TXA_AM_LSB, TXA_AM_LSB,
TXA_AM_USB TXA_AM_USB,
TXA_WFM
}; };
enum txaMeterType enum txaMeterType
@@ -135,6 +136,10 @@ struct _txa
FMMOD p; FMMOD p;
} fmmod; } fmmod;
struct struct
{
WFMMOD p;
} wfmmod;
struct
{ {
SIPHON p; SIPHON p;
} sip1; } sip1;
analyzer.c
+19 -1
View File
@@ -1087,7 +1087,15 @@ void __cdecl sendbuf(void *arg)
LeaveCriticalSection(&(a->BufferControlSection[a->ss][a->LO])); 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); InterlockedBitTestAndReset(&a->dispatcher, 0);
_endthread(); _endthread();
@@ -1202,6 +1210,9 @@ void SetAnalyzer ( int disp, // display identifier
EnterCriticalSection(&a->SetAnalyzerSection); EnterCriticalSection(&a->SetAnalyzerSection);
a->end_dispatcher = 1; 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)) while (InterlockedAnd(&a->dispatcher, 1))
Sleep(1); Sleep(1);
a->stop = 1; a->stop = 1;
@@ -1347,6 +1358,7 @@ void XCreateAnalyzer( int disp,
a->hSnapEvent[i][j] = CreateEvent(NULL, FALSE, FALSE, TEXT("snap")); a->hSnapEvent[i][j] = CreateEvent(NULL, FALSE, FALSE, TEXT("snap"));
a->snap[i][j] = 0; a->snap[i][j] = 0;
} }
a->Sem_BuffReady = CreateSemaphore(0, 0, 1000, 0);
InitializeCriticalSectionAndSpinCount(&a->ResampleSection, 0); InitializeCriticalSectionAndSpinCount(&a->ResampleSection, 0);
InitializeCriticalSectionAndSpinCount(&a->SetAnalyzerSection, 0); InitializeCriticalSectionAndSpinCount(&a->SetAnalyzerSection, 0);
InitializeCriticalSectionAndSpinCount(&a->StitchSection, 0); InitializeCriticalSectionAndSpinCount(&a->StitchSection, 0);
@@ -1434,6 +1446,7 @@ void DestroyAnalyzer(int disp)
int i, j; int i, j;
a->end_dispatcher = 1; a->end_dispatcher = 1;
ReleaseSemaphore(a->Sem_BuffReady, 1, 0);
while (InterlockedAnd(&a->dispatcher, 1)) while (InterlockedAnd(&a->dispatcher, 1))
Sleep(1); Sleep(1);
@@ -1497,6 +1510,7 @@ void DestroyAnalyzer(int disp)
for (i = 0; i < a->max_stitch; i++) for (i = 0; i < a->max_stitch; i++)
for (j = 0; j < a->max_num_fft; j++) for (j = 0; j < a->max_num_fft; j++)
CloseHandle(a->hSnapEvent[i][j]); CloseHandle(a->hSnapEvent[i][j]);
CloseHandle(a->Sem_BuffReady);
_aligned_free ((void *) a->pnum_threads); _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! 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; a->IQin_index[ss][LO] = 0;
ReleaseSemaphore(a->Sem_BuffReady, 1, 0); // new samples: let the dispatcher run
if (!InterlockedAnd(&a->dispatcher, 1)) if (!InterlockedAnd(&a->dispatcher, 1))
{ {
InterlockedBitTestAndSet (&a->dispatcher, 0); 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! 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; a->IQin_index[ss][LO] = 0;
ReleaseSemaphore(a->Sem_BuffReady, 1, 0); // new samples: let the dispatcher run
if (!InterlockedAnd(&a->dispatcher, 1)) if (!InterlockedAnd(&a->dispatcher, 1))
{ {
InterlockedBitTestAndSet(&a->dispatcher, 0); 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! 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; a->IQin_index[ss][LO] = 0;
ReleaseSemaphore(a->Sem_BuffReady, 1, 0); // new samples: let the dispatcher run
if (!InterlockedAnd(&a->dispatcher, 1)) if (!InterlockedAnd(&a->dispatcher, 1))
{ {
InterlockedBitTestAndSet(&a->dispatcher, 0); 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! 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; a->IQin_index[ss][LO] = 0;
ReleaseSemaphore(a->Sem_BuffReady, 1, 0); // new samples: let the dispatcher run
if (!InterlockedAnd(&a->dispatcher, 1)) if (!InterlockedAnd(&a->dispatcher, 1))
{ {
InterlockedBitTestAndSet(&a->dispatcher, 0); InterlockedBitTestAndSet(&a->dispatcher, 0);
analyzer.h
+3
View File
@@ -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 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 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 PB_ControlsSection[dMAX_PIXOUTS];
CRITICAL_SECTION SetAnalyzerSection; CRITICAL_SECTION SetAnalyzerSection;
CRITICAL_SECTION BufferControlSection[dMAX_STITCH][dMAX_NUM_FFT]; CRITICAL_SECTION BufferControlSection[dMAX_STITCH][dMAX_NUM_FFT];
anf.c
+94 -27
View File
@@ -26,6 +26,45 @@ warren@wpratt.com
#include "comm.h" #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 ( ANF create_anf (
int run, int run,
int position, int position,
@@ -81,53 +120,81 @@ void destroy_anf (ANF a)
void xanf(ANF a, int position) void xanf(ANF a, int position)
{ {
int i, j, idx; int i;
double c0, c1; double c0, c1;
double y, error, sigma, inv_sigp; double y, error, sigma, inv_sigp;
double nel, nev; double nel, nev;
if (a->run && (a->position == position)) 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; dsamp = in_buff[2 * i + 0];
sigma = 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); inv_sigp = 1.0 / (sigma + 1e-10);
error = a->d[a->in_idx] - y; error = dsamp - y;
a->out_buff[2 * i + 0] = error; out_buff[2 * i + 0] = error;
a->out_buff[2 * i + 1] = 0.0; out_buff[2 * i + 1] = 0.0;
if((nel = error * (1.0 - a->two_mu * sigma * inv_sigp)) < 0.0) nel = -nel; if((nel = error * (1.0 - 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((nev = dsamp - (1.0 - two_mu * ngamma) * y - two_mu * error * sigma * inv_sigp) < 0.0) nev = -nev;
if (nev < nel) 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 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; c0 = 1.0 - two_mu * ngamma;
c1 = a->two_mu * error * inv_sigp; c1 = two_mu * error * inv_sigp;
for (j = 0; j < a->n_taps; j++) anf_update (w, d + base, n1, c0, c1);
{ if (n1 < n_taps)
idx = (a->in_idx + j + a->delay) & a->mask; anf_update (w + n1, d, n_taps - n1, c0, c1);
a->w[j] = c0 * a->w[j] + c1 * a->d[idx];
} in_idx = (in_idx + mask) & mask;
a->in_idx = (a->in_idx + a->mask) & a->mask;
} }
a->in_idx = in_idx;
a->lidx = lidx;
a->ngamma = ngamma;
} }
else if (a->in_buff != a->out_buff) else if (a->in_buff != a->out_buff)
memcpy (a->out_buff, a->in_buff, a->buff_size * sizeof (complex)); memcpy (a->out_buff, a->in_buff, a->buff_size * sizeof (complex));
anr.c
+94 -27
View File
@@ -26,6 +26,45 @@ warren@wpratt.com
#include "comm.h" #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 ( ANR create_anr (
int run, int run,
int position, int position,
@@ -81,53 +120,81 @@ void destroy_anr (ANR a)
void xanr (ANR a, int position) void xanr (ANR a, int position)
{ {
int i, j, idx; int i;
double c0, c1; double c0, c1;
double y, error, sigma, inv_sigp; double y, error, sigma, inv_sigp;
double nel, nev; double nel, nev;
if (a->run && (a->position == position)) 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; dsamp = in_buff[2 * i + 0];
sigma = 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); inv_sigp = 1.0 / (sigma + 1e-10);
error = a->d[a->in_idx] - y; error = dsamp - y;
a->out_buff[2 * i + 0] = y; out_buff[2 * i + 0] = y;
a->out_buff[2 * i + 1] = 0.0; out_buff[2 * i + 1] = 0.0;
if((nel = error * (1.0 - a->two_mu * sigma * inv_sigp)) < 0.0) nel = -nel; if((nel = error * (1.0 - 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((nev = dsamp - (1.0 - two_mu * ngamma) * y - two_mu * error * sigma * inv_sigp) < 0.0) nev = -nev;
if (nev < nel) 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 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; c0 = 1.0 - two_mu * ngamma;
c1 = a->two_mu * error * inv_sigp; c1 = two_mu * error * inv_sigp;
for (j = 0; j < a->n_taps; j++) anr_update (w, d + base, n1, c0, c1);
{ if (n1 < n_taps)
idx = (a->in_idx + j + a->delay) & a->mask; anr_update (w + n1, d, n_taps - n1, c0, c1);
a->w[j] = c0 * a->w[j] + c1 * a->d[idx];
} in_idx = (in_idx + mask) & mask;
a->in_idx = (a->in_idx + a->mask) & a->mask;
} }
a->in_idx = in_idx;
a->lidx = lidx;
a->ngamma = ngamma;
} }
else if (a->in_buff != a->out_buff) else if (a->in_buff != a->out_buff)
memcpy (a->out_buff, a->in_buff, a->buff_size * sizeof (complex)); memcpy (a->out_buff, a->in_buff, a->buff_size * sizeof (complex));
calcc.c
+171 -2
View File
@@ -138,6 +138,7 @@ CALCC create_calcc (int channel, int runcal, int size, int rate, int ints, int s
a->stbl = stbl; a->stbl = stbl;
a->npsamps = npsamps; a->npsamps = npsamps;
a->alpha = alpha; a->alpha = alpha;
a->outlier_sigma = 0.0;
a->info = (int *) malloc0 (16 * sizeof (int)); a->info = (int *) malloc0 (16 * sizeof (int));
a->binfo = (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; 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) void calc (CALCC a)
{ {
int i; int i;
@@ -336,21 +456,60 @@ void calc (CALCC a)
double tvec[3]; double tvec[3];
double txrxcoefs[4 * 2]; double txrxcoefs[4 * 2];
double rx_scale; double rx_scale;
int xb_ok = 0;
double* tx_filt;
double* rx_filt;
int n_filt = 0;
if (a->ints < 16) rints = 1; if (a->ints < 16) rints = 1;
else rints = 2; else rints = 2;
ix = rints - 1; ix = rints - 1;
for (i = 0; i <= rints; i++) for (i = 0; i <= rints; i++)
tvec[i] = (double)i / (double)rints / a->hw_scale; tvec[i] = (double)i / (double)rints / a->hw_scale;
dx = tvec[rints] - tvec[rints - 1]; 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]); rxscheck (rints, tvec, txrxcoefs, &a->binfo[7]);
if ((a->binfo[0] == 0) && (a->binfo[7] == 0)) 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]))); 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; a->scOK = 0;
goto cleanup; goto cleanup;
} }
if (a->stbl && _InterlockedAnd (&a->ctrl.running, 1)) if (a->stbl && _InterlockedAnd (&a->ctrl.running, 1))
a->rx_scale = a->alpha * a->rx_scale + (1.0 - a->alpha) * rx_scale; a->rx_scale = a->alpha * a->rx_scale + (1.0 - a->alpha) * rx_scale;
else else
@@ -1046,6 +1205,16 @@ void SetPSPtol (int channel, double ptol)
LeaveCriticalSection (&txa[channel].calcc.cs_update); 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 PORT
void GetPSDisp (int channel, double* x, double* ym, double* yc, double* ys, double* cm, double* cc, double* cs) void GetPSDisp (int channel, double* x, double* ym, double* yc, double* ys, double* cm, double* cc, double* cs)
{ {
calcc.h
+1
View File
@@ -48,6 +48,7 @@ typedef struct _calcc
double hw_scale; double hw_scale;
double rx_scale; double rx_scale;
double alpha; double alpha;
double outlier_sigma;
int tsamps; int tsamps;
double* env_TX; double* env_TX;
cfcomp.c
+53 -24
View File
@@ -336,18 +336,36 @@ void xcfcomp (CFCOMP a, int pos)
if (a->run && pos == a->position) if (a->run && pos == a->position)
{ {
int i, j, k, sbuff, sbegin; 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->inaccum[j] = a->in[i];
a->iainidx = (a->iainidx + 1) % a->iasize; if (++j == iasize) j = 0;
} }
a->nsamps += a->bsize; a->iainidx = j;
while (a->nsamps >= a->fsize) a->nsamps += bsize;
while (a->nsamps >= fsize)
{ {
for (i = 0, j = a->iaoutidx; i < a->fsize; i++, j = (j + 1) % a->iasize) int n1 = iasize - a->iaoutidx;
a->forfftin[i] = a->pregain * a->window[i] * a->inaccum[j]; if (n1 > fsize) n1 = fsize;
a->iaoutidx = (a->iaoutidx + a->incr) % a->iasize; for (i = 0; i < n1; i++)
a->nsamps -= a->incr; 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); fftw_execute (a->Rfor);
calc_mask(a); calc_mask(a);
for (i = 0; i < a->msize; i++) 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]; a->revfftin[2 * i + 1] = a->mask[i] * a->forfftout[2 * i + 1];
} }
fftw_execute (a->Rrev); fftw_execute (a->Rrev);
for (i = 0; i < a->fsize; i++) for (i = 0; i < fsize; i++)
a->save[a->saveidx][i] = a->postgain * a->window[i] * a->revfftout[i]; a->save[a->saveidx][i] = postgain * 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; const double* WDSP_RESTRICT sv;
sbegin = a->incr * (a->ovrlp - i); double* WDSP_RESTRICT oa = a->outaccum;
for (j = sbegin, k = a->oainidx; j < a->incr + sbegin; j++, k = (k + 1) % a->oasize) 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) for (j = 0; j < m1; j++) oa[k + j] = sv[j];
a->outaccum[k] = a->save[sbuff][j]; for (; j < incr; j++) oa[j - m1] = sv[j];
else }
a->outaccum[k] += a->save[sbuff][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; if (++a->saveidx == ovrlp) a->saveidx = 0;
a->oainidx = (a->oainidx + a->incr) % a->oasize; 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->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) else if (a->out != a->in)
memcpy (a->out, a->in, a->bsize * sizeof (complex)); memcpy (a->out, a->in, a->bsize * sizeof (complex));
comm.h
+17 -1
View File
@@ -33,7 +33,7 @@ warren@wpratt.com
#endif #endif
#ifdef _WIN32 #ifdef _WIN32
#include <Windows.h> #include <windows.h>
#include <process.h> #include <process.h>
#include <intrin.h> #include <intrin.h>
#endif #endif
@@ -43,6 +43,11 @@ warren@wpratt.com
#ifdef _WIN32 #ifdef _WIN32
#include <avrt.h> #include <avrt.h>
#endif #endif
#ifndef WDSP_FPE_GUARD
#define WDSP_FPE_GUARD ((void)0)
#define WDSP_FPE_RESTORE ((void)0)
#endif
#include "fftw3.h" #include "fftw3.h"
#include "amd.h" #include "amd.h"
@@ -70,12 +75,15 @@ warren@wpratt.com
#include "sbnr.h" // NR3 + NR4 support #include "sbnr.h" // NR3 + NR4 support
#include "emph.h" #include "emph.h"
#include "eq.h" #include "eq.h"
#include "fastmath.h"
#include "fcurve.h" #include "fcurve.h"
#include "fir.h" #include "fir.h"
#include "firmin.h" #include "firmin.h"
#include "fmd.h" #include "fmd.h"
#include "fmmod.h" #include "fmmod.h"
#include "fmsq.h" #include "fmsq.h"
#include "wfmd.h"
#include "wfmmod.h"
#include "gain.h" #include "gain.h"
#include "gaussian.h" #include "gaussian.h"
#include "gen.h" #include "gen.h"
@@ -145,6 +153,14 @@ warren@wpratt.com
#define PI 3.1415926535897932 #define PI 3.1415926535897932
#define TWOPI 6.2831853071795864 #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 // miscellaneous
typedef double complex[2]; typedef double complex[2];
#define PORT __declspec( dllexport ) #define PORT __declspec( dllexport )
emnr.c
+111 -73
View File
@@ -557,6 +557,7 @@ void calc_emnr(EMNR a)
a->ae.psi = 20.0; a->ae.psi = 20.0;
a->ae.t2 = 0.20; a->ae.t2 = 0.20;
a->ae.nmask = (double *)malloc0(a->ae.msize * sizeof(double)); a->ae.nmask = (double *)malloc0(a->ae.msize * sizeof(double));
a->ae.csum = (double *)malloc0((a->ae.msize + 1) * sizeof(double));
// //
// post2 // post2
a->post2.run = 0; a->post2.run = 0;
@@ -580,6 +581,7 @@ void decalc_emnr(EMNR a)
_aligned_free(a->post2.noise_frame); _aligned_free(a->post2.noise_frame);
_aligned_free(a->post2.w); _aligned_free(a->post2.w);
// ae // ae
_aligned_free(a->ae.csum);
_aligned_free(a->ae.nmask); _aligned_free(a->ae.nmask);
// npl // npl
_aligned_free(a->npl.D); _aligned_free(a->npl.D);
@@ -868,26 +870,26 @@ void aepf(EMNR a)
else else
N = 1 + 2 * (int)(0.5 + a->ae.psi * (1.0 - zetaT / a->ae.zetaThresh)); N = 1 + 2 * (int)(0.5 + a->ae.psi * (1.0 - zetaT / a->ae.zetaThresh));
n = N / 2; 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; const int msize = a->ae.msize;
for (m = 0; m <= 2 * k; m++) const double* WDSP_RESTRICT mask = a->mask;
a->ae.nmask[k] += a->mask[m]; double* WDSP_RESTRICT nmask = a->ae.nmask;
a->ae.nmask[k] /= (double)(2 * k + 1); double* WDSP_RESTRICT csum = a->ae.csum;
}
for (k = n; k < (a->ae.msize - n); k++) csum[0] = 0.0;
{ for (k = 0; k < msize; k++)
a->ae.nmask[k] = 0.0; csum[k + 1] = csum[k] + mask[k];
for (m = k - n; m <= (k + n); m++)
a->ae.nmask[k] += a->mask[m]; for (k = 0; k < n; k++) // window [0, 2k]
a->ae.nmask[k] /= (double)N; nmask[k] = (csum[2 * k + 1] - csum[0]) / (double)(2 * k + 1);
} for (k = n; k < (msize - n); k++) // window [k-n, k+n]
for (k = a->ae.msize - n; k < a->ae.msize; k++) nmask[k] = (csum[k + n + 1] - csum[k - n]) / (double)N;
{ for (k = msize - n; k < msize; k++) // window [2k+1-msize, msize-1]
a->ae.nmask[k] = 0.0; nmask[k] = (csum[msize] - csum[2 * k + 1 - msize]) / (double)(2 * (msize - k) - 1);
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);
} }
memcpy (a->mask, a->ae.nmask, a->ae.msize * sizeof (double)); memcpy (a->mask, a->ae.nmask, a->ae.msize * sizeof (double));
if (a->g.gain_method == 3 && zetaT < a->ae.t2) if (a->g.gain_method == 3 && zetaT < a->ae.t2)
@@ -998,52 +1000,54 @@ void SetRXAEMNRpost2Rate(int channel, double tc)
* End Post-Processing Functions * * 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 dmin = 0.001;
const double dmax = 1000.0; const double dmax = 1000.0;
if (gamma <= dmin) if (v <= dmin)
{ {
ngamma1 = ngamma2 = 0; *n1 = *n2 = 0;
tg = 0.0; *d = 0.0;
} }
else if (gamma >= dmax) else if (v >= dmax)
{ {
ngamma1 = ngamma2 = 240; *n1 = *n2 = 240;
tg = 60.0; *d = 0.0;
} }
else else
{ {
tg = 10.0 * log10(gamma / dmin); double f = 40.0 * wdsp_log10 (v / dmin); // 4 * (10 * log10)
ngamma1 = (int)(4.0 * tg); int i = (int)f;
ngamma2 = ngamma1 + 1; /* 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; static inline double keyLerp (const double* type, int ngamma1, int ngamma2, double dg,
tx = 0.0; int nxi1, int nxi2, double dx)
} {
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;
return (1.0 - dg) * (1.0 - dx) * type[241 * nxi1 + ngamma1] return (1.0 - dg) * (1.0 - dx) * type[241 * nxi1 + ngamma1]
+ (1.0 - dg) * dx * type[241 * nxi2 + ngamma1] + (1.0 - dg) * dx * type[241 * nxi2 + ngamma1]
+ dg * (1.0 - dx) * type[241 * nxi1 + ngamma2] + dg * (1.0 - dx) * type[241 * nxi1 + ngamma2]
+ dg * dx * type[241 * nxi2 + 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 getZeta( EMNR a, double gamma, double eps, double* zeta)
{ {
int index, i_gamma, i_xi; int index, i_gamma, i_xi;
@@ -1133,13 +1137,20 @@ void calc_gain (EMNR a)
case 2: case 2:
{ {
double gamma, eps_hat, eps_p; 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++) for (k = 0; k < a->g.msize; k++)
{ {
gamma = min(a->g.lambda_y[k] / a->g.lambda_d[k], a->g.gamma_max); 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] 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); + (1.0 - a->g.alpha) * max(gamma - 1.0, a->g.eps_floor);
eps_p = eps_hat / (1.0 - a->g.q); 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_gamma[k] = gamma;
a->g.prev_mask[k] = a->g.mask[k]; 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; int i, j, k, sbuff, sbegin;
double g1; 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->inaccum[j] = a->in[i];
a->iainidx = (a->iainidx + 1) % a->iasize; if (++j == iasize) j = 0;
} }
a->nsamps += a->bsize; a->iainidx = j;
while (a->nsamps >= a->fsize) a->nsamps += bsize;
while (a->nsamps >= fsize)
{ {
for (i = 0, j = a->iaoutidx; i < a->fsize; i++, j = (j + 1) % a->iasize) int n1 = iasize - a->iaoutidx;
a->forfftin[i] = a->window[i] * a->inaccum[j]; if (n1 > fsize) n1 = fsize;
a->iaoutidx = (a->iaoutidx + a->incr) % a->iasize; for (i = 0; i < n1; i++)
a->nsamps -= a->incr; 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); fftw_execute (a->Rfor);
calc_gain(a); calc_gain(a);
for (i = 0; i < a->msize; i++) for (i = 0; i < a->msize; i++)
@@ -1225,29 +1252,40 @@ void xemnr (EMNR a, int pos)
} }
post2(a); post2(a);
fftw_execute (a->Rrev); 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]; 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; const double* WDSP_RESTRICT sv;
sbegin = a->incr * (a->ovrlp - i); double* WDSP_RESTRICT oa = a->outaccum;
for (j = sbegin, k = a->oainidx; j < a->incr + sbegin; j++, k = (k + 1) % a->oasize) 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) for (j = 0; j < m1; j++) oa[k + j] = sv[j];
a->outaccum[k] = a->save[sbuff][j]; for (; j < incr; j++) oa[j - m1] = sv[j];
else }
a->outaccum[k] += a->save[sbuff][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; if (++a->saveidx == ovrlp) a->saveidx = 0;
a->oainidx = (a->oainidx + a->incr) % a->oasize; 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->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) else if (a->out != a->in)
memcpy (a->out, a->in, a->bsize * sizeof (complex)); memcpy (a->out, a->in, a->bsize * sizeof (complex));
emnr.h
+1
View File
@@ -185,6 +185,7 @@ typedef struct _emnr
double zetaThresh; double zetaThresh;
double psi; double psi;
double* nmask; double* nmask;
double* csum; // prefix sums of mask[], msize + 1 entries
double t2; double t2;
} ae; } ae;
struct _post2 struct _post2
fastmath.h
+74
View File
@@ -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
iir.c
+22 -10
View File
@@ -616,19 +616,31 @@ void xphrot (PHROT a)
if (a->run) if (a->run)
{ {
int i, n; 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]; double v = in[2 * i + 0];
for (n = 0; n < a->nstages; n++) for (n = 0; n < nstages; n++)
{ {
if (n > 0) a->x0[n] = a->y0[n - 1]; double y = b0 * v + b1 * x1[n] - a1 * y1[n];
a->y0[n] = a->b0 * a->x0[n] x1[n] = v;
+ a->b1 * a->x1[n] y1[n] = y;
- a->a1 * a->y1[n]; v = y;
a->y1[n] = a->y0[n];
a->x1[n] = a->x0[n];
} }
a->out[2 * i + 0] = a->y0[a->nstages - 1]; out[2 * i + 0] = v;
} }
} }
else if (a->out != a->in) else if (a->out != a->in)
linux_port.c
+11 -3
View File
@@ -39,9 +39,17 @@ john.d.melton@googlemail.com
#if defined(linux) || defined(__APPLE__) #if defined(linux) || defined(__APPLE__)
void QueueUserWorkItem(void *function,void *context,int flags) { void QueueUserWorkItem(void *function,void *context,int flags) {
pthread_t t; //
pthread_create(&t, NULL, function, context); // The Windows call queues the work item on a thread pool and returns at
pthread_join(t, NULL); // 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) { static inline void init_crit_section(pthread_mutex_t *mutex) {
resample.c
+117 -41
View File
@@ -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) void calc_resample (RESAMPLE a)
{ {
int x, y, z; int x, y, z;
@@ -71,7 +102,8 @@ void calc_resample (RESAMPLE a)
for (k = 0; k < a->ncoef; k += a->L) for (k = 0; k < a->ncoef; k += a->L)
a->h[i++] = impulse[j + k]; a->h[i++] = impulse[j + k];
a->ringsize = a->cpp; 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->idx_in = a->ringsize - 1;
a->phnum = 0; a->phnum = 0;
_aligned_free(impulse); _aligned_free(impulse);
@@ -79,7 +111,8 @@ void calc_resample (RESAMPLE a)
void decalc_resample (RESAMPLE a) void decalc_resample (RESAMPLE a)
{ {
_aligned_free(a->ring); _aligned_free(a->ringQ);
_aligned_free(a->ringI);
_aligned_free(a->h); _aligned_free(a->h);
} }
@@ -112,7 +145,8 @@ void destroy_resample (RESAMPLE a)
PORT PORT
void flush_resample (RESAMPLE a) 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->idx_in = a->ringsize - 1;
a->phnum = 0; a->phnum = 0;
} }
@@ -123,40 +157,49 @@ int xresample (RESAMPLE a)
int outsamps = 0; int outsamps = 0;
if (a->run) if (a->run)
{ {
int i, j, n; int i, n1;
int idx_out;
double I, Q; 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 idx_in = a->idx_in;
int ringsize = a->ringsize; int phnum = a->phnum;
double* h = a->h;
double* ring = a->ring;
for (i = 0; i < a->size; i++) for (i = 0; i < size; i++)
{ {
ring[2 * idx_in + 0] = a->in[2 * i + 0]; ringI[idx_in] = in[2 * i + 0];
ring[2 * idx_in + 1] = a->in[2 * i + 1]; ringQ[idx_in] = in[2 * i + 1];
while (a->phnum < a->L) 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; I = 0.0;
Q = 0.0; Q = 0.0;
n = cpp * a->phnum; resample_dot (hp, ringI + idx_in, ringQ + idx_in, n1, &I, &Q);
for (j = 0; j < cpp; j++) if (n1 < cpp)
{ resample_dot (hp + n1, ringI, ringQ, cpp - n1, &I, &Q);
if ((idx_out = idx_in + j) >= ringsize) idx_out -= ringsize; out[2 * outsamps + 0] = I;
I += h[n + j] * ring[2 * idx_out + 0]; out[2 * outsamps + 1] = Q;
Q += h[n + j] * ring[2 * idx_out + 1];
}
a->out[2 * outsamps + 0] = I;
a->out[2 * outsamps + 1] = Q;
outsamps++; outsamps++;
a->phnum += a->M; phnum += M;
} }
a->phnum -= a->L; phnum -= L;
if (--idx_in < 0) idx_in = a->ringsize - 1; if (--idx_in < 0) idx_in = ringsize - 1;
} }
a->idx_in = idx_in; a->idx_in = idx_in;
a->phnum = phnum;
} }
else if (a->in != a->out) else if (a->in != a->out)
memcpy (a->out, a->in, a->size * sizeof (complex)); 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 create_resampleF ( int run, int size, float* in, float* out, int in_rate, int out_rate)
{ {
RESAMPLEF a = (RESAMPLEF) malloc0 (sizeof (resampleF)); RESAMPLEF a = (RESAMPLEF) malloc0 (sizeof (resampleF));
@@ -305,31 +366,46 @@ int xresampleF (RESAMPLEF a)
int outsamps = 0; int outsamps = 0;
if (a->run) if (a->run)
{ {
int i, j, n; int i;
int idx_out;
double I; double I;
for (i = 0; i < a->size; i++) const int cpp = a->cpp;
{ const int ringsize = a->ringsize;
a->ring[a->idx_in] = (double)a->in[i]; 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; I = 0.0;
n = a->cpp * a->phnum; resampleF_dot (hp, ring + idx_in, n1, &I);
for (j = 0; j < a->cpp; j++) if (n1 < cpp)
{ resampleF_dot (hp + n1, ring, cpp - n1, &I);
if ((idx_out = a->idx_in + j) >= a->ringsize) idx_out -= a->ringsize; out[outsamps] = (float)I;
I += a->h[n + j] * a->ring[idx_out];
}
a->out[outsamps] = (float)I;
outsamps++; outsamps++;
a->phnum += a->M; phnum += M;
} }
a->phnum -= a->L; phnum -= L;
if (--a->idx_in < 0) a->idx_in = a->ringsize - 1; if (--idx_in < 0) idx_in = ringsize - 1;
} }
a->idx_in = idx_in;
a->phnum = phnum;
} }
else if (a->in != a->out) else if (a->in != a->out)
memcpy (a->out, a->in, a->size * sizeof (float)); memcpy (a->out, a->in, a->size * sizeof (float));
resample.h
+3 -1
View File
@@ -52,7 +52,9 @@ typedef struct _resample
int M; // decimation factor int M; // decimation factor
double* h; // coefficients double* h; // coefficients
int ringsize; // number of complex pairs the ring buffer holds 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 cpp; // coefficients of the phase
int phnum; // phase number int phnum; // phase number
} resample, *RESAMPLE; } resample, *RESAMPLE;
third_party/libspecbleach/libspecbleach.a
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/processors/adaptivedenoiser/adaptive_denoiser.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/processors/specbleach_adenoiser.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/processors/specbleach_denoiser.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/gain_estimation/gain_estimators.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/noise_estimation/adaptive_noise_estimator.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/post_estimation/noise_floor_manager.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/post_estimation/postfilter.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/post_estimation/spectral_whitening.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/pre_estimation/absolute_hearing_thresholds.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/pre_estimation/critical_bands.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/pre_estimation/masking_estimator.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/pre_estimation/noise_scaling_criterias.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/pre_estimation/spectral_smoother.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/pre_estimation/transient_detector.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/stft/fft_transform.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/stft/stft_buffer.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/stft/stft_processor.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/stft/stft_windows.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/utils/denoise_mixer.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/utils/general_utils.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/utils/spectral_features.o
Vendored
BIN
View File
Binary file not shown.
third_party/libspecbleach/src/shared/utils/spectral_utils.o
Vendored
BIN
View File
Binary file not shown.
third_party/rnnoise/librnnoise.a
Vendored
BIN
View File
Binary file not shown.
third_party/rnnoise/src/celt_lpc.o
Vendored
BIN
View File
Binary file not shown.
third_party/rnnoise/src/denoise.o
Vendored
BIN
View File
Binary file not shown.
third_party/rnnoise/src/kiss_fft.o
Vendored
BIN
View File
Binary file not shown.
third_party/rnnoise/src/nnet.o
Vendored
BIN
View File
Binary file not shown.
third_party/rnnoise/src/nnet_default.o
Vendored
BIN
View File
Binary file not shown.
third_party/rnnoise/src/parse_lpcnet_weights.o
Vendored
BIN
View File
Binary file not shown.
third_party/rnnoise/src/pitch.o
Vendored
BIN
View File
Binary file not shown.
third_party/rnnoise/src/rnn.o
Vendored
BIN
View File
Binary file not shown.
third_party/rnnoise/src/rnnoise_data.o
Vendored
BIN
View File
Binary file not shown.
third_party/rnnoise/src/rnnoise_data_1.o
Vendored
BIN
View File
Binary file not shown.
third_party/rnnoise/src/rnnoise_data_2.o
Vendored
BIN
View File
Binary file not shown.
third_party/rnnoise/src/rnnoise_data_3.o
Vendored
BIN
View File
Binary file not shown.
third_party/rnnoise/src/rnnoise_data_4.o
Vendored
BIN
View File
Binary file not shown.
third_party/rnnoise/src/rnnoise_data_5.o
Vendored
BIN
View File
Binary file not shown.
third_party/rnnoise/src/rnnoise_data_6.o
Vendored
BIN
View File
Binary file not shown.
third_party/rnnoise/src/rnnoise_tables.o
Vendored
BIN
View File
Binary file not shown.
varsamp.c
+88 -23
View File
@@ -59,10 +59,27 @@ void calc_varsamp (VARSAMP a)
fc_norm_low = a->fc_low / norm_rate; fc_norm_low = a->fc_low / norm_rate;
a->rsize = (int)(140.0 * norm_rate / min_rate); a->rsize = (int)(140.0 * norm_rate / min_rate);
a->ncoef = a->rsize + 1; a->ncoef = a->rsize + 1;
a->ncoef += (a->R - 1) * (a->ncoef - 1); a->ncoef += (a->R - 1) * (a->ncoef - 1); // = R * rsize + 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); /* Store the coefficients transposed into phases. hshift() walks
a->ring = (double *)malloc0(a->rsize * sizeof(complex)); 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->idx_in = a->rsize - 1;
a->h_offset = 0.0; a->h_offset = 0.0;
a->hs = (double *)malloc0 (a->rsize * sizeof (double)); a->hs = (double *)malloc0 (a->rsize * sizeof (double));
@@ -72,8 +89,9 @@ void calc_varsamp (VARSAMP a)
void decalc_varsamp (VARSAMP a) void decalc_varsamp (VARSAMP a)
{ {
_aligned_free (a->hs); _aligned_free (a->hs);
_aligned_free (a->ring); _aligned_free (a->ringQ);
_aligned_free (a->h); _aligned_free (a->ringI);
_aligned_free (a->hp);
} }
VARSAMP create_varsamp ( int run, int size, double* in, double* out, 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) 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->idx_in = a->rsize - 1;
a->h_offset = 0.0; a->h_offset = 0.0;
a->isamps = 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) void hshift (VARSAMP a)
{ {
int i, j, k; int m;
int hidx; int hidx;
double frac, pos; 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; pos = (double)a->R * a->h_offset;
hidx = (int)(pos); hidx = (int)(pos);
frac = pos - (double)hidx; frac = pos - (double)hidx;
for (i = a->rsize - 1, j = hidx, k = hidx + 1; i >= 0; i--, j += a->R, k += a->R) h0 = a->hp + (size_t)hidx * rsize;
a->hs[i] = a->h[j] + frac * (a->h[k] - a->h[j]); 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) int xvarsamp (VARSAMP a, double var)
@@ -140,13 +196,21 @@ int xvarsamp (VARSAMP a, double var)
else a->dicvar = 0.0; else a->dicvar = 0.0;
if (a->run) if (a->run)
{ {
int i, j; int i, n1;
int idx_out;
double I, Q; 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++) for (i = 0; i < a->size; i++)
{ {
a->ring[2 * a->idx_in + 0] = a->in[2 * i + 0]; ringI[idx_in] = in[2 * i + 0];
a->ring[2 * a->idx_in + 1] = a->in[2 * i + 1]; ringQ[idx_in] = in[2 * i + 1];
a->inv_cvar += a->dicvar; a->inv_cvar += a->dicvar;
picvar = (uint64_t*)(&a->inv_cvar); picvar = (uint64_t*)(&a->inv_cvar);
N = *picvar & 0xffffffffffff0000; N = *picvar & 0xffffffffffff0000;
@@ -160,20 +224,21 @@ int xvarsamp (VARSAMP a, double var)
a->h_offset += a->delta; a->h_offset += a->delta;
while (a->h_offset >= 1.0) a->h_offset -= 1.0; while (a->h_offset >= 1.0) a->h_offset -= 1.0;
while (a->h_offset < 0.0) a->h_offset += 1.0; while (a->h_offset < 0.0) a->h_offset += 1.0;
for (j = 0; j < a->rsize; j++) /* the ring wraps at most once over rsize taps; split it so both
{ halves are unit-stride */
if ((idx_out = a->idx_in + j) >= a->rsize) idx_out -= a->rsize; n1 = rsize - idx_in;
I += a->hs[j] * a->ring[2 * idx_out + 0]; varsamp_dot (a->hs, ringI + idx_in, ringQ + idx_in, n1, &I, &Q);
Q += a->hs[j] * a->ring[2 * idx_out + 1]; if (n1 < rsize)
} varsamp_dot (a->hs + n1, ringI, ringQ, rsize - n1, &I, &Q);
a->out[2 * outsamps + 0] = I; out[2 * outsamps + 0] = I;
a->out[2 * outsamps + 1] = Q; out[2 * outsamps + 1] = Q;
outsamps++; outsamps++;
a->isamps += a->inv_cvar; a->isamps += a->inv_cvar;
} }
a->isamps -= 1.0; 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) else if (a->in != a->out)
memcpy (a->out, a->in, a->size * sizeof (complex)); memcpy (a->out, a->in, a->size * sizeof (complex));
varsamp.h
+6 -2
View File
@@ -41,9 +41,13 @@ typedef struct _varsamp
double gain; double gain;
int idx_in; int idx_in;
int ncoef; 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; 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; double var;
int varmode; int varmode;
double cvar; double cvar;
wdsp.h
+25
View File
@@ -277,6 +277,7 @@ extern void SetPSHWPeak (int channel, double peak);
extern void GetPSHWPeak (int channel, double* peak); extern void GetPSHWPeak (int channel, double* peak);
extern void GetPSMaxTX (int channel, double* maxtx); extern void GetPSMaxTX (int channel, double* maxtx);
extern void SetPSPtol (int channel, double ptol); 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 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 SetPSFeedbackRate (int channel, int rate);
extern void SetPSPinMode (int channel, int pin); 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 SetTXAFMMP (int channel, int mp);
extern void SetTXAFMAFFreqs (int channel, double low, double high); 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 // Interfaces from fmsq.c
// //
wfmd.c
+322
View File
@@ -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);
}
wfmd.h
+103
View File
@@ -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
wfmmod.c
+250
View File
@@ -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);
}
wfmmod.h
+86
View File
@@ -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