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

Compare commits

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

13 Commits
4f44118299 ... wfm

Author SHA1 Message Date
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 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
30 changed files with 1778 additions and 246 deletions
Expand all files Collapse all files
.gitignore
+2
View File
@@ -9,6 +9,8 @@ java/build/
third_party/fftw/
third_party/**/*.o
third_party/**/*.a
third_party/fftw-3.3.11.tar.gz
third_party/fftw-3.3.11/
# macOS
.DS_Store
Makefile
+2
View File
@@ -144,6 +144,8 @@ utilities.c \
varsamp.c \
version.c \
wcpAGC.c \
wfmd.c \
wfmmod.c \
wisdom.c \
zetaHat.c
Makefile.android
+28 -3
View File
@@ -8,13 +8,36 @@
# FFTW source is downloaded and extracted automatically.
# Override FFTW_SRC to use an existing directory.
ANDROID_NDK ?= /home/vladimir/Android/Sdk/ndk/29.0.14206865
# Auto-detect NDK path (can always be overridden by setting ANDROID_NDK explicitly):
# 1. ANDROID_NDK_HOME environment variable
# 2. Latest NDK under ANDROID_SDK_ROOT or ANDROID_HOME
# 3. Default SDK locations: ~/Android/Sdk (Linux) or ~/Library/Android/sdk (macOS)
ifeq ($(origin ANDROID_NDK),undefined)
ifdef ANDROID_NDK_HOME
ANDROID_NDK := $(ANDROID_NDK_HOME)
else
_SDK_ROOT := $(or $(ANDROID_SDK_ROOT),$(ANDROID_HOME),\
$(wildcard $(HOME)/Android/Sdk),\
$(wildcard $(HOME)/Library/Android/sdk))
ifdef _SDK_ROOT
ANDROID_NDK := $(lastword $(sort $(wildcard $(_SDK_ROOT)/ndk/*)))
endif
endif
endif
ifndef ANDROID_NDK
$(error Cannot find Android NDK. Set ANDROID_NDK, ANDROID_NDK_HOME, ANDROID_SDK_ROOT, or ANDROID_HOME)
endif
ANDROID_API ?= 24
ANDROID_ABIS ?= arm64-v8a armeabi-v7a x86_64
ANDROID_HOST_TAG ?= $(shell uname -s | tr '[:upper:]' '[:lower:]' | sed 's/darwin/darwin/;s/linux/linux/')-$(shell uname -m | sed 's/aarch64/arm64/;s/x86_64/x86_64/')
JBR_BIN ?= /opt/android-studio/jbr/bin
JAVAC ?= $(JBR_BIN)/javac
# Auto-detect javac: JAVA_HOME > system PATH > Android Studio JBR (Linux/macOS)
_JAVAC_CANDIDATES := \
$(if $(JAVA_HOME),$(wildcard $(JAVA_HOME)/bin/javac)) \
$(shell command -v javac 2>/dev/null) \
$(wildcard /opt/android-studio/jbr/bin/javac) \
$(wildcard /Applications/Android\ Studio.app/Contents/jbr/Contents/Home/bin/javac)
JAVAC ?= $(firstword $(_JAVAC_CANDIDATES))
TOOLCHAIN := $(ANDROID_NDK)/toolchains/llvm/prebuilt/$(ANDROID_HOST_TAG)
@@ -96,6 +119,8 @@ utilities.c \
varsamp.c \
version.c \
wcpAGC.c \
wfmd.c \
wfmmod.c \
wisdom.c \
zetaHat.c
Makefile.windows
+1 -1
View File
@@ -86,7 +86,7 @@ SOURCES = \
matchedCW.c meter.c meterlog10.c nbp.c nob.c nobII.c osctrl.c \
patchpanel.c resample.c rmatch.c rnnr.c RXA.c sbnr.c sender.c \
shift.c siphon.c slew.c snb.c ssql.c syncbuffs.c TXA.c \
utilities.c varsamp.c version.c wcpAGC.c wisdom.c zetaHat.c
utilities.c varsamp.c version.c wcpAGC.c wfmd.c wfmmod.c wisdom.c zetaHat.c
OBJS := $(addprefix $(OBJDIR)/, $(SOURCES:.c=.o))
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
0); // minimum phase flag
// WFM demod
rxa[channel].wfmd.p = create_wfmd (
0, // run
ch[channel].dsp_size, // buffer size
rxa[channel].midbuff, // pointer to input buffer
rxa[channel].midbuff, // pointer to output buffer
ch[channel].dsp_rate, // sample rate
75000.0, // deviation
20.0, // f_low
15000.0, // f_high
0.02, // tau - for dc removal
1, // run de-emphasis
75.0e-6, // de-emphasis time constant
0.5, // audio gain
max(2048, ch[channel].dsp_size), // # coefs for audio cutoff filter
0); // min phase flag for audio cutoff filter
// snba
rxa[channel].snba.p = create_snba (
0, // run
@@ -579,6 +596,7 @@ void destroy_rxa (int channel)
destroy_anf (rxa[channel].anf.p);
destroy_eqp (rxa[channel].eqp.p);
destroy_snba (rxa[channel].snba.p);
destroy_wfmd (rxa[channel].wfmd.p);
destroy_fmsq (rxa[channel].fmsq.p);
destroy_fmd (rxa[channel].fmd.p);
destroy_amd (rxa[channel].amd.p);
@@ -614,6 +632,7 @@ void flush_rxa (int channel)
flush_amd (rxa[channel].amd.p);
flush_fmd (rxa[channel].fmd.p);
flush_fmsq (rxa[channel].fmsq.p);
flush_wfmd (rxa[channel].wfmd.p);
flush_snba (rxa[channel].snba.p);
flush_eqp (rxa[channel].eqp.p);
flush_anf (rxa[channel].anf.p);
@@ -649,6 +668,7 @@ void xrxa (int channel)
xamd (rxa[channel].amd.p);
xfmd (rxa[channel].fmd.p);
xfmsq (rxa[channel].fmsq.p);
xwfmd (rxa[channel].wfmd.p);
xbpsnbain (rxa[channel].bpsnba.p, 1);
xbpsnbaout (rxa[channel].bpsnba.p, 1);
xsnba (rxa[channel].snba.p);
@@ -733,6 +753,7 @@ void setDSPSamplerate_rxa (int channel)
setSamplerate_fmd (rxa[channel].fmd.p, ch[channel].dsp_rate);
setBuffers_fmsq (rxa[channel].fmsq.p, rxa[channel].midbuff, rxa[channel].midbuff, rxa[channel].fmd.p->audio);
setSamplerate_fmsq (rxa[channel].fmsq.p, ch[channel].dsp_rate);
setSamplerate_wfmd (rxa[channel].wfmd.p, ch[channel].dsp_rate);
setSamplerate_snba (rxa[channel].snba.p, ch[channel].dsp_rate);
setSamplerate_eqp (rxa[channel].eqp.p, ch[channel].dsp_rate);
setSamplerate_anf (rxa[channel].anf.p, ch[channel].dsp_rate);
@@ -794,6 +815,8 @@ void setDSPBuffsize_rxa (int channel)
setSize_fmd (rxa[channel].fmd.p, ch[channel].dsp_size);
setBuffers_fmsq (rxa[channel].fmsq.p, rxa[channel].midbuff, rxa[channel].midbuff, rxa[channel].fmd.p->audio);
setSize_fmsq (rxa[channel].fmsq.p, ch[channel].dsp_size);
setBuffers_wfmd (rxa[channel].wfmd.p, rxa[channel].midbuff, rxa[channel].midbuff);
setSize_wfmd (rxa[channel].wfmd.p, ch[channel].dsp_size);
setBuffers_snba (rxa[channel].snba.p, rxa[channel].midbuff, rxa[channel].midbuff);
setSize_snba (rxa[channel].snba.p, ch[channel].dsp_size);
setBuffers_eqp (rxa[channel].eqp.p, rxa[channel].midbuff, rxa[channel].midbuff);
@@ -857,6 +880,7 @@ void SetRXAMode (int channel, int mode)
rxa[channel].mode = mode;
rxa[channel].amd.p->run = 0;
rxa[channel].fmd.p->run = 0;
rxa[channel].wfmd.p->run = 0;
rxa[channel].agc.p->run = 1;
switch (mode)
{
@@ -875,6 +899,10 @@ void SetRXAMode (int channel, int mode)
rxa[channel].fmd.p->run = 1;
rxa[channel].agc.p->run = 0;
break;
case RXA_WFM:
rxa[channel].wfmd.p->run = 1;
rxa[channel].agc.p->run = 0;
break;
default:
break;
@@ -961,6 +989,7 @@ void RXAbpsnbaCheck (int channel, int mode, int notch_run)
run_notches = 0;
break;
case RXA_FM:
case RXA_WFM:
f_low = +a->abs_low_freq;
f_high = +a->abs_high_freq;
run_notches = 0;
@@ -1010,6 +1039,7 @@ void RXAbpsnbaSet (int channel)
a->position = 1;
break;
case RXA_FM:
case RXA_WFM:
a->run = rxa[channel].snba.p->run;
a->position = 1;
break;
@@ -1046,6 +1076,7 @@ void RXASetNC (int channel, int nc)
SetRXAFMSQNC (channel, nc);
SetRXAFMNCde (channel, nc);
SetRXAFMNCaud (channel, nc);
SetRXAWFMNCaud (channel, nc);
SetChannelState (channel, oldstate, 0);
}
@@ -1059,4 +1090,5 @@ void RXASetMP (int channel, int mp)
SetRXAFMSQMP (channel, mp);
SetRXAFMMPde (channel, mp);
SetRXAFMMPaud (channel, mp);
SetRXAWFMMPaud (channel, mp);
}
RXA.h
+6 -1
View File
@@ -41,7 +41,8 @@ enum rxaMode
RXA_SPEC,
RXA_DIGL,
RXA_SAM,
RXA_DRM
RXA_DRM,
RXA_WFM
};
enum rxaMeterType
@@ -121,6 +122,10 @@ struct _rxa
FMSQ p;
} fmsq;
struct
{
WFMD p;
} wfmd;
struct
{
EQP p;
} eqp;
TXA.c
+30
View File
@@ -358,6 +358,21 @@ void create_txa (int channel)
max(2048, ch[channel].dsp_size), // number coefficients for bandpass filter
0); // minimum phase flag
txa[channel].wfmmod.p = create_wfmmod (
0, // run - OFF by default
ch[channel].dsp_size, // size
txa[channel].midbuff, // pointer to input buffer
txa[channel].midbuff, // pointer to output buffer
ch[channel].dsp_rate, // samplerate
75000.0, // deviation
20.0, // low cutoff frequency
15000.0, // high cutoff frequency
1, // run pre-emphasis
75.0e-6, // pre-emphasis time constant
1, // run bandpass filter
max(2048, ch[channel].dsp_size), // number coefficients for bandpass filter
0); // minimum phase flag
txa[channel].gen1.p = create_gen (
0, // run
ch[channel].dsp_size, // buffer size
@@ -490,6 +505,7 @@ void destroy_txa (int channel)
destroy_meter (txa[channel].alcmeter.p);
destroy_uslew (txa[channel].uslew.p);
destroy_gen (txa[channel].gen1.p);
destroy_wfmmod (txa[channel].wfmmod.p);
destroy_fmmod (txa[channel].fmmod.p);
destroy_ammod (txa[channel].ammod.p);
destroy_wcpagc (txa[channel].alc.p);
@@ -544,6 +560,7 @@ void flush_txa (int channel)
flush_wcpagc (txa[channel].alc.p);
flush_ammod (txa[channel].ammod.p);
flush_fmmod (txa[channel].fmmod.p);
flush_wfmmod (txa[channel].wfmmod.p);
flush_gen (txa[channel].gen1.p);
flush_uslew (txa[channel].uslew.p);
flush_meter (txa[channel].alcmeter.p);
@@ -580,6 +597,7 @@ void xtxa (int channel)
xammod (txa[channel].ammod.p); // AM Modulator
xemphp (txa[channel].preemph.p, 1); // FM pre-emphasis (second option)
xfmmod (txa[channel].fmmod.p); // FM Modulator
xwfmmod (txa[channel].wfmmod.p); // WFM Modulator (pre-emphasis is internal)
xgen (txa[channel].gen1.p); // output signal generator (TUN and Two-tone)
xuslew (txa[channel].uslew.p); // up-slew for AM, FM, and gens
xmeter (txa[channel].alcmeter.p); // ALC Meter
@@ -653,6 +671,7 @@ void setDSPSamplerate_txa (int channel)
setSamplerate_wcpagc (txa[channel].alc.p, ch[channel].dsp_rate);
setSamplerate_ammod (txa[channel].ammod.p, ch[channel].dsp_rate);
setSamplerate_fmmod (txa[channel].fmmod.p, ch[channel].dsp_rate);
setSamplerate_wfmmod (txa[channel].wfmmod.p, ch[channel].dsp_rate);
setSamplerate_gen (txa[channel].gen1.p, ch[channel].dsp_rate);
setSamplerate_uslew (txa[channel].uslew.p, ch[channel].dsp_rate);
setSamplerate_meter (txa[channel].alcmeter.p, ch[channel].dsp_rate);
@@ -723,6 +742,8 @@ void setDSPBuffsize_txa (int channel)
setSize_ammod (txa[channel].ammod.p, ch[channel].dsp_size);
setBuffers_fmmod (txa[channel].fmmod.p, txa[channel].midbuff, txa[channel].midbuff);
setSize_fmmod (txa[channel].fmmod.p, ch[channel].dsp_size);
setBuffers_wfmmod (txa[channel].wfmmod.p, txa[channel].midbuff, txa[channel].midbuff);
setSize_wfmmod (txa[channel].wfmmod.p, ch[channel].dsp_size);
setBuffers_gen (txa[channel].gen1.p, txa[channel].midbuff, txa[channel].midbuff);
setSize_gen (txa[channel].gen1.p, ch[channel].dsp_size);
setBuffers_uslew (txa[channel].uslew.p, txa[channel].midbuff, txa[channel].midbuff);
@@ -758,6 +779,7 @@ void SetTXAMode (int channel, int mode)
txa[channel].mode = mode;
txa[channel].ammod.p->run = 0;
txa[channel].fmmod.p->run = 0;
txa[channel].wfmmod.p->run = 0;
txa[channel].preemph.p->run = 0;
switch (mode)
{
@@ -779,6 +801,10 @@ void SetTXAMode (int channel, int mode)
txa[channel].fmmod.p->run = 1;
txa[channel].preemph.p->run = 1;
break;
case TXA_WFM:
// wfmmod carries its own RC pre-emphasis; the shared emphp stays off
txa[channel].wfmmod.p->run = 1;
break;
default:
break;
@@ -820,6 +846,7 @@ int TXAUslewCheck (int channel)
{
return (txa[channel].ammod.p->run == 1) ||
(txa[channel].fmmod.p->run == 1) ||
(txa[channel].wfmmod.p->run == 1) ||
(txa[channel].gen0.p->run == 1) ||
(txa[channel].gen1.p->run == 1);
}
@@ -855,6 +882,7 @@ void TXASetupBPFilters (int channel)
case TXA_AM:
case TXA_SAM:
case TXA_FM:
case TXA_WFM:
if (txa[channel].compressor.p->run)
{
CalcBandpassFilter (txa[channel].bp0.p, 0.0, txa[channel].f_high, 2.0);
@@ -914,6 +942,7 @@ void TXASetNC (int channel, int nc)
SetTXAFMEmphNC (channel, nc);
SetTXAEQNC (channel, nc);
SetTXAFMNC (channel, nc);
SetTXAWFMNC (channel, nc);
SetTXACFIRNC (channel, nc);
SetChannelState (channel, oldstate, 0);
}
@@ -925,6 +954,7 @@ void TXASetMP (int channel, int mp)
SetTXAFMEmphMP (channel, mp);
SetTXAEQMP (channel, mp);
SetTXAFMMP (channel, mp);
SetTXAWFMMP (channel, mp);
}
PORT
TXA.h
+6 -1
View File
@@ -43,7 +43,8 @@ enum txaMode
TXA_SAM,
TXA_DRM,
TXA_AM_LSB,
TXA_AM_USB
TXA_AM_USB,
TXA_WFM
};
enum txaMeterType
@@ -135,6 +136,10 @@ struct _txa
FMMOD p;
} fmmod;
struct
{
WFMMOD p;
} wfmmod;
struct
{
SIPHON p;
} sip1;
analyzer.c
+19 -1
View File
@@ -1087,7 +1087,15 @@ void __cdecl sendbuf(void *arg)
LeaveCriticalSection(&(a->BufferControlSection[a->ss][a->LO]));
}
}
Sleep(1);
//
// Block until a Spectrum*() call announces new samples. This used to be
// Sleep(1), i.e. 1000 wakeups per second spent re-reading the same flags
// -- 0.67% of a core even with no data arriving at all.
//
// Whoever sets end_dispatcher also signals the semaphore, so the wait
// below always has a way out.
//
WaitForSingleObject(a->Sem_BuffReady, INFINITE);
}
InterlockedBitTestAndReset(&a->dispatcher, 0);
_endthread();
@@ -1202,6 +1210,9 @@ void SetAnalyzer ( int disp, // display identifier
EnterCriticalSection(&a->SetAnalyzerSection);
a->end_dispatcher = 1;
// wake the dispatcher out of its blocking wait so it can observe the flag;
// it does not take SetAnalyzerSection, so holding it here is safe
ReleaseSemaphore(a->Sem_BuffReady, 1, 0);
while (InterlockedAnd(&a->dispatcher, 1))
Sleep(1);
a->stop = 1;
@@ -1347,6 +1358,7 @@ void XCreateAnalyzer( int disp,
a->hSnapEvent[i][j] = CreateEvent(NULL, FALSE, FALSE, TEXT("snap"));
a->snap[i][j] = 0;
}
a->Sem_BuffReady = CreateSemaphore(0, 0, 1000, 0);
InitializeCriticalSectionAndSpinCount(&a->ResampleSection, 0);
InitializeCriticalSectionAndSpinCount(&a->SetAnalyzerSection, 0);
InitializeCriticalSectionAndSpinCount(&a->StitchSection, 0);
@@ -1434,6 +1446,7 @@ void DestroyAnalyzer(int disp)
int i, j;
a->end_dispatcher = 1;
ReleaseSemaphore(a->Sem_BuffReady, 1, 0);
while (InterlockedAnd(&a->dispatcher, 1))
Sleep(1);
@@ -1497,6 +1510,7 @@ void DestroyAnalyzer(int disp)
for (i = 0; i < a->max_stitch; i++)
for (j = 0; j < a->max_num_fft; j++)
CloseHandle(a->hSnapEvent[i][j]);
CloseHandle(a->Sem_BuffReady);
_aligned_free ((void *) a->pnum_threads);
@@ -1633,6 +1647,7 @@ void CloseBuffer(int disp, int ss, int LO)
if((a->IQin_index[ss][LO] += a->buff_size) >= a->bsize) //REQUIRES buff_size IS A SUB-MULTIPLE OF SIZE OF INPUT SAMPLE BUFFS!
a->IQin_index[ss][LO] = 0;
ReleaseSemaphore(a->Sem_BuffReady, 1, 0); // new samples: let the dispatcher run
if (!InterlockedAnd(&a->dispatcher, 1))
{
InterlockedBitTestAndSet (&a->dispatcher, 0);
@@ -1672,6 +1687,7 @@ void Spectrum(int disp, int ss, int LO, dINREAL* pI, dINREAL* pQ)
if((a->IQin_index[ss][LO] += a->buff_size) >= a->bsize) //REQUIRES buff_size IS A SUB-MULTIPLE OF SIZE OF INPUT SAMPLE BUFFS!
a->IQin_index[ss][LO] = 0;
ReleaseSemaphore(a->Sem_BuffReady, 1, 0); // new samples: let the dispatcher run
if (!InterlockedAnd(&a->dispatcher, 1))
{
InterlockedBitTestAndSet(&a->dispatcher, 0);
@@ -1717,6 +1733,7 @@ void Spectrum2(int run, int disp, int ss, int LO, dINREAL* pbuff)
if((a->IQin_index[ss][LO] += a->buff_size) >= a->bsize) //REQUIRES buff_size IS A SUB-MULTIPLE OF SIZE OF INPUT SAMPLE BUFFS!
a->IQin_index[ss][LO] = 0;
ReleaseSemaphore(a->Sem_BuffReady, 1, 0); // new samples: let the dispatcher run
if (!InterlockedAnd(&a->dispatcher, 1))
{
InterlockedBitTestAndSet(&a->dispatcher, 0);
@@ -1763,6 +1780,7 @@ void Spectrum0(int run, int disp, int ss, int LO, double* pbuff)
if((a->IQin_index[ss][LO] += a->buff_size) >= a->bsize) //REQUIRES buff_size IS A SUB-MULTIPLE OF SIZE OF INPUT SAMPLE BUFFS!
a->IQin_index[ss][LO] = 0;
ReleaseSemaphore(a->Sem_BuffReady, 1, 0); // new samples: let the dispatcher run
if (!InterlockedAnd(&a->dispatcher, 1))
{
InterlockedBitTestAndSet(&a->dispatcher, 0);
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
double *snap_buff[dMAX_STITCH][dMAX_NUM_FFT]; // pointers to buffers for the snap
HANDLE Sem_BuffReady; // signalled when input samples arrive, so the
// dispatcher can block instead of polling
CRITICAL_SECTION PB_ControlsSection[dMAX_PIXOUTS];
CRITICAL_SECTION SetAnalyzerSection;
CRITICAL_SECTION BufferControlSection[dMAX_STITCH][dMAX_NUM_FFT];
anf.c
+95 -28
View File
@@ -26,6 +26,45 @@ warren@wpratt.com
#include "comm.h"
/* Filter output and tap-window energy over a unit-stride run of the delay line.
The delay line is indexed (in_idx + j + delay) & mask, which wraps at most
once across the tap window; xanf() splits the window at the wrap so both
halves are contiguous here. Four independent accumulator pairs keep the FMAs
off a single dependency chain and let the vectorizer in -- a 'y += w[j]*x[j]'
reduction cannot be reassociated without -ffast-math, which this library must
not enable (it relies on IEEE semantics for 0/0 = NaN and x/0 = Inf). */
static inline void anf_dot (const double* WDSP_RESTRICT w,
const double* WDSP_RESTRICT x, int n, double* py, double* psigma)
{
double y0 = 0.0, y1 = 0.0, y2 = 0.0, y3 = 0.0;
double s0 = 0.0, s1 = 0.0, s2 = 0.0, s3 = 0.0;
int j = 0;
for (; j <= n - 4; j += 4)
{
y0 += w[j + 0] * x[j + 0]; s0 += x[j + 0] * x[j + 0];
y1 += w[j + 1] * x[j + 1]; s1 += x[j + 1] * x[j + 1];
y2 += w[j + 2] * x[j + 2]; s2 += x[j + 2] * x[j + 2];
y3 += w[j + 3] * x[j + 3]; s3 += x[j + 3] * x[j + 3];
}
for (; j < n; j++)
{
y0 += w[j] * x[j];
s0 += x[j] * x[j];
}
*py += (y0 + y1) + (y2 + y3);
*psigma += (s0 + s1) + (s2 + s3);
}
/* Leaky-LMS tap update over the same unit-stride run. */
static inline void anf_update (double* WDSP_RESTRICT w,
const double* WDSP_RESTRICT x, int n, double c0, double c1)
{
int j;
for (j = 0; j < n; j++)
w[j] = c0 * w[j] + c1 * x[j];
}
ANF create_anf (
int run,
int position,
@@ -81,53 +120,81 @@ void destroy_anf (ANF a)
void xanf(ANF a, int position)
{
int i, j, idx;
int i;
double c0, c1;
double y, error, sigma, inv_sigp;
double nel, nev;
if (a->run && (a->position == position))
{
for (i = 0; i < a->buff_size; i++)
{
a->d[a->in_idx] = a->in_buff[2 * i + 0];
const int n_taps = a->n_taps;
const int dline_size = a->dline_size;
const int mask = a->mask;
const int delay = a->delay;
const int buff_size = a->buff_size;
const double two_mu = a->two_mu;
const double gamma = a->gamma;
const double den_mult = a->den_mult;
const double lincr = a->lincr;
const double ldecr = a->ldecr;
const double lidx_min = a->lidx_min;
const double lidx_max = a->lidx_max;
/* in_buff and out_buff are the same buffer in RXA, so neither may be
marked restrict; d and w are private to the struct. */
const double* in_buff = a->in_buff;
double* out_buff = a->out_buff;
double* WDSP_RESTRICT d = a->d;
double* WDSP_RESTRICT w = a->w;
int in_idx = a->in_idx;
double lidx = a->lidx;
double ngamma = a->ngamma;
y = 0;
sigma = 0;
for (j = 0; j < a->n_taps; j++)
for (i = 0; i < buff_size; i++)
{
idx = (a->in_idx + j + a->delay) & a->mask;
y += a->w[j] * a->d[idx];
sigma += a->d[idx] * a->d[idx];
}
double dsamp;
int base, n1;
dsamp = in_buff[2 * i + 0];
d[in_idx] = dsamp;
base = (in_idx + delay) & mask;
if ((n1 = dline_size - base) > n_taps) n1 = n_taps;
y = 0.0;
sigma = 0.0;
anf_dot (w, d + base, n1, &y, &sigma);
if (n1 < n_taps)
anf_dot (w + n1, d, n_taps - n1, &y, &sigma);
inv_sigp = 1.0 / (sigma + 1e-10);
error = a->d[a->in_idx] - y;
error = dsamp - y;
a->out_buff[2 * i + 0] = error;
a->out_buff[2 * i + 1] = 0.0;
out_buff[2 * i + 0] = error;
out_buff[2 * i + 1] = 0.0;
if((nel = error * (1.0 - a->two_mu * sigma * inv_sigp)) < 0.0) nel = -nel;
if((nev = a->d[a->in_idx] - (1.0 - a->two_mu * a->ngamma) * y - a->two_mu * error * sigma * inv_sigp) < 0.0) nev = -nev;
if((nel = error * (1.0 - two_mu * sigma * inv_sigp)) < 0.0) nel = -nel;
if((nev = dsamp - (1.0 - two_mu * ngamma) * y - two_mu * error * sigma * inv_sigp) < 0.0) nev = -nev;
if (nev < nel)
{
if ((a->lidx += a->lincr) > a->lidx_max) a->lidx = a->lidx_max;
if ((lidx += lincr) > lidx_max) lidx = lidx_max;
}
else
{
if ((a->lidx -= a->ldecr) < a->lidx_min) a->lidx = a->lidx_min;
if ((lidx -= ldecr) < lidx_min) lidx = lidx_min;
}
a->ngamma = a->gamma * (a->lidx * a->lidx) * (a->lidx * a->lidx) * a->den_mult;
ngamma = gamma * (lidx * lidx) * (lidx * lidx) * den_mult;
c0 = 1.0 - a->two_mu * a->ngamma;
c1 = a->two_mu * error * inv_sigp;
c0 = 1.0 - two_mu * ngamma;
c1 = two_mu * error * inv_sigp;
for (j = 0; j < a->n_taps; j++)
{
idx = (a->in_idx + j + a->delay) & a->mask;
a->w[j] = c0 * a->w[j] + c1 * a->d[idx];
}
a->in_idx = (a->in_idx + a->mask) & a->mask;
anf_update (w, d + base, n1, c0, c1);
if (n1 < n_taps)
anf_update (w + n1, d, n_taps - n1, c0, c1);
in_idx = (in_idx + mask) & mask;
}
a->in_idx = in_idx;
a->lidx = lidx;
a->ngamma = ngamma;
}
else if (a->in_buff != a->out_buff)
memcpy (a->out_buff, a->in_buff, a->buff_size * sizeof (complex));
anr.c
+95 -28
View File
@@ -26,6 +26,45 @@ warren@wpratt.com
#include "comm.h"
/* Filter output and tap-window energy over a unit-stride run of the delay line.
The delay line is indexed (in_idx + j + delay) & mask, which wraps at most
once across the tap window; xanr() splits the window at the wrap so both
halves are contiguous here. Four independent accumulator pairs keep the FMAs
off a single dependency chain and let the vectorizer in -- an 'y += w[j]*x[j]'
reduction cannot be reassociated without -ffast-math, which this library must
not enable (it relies on IEEE semantics for 0/0 = NaN and x/0 = Inf). */
static inline void anr_dot (const double* WDSP_RESTRICT w,
const double* WDSP_RESTRICT x, int n, double* py, double* psigma)
{
double y0 = 0.0, y1 = 0.0, y2 = 0.0, y3 = 0.0;
double s0 = 0.0, s1 = 0.0, s2 = 0.0, s3 = 0.0;
int j = 0;
for (; j <= n - 4; j += 4)
{
y0 += w[j + 0] * x[j + 0]; s0 += x[j + 0] * x[j + 0];
y1 += w[j + 1] * x[j + 1]; s1 += x[j + 1] * x[j + 1];
y2 += w[j + 2] * x[j + 2]; s2 += x[j + 2] * x[j + 2];
y3 += w[j + 3] * x[j + 3]; s3 += x[j + 3] * x[j + 3];
}
for (; j < n; j++)
{
y0 += w[j] * x[j];
s0 += x[j] * x[j];
}
*py += (y0 + y1) + (y2 + y3);
*psigma += (s0 + s1) + (s2 + s3);
}
/* Leaky-LMS tap update over the same unit-stride run. */
static inline void anr_update (double* WDSP_RESTRICT w,
const double* WDSP_RESTRICT x, int n, double c0, double c1)
{
int j;
for (j = 0; j < n; j++)
w[j] = c0 * w[j] + c1 * x[j];
}
ANR create_anr (
int run,
int position,
@@ -81,53 +120,81 @@ void destroy_anr (ANR a)
void xanr (ANR a, int position)
{
int i, j, idx;
int i;
double c0, c1;
double y, error, sigma, inv_sigp;
double nel, nev;
if (a->run && (a->position == position))
{
for (i = 0; i < a->buff_size; i++)
{
a->d[a->in_idx] = a->in_buff[2 * i + 0];
const int n_taps = a->n_taps;
const int dline_size = a->dline_size;
const int mask = a->mask;
const int delay = a->delay;
const int buff_size = a->buff_size;
const double two_mu = a->two_mu;
const double gamma = a->gamma;
const double den_mult = a->den_mult;
const double lincr = a->lincr;
const double ldecr = a->ldecr;
const double lidx_min = a->lidx_min;
const double lidx_max = a->lidx_max;
/* in_buff and out_buff are the same buffer in RXA, so neither may be
marked restrict; d and w are private to the struct. */
const double* in_buff = a->in_buff;
double* out_buff = a->out_buff;
double* WDSP_RESTRICT d = a->d;
double* WDSP_RESTRICT w = a->w;
int in_idx = a->in_idx;
double lidx = a->lidx;
double ngamma = a->ngamma;
y = 0;
sigma = 0;
for (j = 0; j < a->n_taps; j++)
for (i = 0; i < buff_size; i++)
{
idx = (a->in_idx + j + a->delay) & a->mask;
y += a->w[j] * a->d[idx];
sigma += a->d[idx] * a->d[idx];
}
double dsamp;
int base, n1;
dsamp = in_buff[2 * i + 0];
d[in_idx] = dsamp;
base = (in_idx + delay) & mask;
if ((n1 = dline_size - base) > n_taps) n1 = n_taps;
y = 0.0;
sigma = 0.0;
anr_dot (w, d + base, n1, &y, &sigma);
if (n1 < n_taps)
anr_dot (w + n1, d, n_taps - n1, &y, &sigma);
inv_sigp = 1.0 / (sigma + 1e-10);
error = a->d[a->in_idx] - y;
error = dsamp - y;
a->out_buff[2 * i + 0] = y;
a->out_buff[2 * i + 1] = 0.0;
out_buff[2 * i + 0] = y;
out_buff[2 * i + 1] = 0.0;
if((nel = error * (1.0 - a->two_mu * sigma * inv_sigp)) < 0.0) nel = -nel;
if((nev = a->d[a->in_idx] - (1.0 - a->two_mu * a->ngamma) * y - a->two_mu * error * sigma * inv_sigp) < 0.0) nev = -nev;
if((nel = error * (1.0 - two_mu * sigma * inv_sigp)) < 0.0) nel = -nel;
if((nev = dsamp - (1.0 - two_mu * ngamma) * y - two_mu * error * sigma * inv_sigp) < 0.0) nev = -nev;
if (nev < nel)
{
if ((a->lidx += a->lincr) > a->lidx_max) a->lidx = a->lidx_max;
if ((lidx += lincr) > lidx_max) lidx = lidx_max;
}
else
{
if ((a->lidx -= a->ldecr) < a->lidx_min) a->lidx = a->lidx_min;
if ((lidx -= ldecr) < lidx_min) lidx = lidx_min;
}
a->ngamma = a->gamma * (a->lidx * a->lidx) * (a->lidx * a->lidx) * a->den_mult;
ngamma = gamma * (lidx * lidx) * (lidx * lidx) * den_mult;
c0 = 1.0 - a->two_mu * a->ngamma;
c1 = a->two_mu * error * inv_sigp;
c0 = 1.0 - two_mu * ngamma;
c1 = two_mu * error * inv_sigp;
for (j = 0; j < a->n_taps; j++)
{
idx = (a->in_idx + j + a->delay) & a->mask;
a->w[j] = c0 * a->w[j] + c1 * a->d[idx];
}
a->in_idx = (a->in_idx + a->mask) & a->mask;
anr_update (w, d + base, n1, c0, c1);
if (n1 < n_taps)
anr_update (w + n1, d, n_taps - n1, c0, c1);
in_idx = (in_idx + mask) & mask;
}
a->in_idx = in_idx;
a->lidx = lidx;
a->ngamma = ngamma;
}
else if (a->in_buff != a->out_buff)
memcpy (a->out_buff, a->in_buff, a->buff_size * sizeof (complex));
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->npsamps = npsamps;
a->alpha = alpha;
a->outlier_sigma = 0.0;
a->info = (int *) malloc0 (16 * sizeof (int));
a->binfo = (int *) malloc0 (16 * sizeof (int));
@@ -321,6 +322,125 @@ void rxscheck (int rints, double* tvec, double* coef, int* info)
if (out < 0.00) *info |= 0x0020;
}
// Yurij_eu2av: fallback rx_scale estimator. It averages the top few
// amplitude intervals (ignoring overrange samples) and linearly extrapolates
// to full TX scale (env_TX = 1/hw_scale). Used only if the cubic xbuilder
// fit fails or is rejected by rxscheck.
static int estimate_rx_scale_from_top_intervals(CALCC a, double* rx_scale_out)
{
const int n_top = 4;
double sx[4], sy[4], sw[4];
int valid = 0;
int b, j;
for (b = a->ints - 1; b >= 0 && valid < n_top; b--)
{
int base = b * a->spi;
double sum_x = 0.0, sum_y = 0.0;
int n = 0;
for (j = 0; j < a->spi; j++)
{
int k = base + j;
double nx = a->env_TX[k] * a->hw_scale;
if (nx > 1.0 || nx < 0.0) continue;
if (a->env_TX[k] < 1.0e-30 || a->env_RX[k] < 1.0e-30) continue;
sum_x += a->env_TX[k];
sum_y += a->env_RX[k];
n++;
}
if (n == 0) continue;
sx[valid] = sum_x / (double)n;
sy[valid] = sum_y / (double)n;
sw[valid] = (double)n;
valid++;
}
if (valid < 2) return -1;
{
double s_w = 0.0, s_x = 0.0, s_y = 0.0, s_xx = 0.0, s_xy = 0.0;
double det, aa, bb, target_x, y_at_target;
int i;
for (i = 0; i < valid; i++)
{
double w = sw[i];
s_w += w;
s_x += w * sx[i];
s_y += w * sy[i];
s_xx += w * sx[i] * sx[i];
s_xy += w * sx[i] * sy[i];
}
det = s_w * s_xx - s_x * s_x;
if (fabs(det) < 1e-30) return -1;
bb = (s_w * s_xy - s_x * s_y) / det;
aa = (s_y - bb * s_x) / s_w;
target_x = 1.0 / a->hw_scale;
y_at_target = aa + bb * target_x;
if (y_at_target <= 1e-15) return -1;
*rx_scale_out = 1.0 / y_at_target;
}
return 0;
}
// Yurij_eu2av: robust outlier rejection for the cubic-spline xbuilder.
// Fits rx = k*tx through the origin via median ratio, then rejects points
// whose residual exceeds sigma * MAD.
static int cmp_double(const void* a, const void* b)
{
double da = *(const double*)a;
double db = *(const double*)b;
if (da < db) return -1;
if (da > db) return 1;
return 0;
}
static double median_double(double* v, int n)
{
if (n <= 0) return 0.0;
if (n % 2 == 1)
return v[n / 2];
else
return 0.5 * (v[n / 2 - 1] + v[n / 2]);
}
static int reject_outliers(double* tx, double* rx, int n, double sigma)
{
const int min_points = 32;
int i, keep = 0;
double* ratios;
double* absres;
double med_ratio, med_absres, thr;
if (n < min_points || sigma <= 0.0) return n;
ratios = (double*)malloc0(n * sizeof(double));
for (i = 0; i < n; i++)
ratios[i] = (tx[i] > 1.0e-30) ? rx[i] / tx[i] : 0.0;
qsort(ratios, n, sizeof(double), cmp_double);
med_ratio = median_double(ratios, n);
_aligned_free(ratios);
if (fabs(med_ratio) < 1.0e-30) return n;
absres = (double*)malloc0(n * sizeof(double));
for (i = 0; i < n; i++)
absres[i] = fabs(rx[i] - med_ratio * tx[i]);
qsort(absres, n, sizeof(double), cmp_double);
med_absres = median_double(absres, n);
_aligned_free(absres);
if (med_absres < 1.0e-30) return n;
thr = sigma * med_absres;
for (i = 0; i < n; i++)
{
if (fabs(rx[i] - med_ratio * tx[i]) <= thr)
{
tx[keep] = tx[i];
rx[keep] = rx[i];
keep++;
}
}
return (keep >= min_points) ? keep : n;
}
void calc (CALCC a)
{
int i;
@@ -336,21 +456,60 @@ void calc (CALCC a)
double tvec[3];
double txrxcoefs[4 * 2];
double rx_scale;
int xb_ok = 0;
double* tx_filt;
double* rx_filt;
int n_filt = 0;
if (a->ints < 16) rints = 1;
else rints = 2;
ix = rints - 1;
for (i = 0; i <= rints; i++)
tvec[i] = (double)i / (double)rints / a->hw_scale;
dx = tvec[rints] - tvec[rints - 1];
xbuilder(a->ccbld, a->nsamps, a->env_TX, a->env_RX, rints, tvec, &(a->binfo[0]), txrxcoefs, a->ptol);
// Yurij_eu2av: build a filtered dataset with overrange samples removed
// before running xbuilder. Overrange env_TX*hw_scale > 1.0 can distort
// the cubic fit and produce an incorrect rx_scale.
tx_filt = (double*)malloc0(a->nsamps * sizeof(double));
rx_filt = (double*)malloc0(a->nsamps * sizeof(double));
for (i = 0; i < a->nsamps; i++)
{
double nx = a->env_TX[i] * a->hw_scale;
if (nx > 1.0 || nx < 0.0) continue;
if (a->env_TX[i] < 1.0e-30 || a->env_RX[i] < 1.0e-30) continue;
tx_filt[n_filt] = a->env_TX[i];
rx_filt[n_filt] = a->env_RX[i];
n_filt++;
}
// Yurij_eu2av: optional outlier rejection before cubic-spline fit.
if (a->outlier_sigma > 0.0)
n_filt = reject_outliers(tx_filt, rx_filt, n_filt, a->outlier_sigma);
xbuilder(a->ccbld, n_filt, tx_filt, rx_filt, rints, tvec, &(a->binfo[0]), txrxcoefs, a->ptol);
rxscheck (rints, tvec, txrxcoefs, &a->binfo[7]);
if ((a->binfo[0] == 0) && (a->binfo[7] == 0))
{
rx_scale = 1.0 / (txrxcoefs[4 * ix + 0] + dx * (txrxcoefs[4 * ix + 1] + dx * (txrxcoefs[4 * ix + 2] + dx * txrxcoefs[4 * ix + 3])));
else
xb_ok = 1;
}
else if (estimate_rx_scale_from_top_intervals(a, &rx_scale) == 0)
{
// Yurij_eu2av: xbuilder failed, but the bucket-average fallback
// gave a usable rx_scale. Keep binfo[0] bit 0 set for diagnostics.
a->binfo[0] |= 0x0001;
xb_ok = 1;
}
_aligned_free(tx_filt);
_aligned_free(rx_filt);
if (!xb_ok)
{
a->scOK = 0;
goto cleanup;
}
if (a->stbl && _InterlockedAnd (&a->ctrl.running, 1))
a->rx_scale = a->alpha * a->rx_scale + (1.0 - a->alpha) * rx_scale;
else
@@ -1046,6 +1205,16 @@ void SetPSPtol (int channel, double ptol)
LeaveCriticalSection (&txa[channel].calcc.cs_update);
}
PORT
void SetPSOutlierSigma (int channel, double sigma)
{
// Yurij_eu2av: 0.0 disables the pre-xbuilder outlier filter.
if (sigma < 0.0) sigma = 0.0;
EnterCriticalSection (&txa[channel].calcc.cs_update);
txa[channel].calcc.p->outlier_sigma = sigma;
LeaveCriticalSection (&txa[channel].calcc.cs_update);
}
PORT
void GetPSDisp (int channel, double* x, double* ym, double* yc, double* ys, double* cm, double* cc, double* cs)
{
calcc.h
+1
View File
@@ -48,6 +48,7 @@ typedef struct _calcc
double hw_scale;
double rx_scale;
double alpha;
double outlier_sigma;
int tsamps;
double* env_TX;
cfcomp.c
+56 -27
View File
@@ -336,18 +336,36 @@ void xcfcomp (CFCOMP a, int pos)
if (a->run && pos == a->position)
{
int i, j, k, sbuff, sbegin;
for (i = 0; i < 2 * a->bsize; i += 2)
/* Each ring index below steps by one and, since iasize >= fsize and
oasize >= incr always hold, wraps at most once per loop. The '% size'
per step was therefore an integer division for nothing -- about 5100
of them per call at fsize = 2048. Walk contiguous runs instead. */
const int iasize = a->iasize;
const int oasize = a->oasize;
const int fsize = a->fsize;
const int incr = a->incr;
const int bsize = a->bsize;
const int ovrlp = a->ovrlp;
const double pregain = a->pregain;
const double postgain = a->postgain;
for (i = 0, j = a->iainidx; i < 2 * bsize; i += 2)
{
a->inaccum[a->iainidx] = a->in[i];
a->iainidx = (a->iainidx + 1) % a->iasize;
a->inaccum[j] = a->in[i];
if (++j == iasize) j = 0;
}
a->nsamps += a->bsize;
while (a->nsamps >= a->fsize)
a->iainidx = j;
a->nsamps += bsize;
while (a->nsamps >= fsize)
{
for (i = 0, j = a->iaoutidx; i < a->fsize; i++, j = (j + 1) % a->iasize)
a->forfftin[i] = a->pregain * a->window[i] * a->inaccum[j];
a->iaoutidx = (a->iaoutidx + a->incr) % a->iasize;
a->nsamps -= a->incr;
int n1 = iasize - a->iaoutidx;
if (n1 > fsize) n1 = fsize;
for (i = 0; i < n1; i++)
a->forfftin[i] = pregain * a->window[i] * a->inaccum[a->iaoutidx + i];
for (; i < fsize; i++)
a->forfftin[i] = pregain * a->window[i] * a->inaccum[i - n1];
if ((a->iaoutidx += incr) >= iasize) a->iaoutidx -= iasize;
a->nsamps -= incr;
fftw_execute (a->Rfor);
calc_mask(a);
for (i = 0; i < a->msize; i++)
@@ -356,30 +374,41 @@ void xcfcomp (CFCOMP a, int pos)
a->revfftin[2 * i + 1] = a->mask[i] * a->forfftout[2 * i + 1];
}
fftw_execute (a->Rrev);
for (i = 0; i < a->fsize; i++)
a->save[a->saveidx][i] = a->postgain * a->window[i] * a->revfftout[i];
for (i = a->ovrlp; i > 0; i--)
for (i = 0; i < fsize; i++)
a->save[a->saveidx][i] = postgain * a->window[i] * a->revfftout[i];
for (i = ovrlp; i > 0; i--)
{
sbuff = (a->saveidx + i) % a->ovrlp;
sbegin = a->incr * (a->ovrlp - i);
for (j = sbegin, k = a->oainidx; j < a->incr + sbegin; j++, k = (k + 1) % a->oasize)
const double* WDSP_RESTRICT sv;
double* WDSP_RESTRICT oa = a->outaccum;
int m1;
sbuff = (a->saveidx + i) % ovrlp;
sbegin = incr * (ovrlp - i);
sv = a->save[sbuff] + sbegin;
m1 = oasize - a->oainidx;
if (m1 > incr) m1 = incr;
k = a->oainidx;
if (i == ovrlp)
{
if ( i == a->ovrlp)
a->outaccum[k] = a->save[sbuff][j];
for (j = 0; j < m1; j++) oa[k + j] = sv[j];
for (; j < incr; j++) oa[j - m1] = sv[j];
}
else
a->outaccum[k] += a->save[sbuff][j];
}
}
a->saveidx = (a->saveidx + 1) % a->ovrlp;
a->oainidx = (a->oainidx + a->incr) % a->oasize;
}
for (i = 0; i < a->bsize; i++)
{
a->out[2 * i + 0] = a->outaccum[a->oaoutidx];
a->out[2 * i + 1] = 0.0;
a->oaoutidx = (a->oaoutidx + 1) % a->oasize;
for (j = 0; j < m1; j++) oa[k + j] += sv[j];
for (; j < incr; j++) oa[j - m1] += sv[j];
}
}
if (++a->saveidx == ovrlp) a->saveidx = 0;
if ((a->oainidx += incr) >= oasize) a->oainidx -= oasize;
}
for (i = 0, k = a->oaoutidx; i < bsize; i++)
{
a->out[2 * i + 0] = a->outaccum[k];
a->out[2 * i + 1] = 0.0;
if (++k == oasize) k = 0;
}
a->oaoutidx = k;
}
else if (a->out != a->in)
memcpy (a->out, a->in, a->bsize * sizeof (complex));
}
comm.h
+11
View File
@@ -75,12 +75,15 @@ warren@wpratt.com
#include "sbnr.h" // NR3 + NR4 support
#include "emph.h"
#include "eq.h"
#include "fastmath.h"
#include "fcurve.h"
#include "fir.h"
#include "firmin.h"
#include "fmd.h"
#include "fmmod.h"
#include "fmsq.h"
#include "wfmd.h"
#include "wfmmod.h"
#include "gain.h"
#include "gaussian.h"
#include "gen.h"
@@ -150,6 +153,14 @@ warren@wpratt.com
#define PI 3.1415926535897932
#define TWOPI 6.2831853071795864
// Non-aliasing qualifier for DSP buffers. Spelled __restrict rather than
// restrict because the JNI translation unit is compiled as -std=gnu89.
#if defined(__GNUC__) || defined(__clang__) || defined(_MSC_VER)
#define WDSP_RESTRICT __restrict
#else
#define WDSP_RESTRICT
#endif
// miscellaneous
typedef double complex[2];
#define PORT __declspec( dllexport )
emnr.c
+113 -75
View File
@@ -557,6 +557,7 @@ void calc_emnr(EMNR a)
a->ae.psi = 20.0;
a->ae.t2 = 0.20;
a->ae.nmask = (double *)malloc0(a->ae.msize * sizeof(double));
a->ae.csum = (double *)malloc0((a->ae.msize + 1) * sizeof(double));
//
// post2
a->post2.run = 0;
@@ -580,6 +581,7 @@ void decalc_emnr(EMNR a)
_aligned_free(a->post2.noise_frame);
_aligned_free(a->post2.w);
// ae
_aligned_free(a->ae.csum);
_aligned_free(a->ae.nmask);
// npl
_aligned_free(a->npl.D);
@@ -868,26 +870,26 @@ void aepf(EMNR a)
else
N = 1 + 2 * (int)(0.5 + a->ae.psi * (1.0 - zetaT / a->ae.zetaThresh));
n = N / 2;
for (k = 0; k < n; k++)
/* Each of the three spans below averages mask[] over a window that is
symmetric about k and clipped at the array ends. Taking them straight
from a prefix sum makes each output one subtraction rather than a walk of
up to N = 2*psi + 1 taps, so the pass is O(msize) instead of O(msize*N). */
{
a->ae.nmask[k] = 0.0;
for (m = 0; m <= 2 * k; m++)
a->ae.nmask[k] += a->mask[m];
a->ae.nmask[k] /= (double)(2 * k + 1);
}
for (k = n; k < (a->ae.msize - n); k++)
{
a->ae.nmask[k] = 0.0;
for (m = k - n; m <= (k + n); m++)
a->ae.nmask[k] += a->mask[m];
a->ae.nmask[k] /= (double)N;
}
for (k = a->ae.msize - n; k < a->ae.msize; k++)
{
a->ae.nmask[k] = 0.0;
for (m = (a->ae.msize - 1); m >= (-a->ae.msize + 2 * k + 1); m--)
a->ae.nmask[k] += a->mask[m];
a->ae.nmask[k] /= (double)(2 * (a->ae.msize - k) - 1);
const int msize = a->ae.msize;
const double* WDSP_RESTRICT mask = a->mask;
double* WDSP_RESTRICT nmask = a->ae.nmask;
double* WDSP_RESTRICT csum = a->ae.csum;
csum[0] = 0.0;
for (k = 0; k < msize; k++)
csum[k + 1] = csum[k] + mask[k];
for (k = 0; k < n; k++) // window [0, 2k]
nmask[k] = (csum[2 * k + 1] - csum[0]) / (double)(2 * k + 1);
for (k = n; k < (msize - n); k++) // window [k-n, k+n]
nmask[k] = (csum[k + n + 1] - csum[k - n]) / (double)N;
for (k = msize - n; k < msize; k++) // window [2k+1-msize, msize-1]
nmask[k] = (csum[msize] - csum[2 * k + 1 - msize]) / (double)(2 * (msize - k) - 1);
}
memcpy (a->mask, a->ae.nmask, a->ae.msize * sizeof (double));
if (a->g.gain_method == 3 && zetaT < a->ae.t2)
@@ -998,52 +1000,54 @@ void SetRXAEMNRpost2Rate(int channel, double tc)
* End Post-Processing Functions *
********************************************************************************************************/
double getKey(double* type, double gamma, double xi)
/* Locate v on the table's 0.25 dB grid: n1/n2 bracket it, d is the fraction.
v is compared against the table's dB limits first, so on the interpolating
path the log argument is positive and normal and wdsp_log10 applies. */
static inline void keyIndex (double v, int* n1, int* n2, double* d)
{
int ngamma1, ngamma2, nxi1, nxi2;
double tg, tx, dg, dx;
const double dmin = 0.001;
const double dmax = 1000.0;
if (gamma <= dmin)
if (v <= dmin)
{
ngamma1 = ngamma2 = 0;
tg = 0.0;
*n1 = *n2 = 0;
*d = 0.0;
}
else if (gamma >= dmax)
else if (v >= dmax)
{
ngamma1 = ngamma2 = 240;
tg = 60.0;
*n1 = *n2 = 240;
*d = 0.0;
}
else
{
tg = 10.0 * log10(gamma / dmin);
ngamma1 = (int)(4.0 * tg);
ngamma2 = ngamma1 + 1;
double f = 40.0 * wdsp_log10 (v / dmin); // 4 * (10 * log10)
int i = (int)f;
/* clamp so n2 cannot address the next row of the 241x241 table */
if (i > 239) i = 239;
*n1 = i;
*n2 = i + 1;
*d = f - (double)i;
}
if (xi <= dmin)
}
static inline double keyLerp (const double* type, int ngamma1, int ngamma2, double dg,
int nxi1, int nxi2, double dx)
{
nxi1 = nxi2 = 0;
tx = 0.0;
}
else if (xi >= dmax)
{
nxi1 = nxi2 = 240;
tx = 60.0;
}
else
{
tx = 10.0 * log10(xi / dmin);
nxi1 = (int)(4.0 * tx);
nxi2 = nxi1 + 1;
}
dg = (tg - 0.25 * ngamma1) / 0.25;
dx = (tx - 0.25 * nxi1) / 0.25;
return (1.0 - dg) * (1.0 - dx) * type[241 * nxi1 + ngamma1]
+ (1.0 - dg) * dx * type[241 * nxi2 + ngamma1]
+ dg * (1.0 - dx) * type[241 * nxi1 + ngamma2]
+ dg * dx * type[241 * nxi2 + ngamma2];
}
double getKey(double* type, double gamma, double xi)
{
int ngamma1, ngamma2, nxi1, nxi2;
double dg, dx;
keyIndex (gamma, &ngamma1, &ngamma2, &dg);
keyIndex (xi, &nxi1, &nxi2, &dx);
return keyLerp (type, ngamma1, ngamma2, dg, nxi1, nxi2, dx);
}
int getZeta( EMNR a, double gamma, double eps, double* zeta)
{
int index, i_gamma, i_xi;
@@ -1133,13 +1137,20 @@ void calc_gain (EMNR a)
case 2:
{
double gamma, eps_hat, eps_p;
int ngamma1, ngamma2, nxi1, nxi2, npi1, npi2;
double dg, dx, dp;
for (k = 0; k < a->g.msize; k++)
{
gamma = min(a->g.lambda_y[k] / a->g.lambda_d[k], a->g.gamma_max);
eps_hat = a->g.alpha * a->g.prev_mask[k] * a->g.prev_mask[k] * a->g.prev_gamma[k]
+ (1.0 - a->g.alpha) * max(gamma - 1.0, a->g.eps_floor);
eps_p = eps_hat / (1.0 - a->g.q);
a->g.mask[k] = getKey(a->g.GG, gamma, eps_hat) * getKey(a->g.GGS, gamma, eps_p);
/* both lookups share gamma, so locate it once */
keyIndex (gamma, &ngamma1, &ngamma2, &dg);
keyIndex (eps_hat, &nxi1, &nxi2, &dx);
keyIndex (eps_p, &npi1, &npi2, &dp);
a->g.mask[k] = keyLerp (a->g.GG, ngamma1, ngamma2, dg, nxi1, nxi2, dx)
* keyLerp (a->g.GGS, ngamma1, ngamma2, dg, npi1, npi2, dp);
a->g.prev_gamma[k] = gamma;
a->g.prev_mask[k] = a->g.mask[k];
}
@@ -1203,18 +1214,34 @@ void xemnr (EMNR a, int pos)
{
int i, j, k, sbuff, sbegin;
double g1;
for (i = 0; i < 2 * a->bsize; i += 2)
/* The ring indices below advance by one per iteration and wrap at most
once per loop, so a '% size' each step is an integer division for
nothing (iasize = 3584 here, not a power of two). Walk contiguous runs
and wrap between them instead. */
const int iasize = a->iasize;
const int oasize = a->oasize;
const int fsize = a->fsize;
const int incr = a->incr;
const int bsize = a->bsize;
const int ovrlp = a->ovrlp;
for (i = 0, j = a->iainidx; i < 2 * bsize; i += 2)
{
a->inaccum[a->iainidx] = a->in[i];
a->iainidx = (a->iainidx + 1) % a->iasize;
a->inaccum[j] = a->in[i];
if (++j == iasize) j = 0;
}
a->nsamps += a->bsize;
while (a->nsamps >= a->fsize)
a->iainidx = j;
a->nsamps += bsize;
while (a->nsamps >= fsize)
{
for (i = 0, j = a->iaoutidx; i < a->fsize; i++, j = (j + 1) % a->iasize)
a->forfftin[i] = a->window[i] * a->inaccum[j];
a->iaoutidx = (a->iaoutidx + a->incr) % a->iasize;
a->nsamps -= a->incr;
int n1 = iasize - a->iaoutidx;
if (n1 > fsize) n1 = fsize;
for (i = 0; i < n1; i++)
a->forfftin[i] = a->window[i] * a->inaccum[a->iaoutidx + i];
for (; i < fsize; i++)
a->forfftin[i] = a->window[i] * a->inaccum[i - n1];
if ((a->iaoutidx += incr) >= iasize) a->iaoutidx -= iasize;
a->nsamps -= incr;
fftw_execute (a->Rfor);
calc_gain(a);
for (i = 0; i < a->msize; i++)
@@ -1225,30 +1252,41 @@ void xemnr (EMNR a, int pos)
}
post2(a);
fftw_execute (a->Rrev);
for (i = 0; i < a->fsize; i++)
for (i = 0; i < fsize; i++)
a->save[a->saveidx][i] = a->window[i] * a->revfftout[i];
for (i = a->ovrlp; i > 0; i--)
for (i = ovrlp; i > 0; i--)
{
sbuff = (a->saveidx + i) % a->ovrlp;
sbegin = a->incr * (a->ovrlp - i);
for (j = sbegin, k = a->oainidx; j < a->incr + sbegin; j++, k = (k + 1) % a->oasize)
const double* WDSP_RESTRICT sv;
double* WDSP_RESTRICT oa = a->outaccum;
int m1;
sbuff = (a->saveidx + i) % ovrlp;
sbegin = incr * (ovrlp - i);
sv = a->save[sbuff] + sbegin;
m1 = oasize - a->oainidx;
if (m1 > incr) m1 = incr;
k = a->oainidx;
if (i == ovrlp)
{
if ( i == a->ovrlp)
a->outaccum[k] = a->save[sbuff][j];
for (j = 0; j < m1; j++) oa[k + j] = sv[j];
for (; j < incr; j++) oa[j - m1] = sv[j];
}
else
a->outaccum[k] += a->save[sbuff][j];
}
}
a->saveidx = (a->saveidx + 1) % a->ovrlp;
a->oainidx = (a->oainidx + a->incr) % a->oasize;
}
for (i = 0; i < a->bsize; i++)
{
a->out[2 * i + 0] = a->outaccum[a->oaoutidx];
a->out[2 * i + 1] = 0.0;
a->oaoutidx = (a->oaoutidx + 1) % a->oasize;
for (j = 0; j < m1; j++) oa[k + j] += sv[j];
for (; j < incr; j++) oa[j - m1] += sv[j];
}
}
if (++a->saveidx == ovrlp) a->saveidx = 0;
if ((a->oainidx += incr) >= oasize) a->oainidx -= oasize;
}
for (i = 0, k = a->oaoutidx; i < bsize; i++)
{
a->out[2 * i + 0] = a->outaccum[k];
a->out[2 * i + 1] = 0.0;
if (++k == oasize) k = 0;
}
a->oaoutidx = k;
}
else if (a->out != a->in)
memcpy (a->out, a->in, a->bsize * sizeof (complex));
}
emnr.h
+1
View File
@@ -185,6 +185,7 @@ typedef struct _emnr
double zetaThresh;
double psi;
double* nmask;
double* csum; // prefix sums of mask[], msize + 1 entries
double t2;
} ae;
struct _post2
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)
{
int i, n;
for (i = 0; i < a->size; i++)
const int size = a->size;
const int nstages = a->nstages;
const double b0 = a->b0, b1 = a->b1, a1 = a->a1;
/* in and out are the same buffer in TXA, so neither may be restrict */
const double* in = a->in;
double* out = a->out;
/* x0[]/y0[] never carried state between samples: x0[n] was only ever the
previous stage's output and y0[n] this stage's. Keep that single value
in a register and cascade it, leaving x1/y1 as the actual filter state.
Stores through out[] could alias the struct's doubles, so hoist the
coefficients too. */
double* WDSP_RESTRICT x1 = a->x1;
double* WDSP_RESTRICT y1 = a->y1;
for (i = 0; i < size; i++)
{
a->x0[0] = a->in[2 * i + 0];
for (n = 0; n < a->nstages; n++)
double v = in[2 * i + 0];
for (n = 0; n < nstages; n++)
{
if (n > 0) a->x0[n] = a->y0[n - 1];
a->y0[n] = a->b0 * a->x0[n]
+ a->b1 * a->x1[n]
- a->a1 * a->y1[n];
a->y1[n] = a->y0[n];
a->x1[n] = a->x0[n];
double y = b0 * v + b1 * x1[n] - a1 * y1[n];
x1[n] = v;
y1[n] = y;
v = y;
}
a->out[2 * i + 0] = a->y0[a->nstages - 1];
out[2 * i + 0] = v;
}
}
else if (a->out != a->in)
linux_port.c
+11 -3
View File
@@ -39,9 +39,17 @@ john.d.melton@googlemail.com
#if defined(linux) || defined(__APPLE__)
void QueueUserWorkItem(void *function,void *context,int flags) {
pthread_t t;
pthread_create(&t, NULL, function, context);
pthread_join(t, NULL);
//
// The Windows call queues the work item on a thread pool and returns at
// once, so callers get their items run in parallel. This shim spawned a
// thread and immediately joined it, which is a plain synchronous call that
// happens to cost a thread creation (~15 us) and delivers no parallelism.
// Call the function directly: same ordering, no thread.
//
// The cast matches the one pthread_create() performed here before.
//
(void)flags;
((void *(*)(void *))function)(context);
}
static inline void init_crit_section(pthread_mutex_t *mutex) {
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)
{
int x, y, z;
@@ -71,7 +102,8 @@ void calc_resample (RESAMPLE a)
for (k = 0; k < a->ncoef; k += a->L)
a->h[i++] = impulse[j + k];
a->ringsize = a->cpp;
a->ring = (double *)malloc0(a->ringsize * sizeof(complex));
a->ringI = (double *)malloc0(a->ringsize * sizeof(double));
a->ringQ = (double *)malloc0(a->ringsize * sizeof(double));
a->idx_in = a->ringsize - 1;
a->phnum = 0;
_aligned_free(impulse);
@@ -79,7 +111,8 @@ void calc_resample (RESAMPLE a)
void decalc_resample (RESAMPLE a)
{
_aligned_free(a->ring);
_aligned_free(a->ringQ);
_aligned_free(a->ringI);
_aligned_free(a->h);
}
@@ -112,7 +145,8 @@ void destroy_resample (RESAMPLE a)
PORT
void flush_resample (RESAMPLE a)
{
memset (a->ring, 0, a->ringsize * sizeof (complex));
memset (a->ringI, 0, a->ringsize * sizeof (double));
memset (a->ringQ, 0, a->ringsize * sizeof (double));
a->idx_in = a->ringsize - 1;
a->phnum = 0;
}
@@ -123,40 +157,49 @@ int xresample (RESAMPLE a)
int outsamps = 0;
if (a->run)
{
int i, j, n;
int idx_out;
int i, n1;
double I, Q;
int cpp = a->cpp;
const int cpp = a->cpp;
const int ringsize = a->ringsize;
const int L = a->L;
const int M = a->M;
const int size = a->size;
const double* WDSP_RESTRICT h = a->h;
const double* WDSP_RESTRICT in = a->in;
double* WDSP_RESTRICT ringI = a->ringI;
double* WDSP_RESTRICT ringQ = a->ringQ;
double* WDSP_RESTRICT out = a->out;
int idx_in = a->idx_in;
int ringsize = a->ringsize;
double* h = a->h;
double* ring = a->ring;
int phnum = a->phnum;
for (i = 0; i < a->size; i++)
for (i = 0; i < size; i++)
{
ring[2 * idx_in + 0] = a->in[2 * i + 0];
ring[2 * idx_in + 1] = a->in[2 * i + 1];
while (a->phnum < a->L)
ringI[idx_in] = in[2 * i + 0];
ringQ[idx_in] = in[2 * i + 1];
while (phnum < L)
{
const double* WDSP_RESTRICT hp = h + cpp * phnum;
/* The tap loop walks the ring forward from idx_in and wraps at
most once. Split it at the wrap point so both halves are
unit-stride: the wrap test that used to sit inside the loop
made the address non-affine and blocked vectorization. */
if ((n1 = ringsize - idx_in) > cpp) n1 = cpp;
I = 0.0;
Q = 0.0;
n = cpp * a->phnum;
for (j = 0; j < cpp; j++)
{
if ((idx_out = idx_in + j) >= ringsize) idx_out -= ringsize;
I += h[n + j] * ring[2 * idx_out + 0];
Q += h[n + j] * ring[2 * idx_out + 1];
}
a->out[2 * outsamps + 0] = I;
a->out[2 * outsamps + 1] = Q;
resample_dot (hp, ringI + idx_in, ringQ + idx_in, n1, &I, &Q);
if (n1 < cpp)
resample_dot (hp + n1, ringI, ringQ, cpp - n1, &I, &Q);
out[2 * outsamps + 0] = I;
out[2 * outsamps + 1] = Q;
outsamps++;
a->phnum += a->M;
phnum += M;
}
a->phnum -= a->L;
if (--idx_in < 0) idx_in = a->ringsize - 1;
phnum -= L;
if (--idx_in < 0) idx_in = ringsize - 1;
}
a->idx_in = idx_in;
a->phnum = phnum;
}
else if (a->in != a->out)
memcpy (a->out, a->in, a->size * sizeof (complex));
@@ -240,6 +283,24 @@ void destroy_resampleV (void* ptr)
* *
************************************************************************************************/
/* Real-valued counterpart of resample_dot(). */
static inline void resampleF_dot (const double* WDSP_RESTRICT hp,
const double* WDSP_RESTRICT x, int n, double* pI)
{
double i0 = 0.0, i1 = 0.0, i2 = 0.0, i3 = 0.0;
int j = 0;
for (; j <= n - 4; j += 4)
{
i0 += hp[j + 0] * x[j + 0];
i1 += hp[j + 1] * x[j + 1];
i2 += hp[j + 2] * x[j + 2];
i3 += hp[j + 3] * x[j + 3];
}
for (; j < n; j++)
i0 += hp[j] * x[j];
*pI += (i0 + i1) + (i2 + i3);
}
RESAMPLEF create_resampleF ( int run, int size, float* in, float* out, int in_rate, int out_rate)
{
RESAMPLEF a = (RESAMPLEF) malloc0 (sizeof (resampleF));
@@ -305,31 +366,46 @@ int xresampleF (RESAMPLEF a)
int outsamps = 0;
if (a->run)
{
int i, j, n;
int idx_out;
int i;
double I;
for (i = 0; i < a->size; i++)
{
a->ring[a->idx_in] = (double)a->in[i];
const int cpp = a->cpp;
const int ringsize = a->ringsize;
const int L = a->L;
const int M = a->M;
const int size = a->size;
const double* WDSP_RESTRICT h = a->h;
const float* WDSP_RESTRICT in = a->in;
double* WDSP_RESTRICT ring = a->ring;
float* WDSP_RESTRICT out = a->out;
int idx_in = a->idx_in;
int phnum = a->phnum;
int n1;
while (a->phnum < a->L)
for (i = 0; i < size; i++)
{
ring[idx_in] = (double)in[i];
while (phnum < L)
{
const double* WDSP_RESTRICT hp = h + cpp * phnum;
/* see resample_dot(): split at the ring wrap so both halves are
unit-stride, and carry independent accumulators */
if ((n1 = ringsize - idx_in) > cpp) n1 = cpp;
I = 0.0;
n = a->cpp * a->phnum;
for (j = 0; j < a->cpp; j++)
{
if ((idx_out = a->idx_in + j) >= a->ringsize) idx_out -= a->ringsize;
I += a->h[n + j] * a->ring[idx_out];
}
a->out[outsamps] = (float)I;
resampleF_dot (hp, ring + idx_in, n1, &I);
if (n1 < cpp)
resampleF_dot (hp + n1, ring, cpp - n1, &I);
out[outsamps] = (float)I;
outsamps++;
a->phnum += a->M;
phnum += M;
}
a->phnum -= a->L;
if (--a->idx_in < 0) a->idx_in = a->ringsize - 1;
phnum -= L;
if (--idx_in < 0) idx_in = ringsize - 1;
}
a->idx_in = idx_in;
a->phnum = phnum;
}
else if (a->in != a->out)
memcpy (a->out, a->in, a->size * sizeof (float));
resample.h
+3 -1
View File
@@ -52,7 +52,9 @@ typedef struct _resample
int M; // decimation factor
double* h; // coefficients
int ringsize; // number of complex pairs the ring buffer holds
double* ring; // ring buffer
double* ringI; // ring buffer, in-phase
double* ringQ; // ring buffer, quadrature (split from I so the tap loop
// reads unit-stride and vectorizes)
int cpp; // coefficients of the phase
int phnum; // phase number
} resample, *RESAMPLE;
varsamp.c
+88 -23
View File
@@ -59,10 +59,27 @@ void calc_varsamp (VARSAMP a)
fc_norm_low = a->fc_low / norm_rate;
a->rsize = (int)(140.0 * norm_rate / min_rate);
a->ncoef = a->rsize + 1;
a->ncoef += (a->R - 1) * (a->ncoef - 1);
a->h = fir_bandpass(a->ncoef, fc_norm_low, fc_norm_high, (double)a->R, 1, 0, (double)a->R * a->gain);
// print_impulse ("imp.txt", a->ncoef, a->h, 0, 0);
a->ring = (double *)malloc0(a->rsize * sizeof(complex));
a->ncoef += (a->R - 1) * (a->ncoef - 1); // = R * rsize + 1
{
/* Store the coefficients transposed into phases. hshift() walks
h[hidx + m*R] for m = 0..rsize-1, which strides by R doubles -- 8 KB
at R = 1024 -- over a 1.1 MB table, so every tap is its own cache
line. Transposing makes each phase contiguous; hshift() interpolates
between phases hidx and hidx+1, hence R+1 of them. */
int p, m;
const int R = a->R, rsize = a->rsize;
double* h = fir_bandpass(a->ncoef, fc_norm_low, fc_norm_high, (double)R, 1, 0, (double)R * a->gain);
// every element is written below, so skip malloc0()'s memset of ~1 MB
a->hp = (double *)_aligned_malloc ((size_t)(R + 1) * rsize * sizeof (double), 16);
// walk h forward (p is its fast axis) so the prefetcher sees a linear
// stream; h is cold here, straight from fir_bandpass()
for (m = 0; m < rsize; m++)
for (p = 0; p <= R; p++)
a->hp[(size_t)p * rsize + m] = h[p + (size_t)m * R];
_aligned_free (h);
}
a->ringI = (double *)malloc0(a->rsize * sizeof(double));
a->ringQ = (double *)malloc0(a->rsize * sizeof(double));
a->idx_in = a->rsize - 1;
a->h_offset = 0.0;
a->hs = (double *)malloc0 (a->rsize * sizeof (double));
@@ -72,8 +89,9 @@ void calc_varsamp (VARSAMP a)
void decalc_varsamp (VARSAMP a)
{
_aligned_free (a->hs);
_aligned_free (a->ring);
_aligned_free (a->h);
_aligned_free (a->ringQ);
_aligned_free (a->ringI);
_aligned_free (a->hp);
}
VARSAMP create_varsamp ( int run, int size, double* in, double* out,
@@ -105,22 +123,60 @@ void destroy_varsamp (VARSAMP a)
void flush_varsamp (VARSAMP a)
{
memset (a->ring, 0, a->rsize * sizeof (complex));
memset (a->ringI, 0, a->rsize * sizeof (double));
memset (a->ringQ, 0, a->rsize * sizeof (double));
a->idx_in = a->rsize - 1;
a->h_offset = 0.0;
a->isamps = 0.0;
}
/* Accumulate n taps of a unit-stride complex dot product into *pI / *pQ.
Four independent accumulator pairs keep the FMAs off a single dependency
chain and let the vectorizer in: an 'I += h[j]*x[j]' reduction cannot be
reassociated without -ffast-math, which this library must not enable (it
relies on IEEE semantics for 0/0 = NaN and x/0 = Inf). */
static inline void varsamp_dot (const double* WDSP_RESTRICT hp,
const double* WDSP_RESTRICT xI, const double* WDSP_RESTRICT xQ,
int n, double* pI, double* pQ)
{
double i0 = 0.0, i1 = 0.0, i2 = 0.0, i3 = 0.0;
double q0 = 0.0, q1 = 0.0, q2 = 0.0, q3 = 0.0;
int j = 0;
for (; j <= n - 4; j += 4)
{
i0 += hp[j + 0] * xI[j + 0]; q0 += hp[j + 0] * xQ[j + 0];
i1 += hp[j + 1] * xI[j + 1]; q1 += hp[j + 1] * xQ[j + 1];
i2 += hp[j + 2] * xI[j + 2]; q2 += hp[j + 2] * xQ[j + 2];
i3 += hp[j + 3] * xI[j + 3]; q3 += hp[j + 3] * xQ[j + 3];
}
for (; j < n; j++)
{
i0 += hp[j] * xI[j];
q0 += hp[j] * xQ[j];
}
*pI += (i0 + i1) + (i2 + i3);
*pQ += (q0 + q1) + (q2 + q3);
}
void hshift (VARSAMP a)
{
int i, j, k;
int m;
int hidx;
double frac, pos;
const int rsize = a->rsize;
const double* WDSP_RESTRICT h0;
const double* WDSP_RESTRICT h1;
double* WDSP_RESTRICT hs = a->hs;
/* h_offset is normalized to [0,1) by the caller, so hidx is in [0, R-1]
and phase hidx+1 <= R exists. */
pos = (double)a->R * a->h_offset;
hidx = (int)(pos);
frac = pos - (double)hidx;
for (i = a->rsize - 1, j = hidx, k = hidx + 1; i >= 0; i--, j += a->R, k += a->R)
a->hs[i] = a->h[j] + frac * (a->h[k] - a->h[j]);
h0 = a->hp + (size_t)hidx * rsize;
h1 = h0 + rsize;
for (m = 0; m < rsize; m++)
hs[rsize - 1 - m] = h0[m] + frac * (h1[m] - h0[m]);
}
int xvarsamp (VARSAMP a, double var)
@@ -140,13 +196,21 @@ int xvarsamp (VARSAMP a, double var)
else a->dicvar = 0.0;
if (a->run)
{
int i, j;
int idx_out;
int i, n1;
double I, Q;
const int rsize = a->rsize;
/* a->hs is rewritten by hshift() on every output sample, so it must not
be hoisted behind a restrict pointer here; varsamp_dot() re-reads it. */
const double* in = a->in;
double* out = a->out;
double* WDSP_RESTRICT ringI = a->ringI;
double* WDSP_RESTRICT ringQ = a->ringQ;
int idx_in = a->idx_in;
for (i = 0; i < a->size; i++)
{
a->ring[2 * a->idx_in + 0] = a->in[2 * i + 0];
a->ring[2 * a->idx_in + 1] = a->in[2 * i + 1];
ringI[idx_in] = in[2 * i + 0];
ringQ[idx_in] = in[2 * i + 1];
a->inv_cvar += a->dicvar;
picvar = (uint64_t*)(&a->inv_cvar);
N = *picvar & 0xffffffffffff0000;
@@ -160,20 +224,21 @@ int xvarsamp (VARSAMP a, double var)
a->h_offset += a->delta;
while (a->h_offset >= 1.0) a->h_offset -= 1.0;
while (a->h_offset < 0.0) a->h_offset += 1.0;
for (j = 0; j < a->rsize; j++)
{
if ((idx_out = a->idx_in + j) >= a->rsize) idx_out -= a->rsize;
I += a->hs[j] * a->ring[2 * idx_out + 0];
Q += a->hs[j] * a->ring[2 * idx_out + 1];
}
a->out[2 * outsamps + 0] = I;
a->out[2 * outsamps + 1] = Q;
/* the ring wraps at most once over rsize taps; split it so both
halves are unit-stride */
n1 = rsize - idx_in;
varsamp_dot (a->hs, ringI + idx_in, ringQ + idx_in, n1, &I, &Q);
if (n1 < rsize)
varsamp_dot (a->hs + n1, ringI, ringQ, rsize - n1, &I, &Q);
out[2 * outsamps + 0] = I;
out[2 * outsamps + 1] = Q;
outsamps++;
a->isamps += a->inv_cvar;
}
a->isamps -= 1.0;
if (--a->idx_in < 0) a->idx_in = a->rsize - 1;
if (--idx_in < 0) idx_in = rsize - 1;
}
a->idx_in = idx_in;
}
else if (a->in != a->out)
memcpy (a->out, a->in, a->size * sizeof (complex));
varsamp.h
+6 -2
View File
@@ -41,9 +41,13 @@ typedef struct _varsamp
double gain;
int idx_in;
int ncoef;
double* h;
double* hp; // coefficients, polyphase: hp[p * rsize + m] = h[p + m * R],
// p = 0..R. hshift() then reads two adjacent phases
// contiguously instead of striding by R.
int rsize;
double* ring;
double* ringI; // ring buffer, in-phase
double* ringQ; // ring buffer, quadrature (split from I so the tap loop
// reads unit-stride and vectorizes)
double var;
int varmode;
double cvar;
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 GetPSMaxTX (int channel, double* maxtx);
extern void SetPSPtol (int channel, double ptol);
extern void SetPSOutlierSigma (int channel, double sigma);
extern void GetPSDisp (int channel, double* x, double* ym, double* yc, double* ys, double* cm, double* cc, double* cs);
extern void SetPSFeedbackRate (int channel, int rate);
extern void SetPSPinMode (int channel, int pin);
@@ -517,6 +518,30 @@ extern void SetTXAFMNC (int channel, int nc);
extern void SetTXAFMMP (int channel, int mp);
extern void SetTXAFMAFFreqs (int channel, double low, double high);
//
// Interfaces from wfmd.c
//
extern void SetRXAWFMDeviation (int channel, double deviation);
extern void SetRXAWFMNCaud (int channel, int nc);
extern void SetRXAWFMMPaud (int channel, int mp);
extern void SetRXAWFMAFFilter (int channel, double low, double high);
extern void SetRXAWFMDeemphRun (int channel, int run);
extern void SetRXAWFMDeemphTau (int channel, double tau);
extern void SetRXAWFMLimRun (int channel, int run);
extern void SetRXAWFMLimGain (int channel, double gaindB);
//
// Interfaces from wfmmod.c
//
extern void SetTXAWFMDeviation (int channel, double deviation);
extern void SetTXAWFMNC (int channel, int nc);
extern void SetTXAWFMMP (int channel, int mp);
extern void SetTXAWFMAFFreqs (int channel, double low, double high);
extern void SetTXAWFMPreEmphRun (int channel, int run);
extern void SetTXAWFMPreEmphTau (int channel, double tau);
//
// Interfaces from fmsq.c
//
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