From e91ead679b5e98f3791c4d9346ec28c4886de66b Mon Sep 17 00:00:00 2001 From: Brad Barnett <127794626+bdbarnett@users.noreply.github.com> Date: Fri, 30 Oct 2026 03:12:09 -0500 Subject: [PATCH] audiopump: DriftResampler, which holds a Ring filled on another clock at its target fill A 32-tap windowed-sinc interpolator with a 32.32 read position, steered by a PI loop on the Ring's smoothed fill. The kernel is shared C; MicroPython, CircuitPython and CPython bind it and render it byte for byte. It lives in audiopump because CircuitPython owns audiospeed. --- ROADMAP.md | 3 + circuitpython.mk | 1 + docs/how-it-works/additions.md | 3 +- docs/how-it-works/drift-resampler.md | 122 ++ docs/how-it-works/index.md | 1 + docs/reference/audiopump.md | 35 +- micropython.cmake | 2 + micropython.mk | 4 +- setup.py | 1 + src/audiopump/audiopump.c | 2 + src/audiopump/audiopump_drift.c | 295 ++++ src/audiopump/audiopump_drift.h | 10 + src/audiopump/audiopump_ring.c | 40 + src/audiopump/audiopump_ring.h | 15 + .../audiopump/DriftResampler.c | 254 ++++ .../audiopump/DriftResampler.h | 7 + .../shared-bindings/audiopump/__init__.c | 6 + .../shared-module/audiopump/DriftResampler.c | 123 ++ .../shared-module/audiopump/DriftResampler.h | 54 + .../shared-module/audiopump/Ring.c | 33 + .../shared-module/audiopump/Ring.h | 10 + src/cpython/_audiodsp.c | 181 +++ src/cpython/audiopump.py | 140 +- src/shared/audiodsp_drift.c | 270 ++++ src/shared/audiodsp_drift.h | 145 ++ src/shared/audiodsp_drift_table.h | 1305 +++++++++++++++++ tests/parity/drift_probe.py | 93 ++ tests/parity/verify_dsp.py | 3 + tests/test_binding_parity.py | 20 +- tests/test_cpython_drift.py | 269 ++++ tools/design_drift_table.py | 131 ++ 31 files changed, 3564 insertions(+), 14 deletions(-) create mode 100644 docs/how-it-works/drift-resampler.md create mode 100644 src/audiopump/audiopump_drift.c create mode 100644 src/audiopump/audiopump_drift.h create mode 100644 src/circuitpython_spike/shared-bindings/audiopump/DriftResampler.c create mode 100644 src/circuitpython_spike/shared-bindings/audiopump/DriftResampler.h create mode 100644 src/circuitpython_spike/shared-module/audiopump/DriftResampler.c create mode 100644 src/circuitpython_spike/shared-module/audiopump/DriftResampler.h create mode 100644 src/shared/audiodsp_drift.c create mode 100644 src/shared/audiodsp_drift.h create mode 100644 src/shared/audiodsp_drift_table.h create mode 100644 tests/parity/drift_probe.py create mode 100644 tests/test_cpython_drift.py create mode 100644 tools/design_drift_table.py diff --git a/ROADMAP.md b/ROADMAP.md index d57a977f..9c340150 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -15,6 +15,9 @@ audio pump holding up under live use such as a stompbox or a stage rig. - A level ramp on `audiomixer` voices: a level move ramps over a set number of frames instead of stepping, so an effect can fade its wet path in or out without a click when you switch it in mid-note. +- `audiopump.DriftResampler` converting down as well as up, with its filter's + cutoff following the output rate, so a 48 kHz stream can play on a + 44.1 kHz output. ## Housekeeping diff --git a/circuitpython.mk b/circuitpython.mk index 59e3f67d..75334883 100644 --- a/circuitpython.mk +++ b/circuitpython.mk @@ -60,6 +60,7 @@ AUDIODSP_CP_SOURCES += $(addprefix $(AUDIODSP_CP_SRC)/shared/, \ audiodsp_vorbis.c \ audiodsp_port.c \ audiodsp_graph.c \ + audiodsp_drift.c \ ) # The graph's OVERDRIVE node needs double-precision libm the RP2040 build diff --git a/docs/how-it-works/additions.md b/docs/how-it-works/additions.md index 5865b38d..d6c82b79 100644 --- a/docs/how-it-works/additions.md +++ b/docs/how-it-works/additions.md @@ -45,7 +45,8 @@ levels for a spectrum meter) and `audiofeatures` (speech features for a wake word model) are the two CircuitPython doesn't get: their sources are a board's pump and a USB sound card, and CircuitPython has neither. `audiovorbis` plays Ogg Vorbis files, which CircuitPython has no module for. And `audiopump.Ring` lets CircuitPython play audio your own code -computes. +computes, with `audiopump.DriftResampler` to play one filled on another +device's clock ([the drift resampler](drift-resampler.md)). ## Why `audioroute.Port` exists diff --git a/docs/how-it-works/drift-resampler.md b/docs/how-it-works/drift-resampler.md new file mode 100644 index 00000000..5609d867 --- /dev/null +++ b/docs/how-it-works/drift-resampler.md @@ -0,0 +1,122 @@ +# The drift resampler + +`audiopump.DriftResampler` plays a `Ring` that something else fills on its +own clock: a USB microphone, a network stream, a second codec. You hand it +the Ring and play it instead of the Ring, and it keeps the Ring at a target +fill for as long as the two clocks run, without a gap or a skip: + +```python +ring = audiopump.Ring(sample_rate=48000, channel_count=2, frames=256, capacity=8) +drift = audiopump.DriftResampler(ring) # holds the Ring at 1024 frames +mixer.voice[0].play(drift) +``` + +This page is why it is built the way it is, and what it measured. + +## Why a fixed ratio isn't enough + +Two devices that both say 48 kHz differ by tens of parts per million, and +up to a few hundred for cheap crystals. At 100 ppm a Ring between them +gains or loses 4.8 frames a second, so a 1024-frame cushion is gone in +under four minutes. The fix is to read the Ring very slightly faster or +slower than nominal, by exactly the clocks' difference, which nobody knows +in advance and which wanders with temperature. + +`audiospeed.SpeedChanger` holds its rate in 16.16 fixed point, a step of +about 15 ppm: the size of the error being corrected. Here the read position +and the ratio are 32.32 fixed point, a step of about 0.0002 ppm, so the +ratio can sit where the clocks are rather than hunting between two steps +either side of it. + +## Where it lives + +The node is in `audiopump`, beside `Ring`, rather than beside +`audiospeed.Resampler`. CircuitPython has its own `audiospeed`, and a module +of ours never shadows one of CircuitPython's +([the correctness standard](correctness-standard.md)). `audiopump` is ours +on all three runtimes, and the node only reads a Ring, so it belongs with +the Ring. + +## Interpolation + +The read position falls between input samples, so the node interpolates +with a 32-tap Kaiser-windowed sinc (beta 9, cutoff 0.45 of the input rate). +The filter is tabulated at 128 fractional positions, plus one more so a +neighbour always exists, and the node blends linearly between the two rows +either side of the read position. Each row is scaled to sum to exactly 1, so +a constant input comes out unchanged at every position. The table is +[`src/shared/audiodsp_drift_table.h`](../../src/shared/audiodsp_drift_table.h), +written by [`tools/design_drift_table.py`](../../tools/design_drift_table.py); +it is a checked-in constant rather than something each runtime computes, +because CircuitPython's floats are single precision and its libm is not +glibc's, and every target has to render the node identically. + +The cutoff is set for the input rate, so the node converts up or across, +never down: an output rate from the Ring's own up to twice it. Converting +down would fold everything between the two Nyquist frequencies back into +the band. + +**Cost.** Per output frame: 32 coefficient blends, then 32 multiply-adds a +channel, all single-precision float. That is about 96 floating-point +multiply-adds a stereo frame, 4.6 million a second at 48 kHz, or roughly 3% +of one ESP32-P4 core at 360 MHz. That figure is an estimate from the +operation count; it has not been timed on a board. The work buffer is the block plus 36 +frames, as floats. + +## Control + +Once a block, the node measures the fill: the frames the Ring holds plus +those already taken into the filter beyond its look-ahead, less the +fractional read position. The ring's level jumps by a whole packet each +time the producer writes, so the measurement goes through two one-pole +smoothers before a PI loop. The loop's output is the ratio's offset from +nominal, clamped to `max_ppm` (1000 by default), and the integrator is +clamped to the same bound so it cannot wind up. + +The loop is a second-order system with damping 0.8 and natural frequency +`bandwidth` (0.01 Hz by default, 0.001 to 0.1 allowed). It starts at eight +times that and narrows with a 15-second time constant, so it locks in +seconds and then holds still. The smoothers sit at four times the loop's +frequency. Above about 0.1 Hz the loop starts following the producer's +packets instead of its clock, and the ratio swings by hundreds of ppm from +one block to the next, so the bound stops there. + +Playback starts when the fill first reaches the target, and anything beyond +it is dropped on that block, so the loop does not start by draining an +excess. If the Ring runs short of a block, the node plays a block of +silence, counts an underrun on itself and on the Ring, and waits for the +target again. It keeps the rate it had learned, since the clocks have not +changed. + +## What it measured + +On the desktop, through the CPython build, with a producer writing a Ring +of 8 blocks of 256 frames (target 1024) on a clock offset from the output's. +Two packet sizes: 48 frames (1 ms, a USB microphone) and 480 (10 ms, a +network or Bluetooth burst). 90 to 120 seconds each, default settings. + +| Producer's clock | Packet | Underruns, overruns | Within 8 frames of target after | Ratio after lock | +| --- | --- | --- | --- | --- | +| +100 ppm | 48 | 0, 0 | 2.3 s | 100.00 ± 0.36 ppm | +| −100 ppm | 48 | 0, 0 | 6.2 s | −99.94 ± 0.36 ppm | +| +300 ppm | 48 | 0, 0 | 5.6 s | — | +| +100 ppm | 480 | 0, 0 | 32 s | 101.7 ± 6.4 ppm | +| −100 ppm | 480 | 0, 0 | 28 s | −98.1 ± 5.2 ppm | +| ±100 ppm over 2 minutes | 48 and 480 | 0, 0 | tracks within 54 frames | follows it | + +The 480-frame cases take longer because the first measurement lands +somewhere in the packet's sawtooth, up to about 100 frames from where the +smoothed fill settles. The ±5 ppm left on the ratio there is the packets +showing through the smoothers, a fraction of a cent of pitch moving below +1 Hz. + +**THD+N**, a 0.9 full-scale sine in, everything from 20 Hz to 20 kHz but the +fundamental out: −93.7 to −95.2 dB at 1 kHz and 10 kHz across every case +above. A sine quantized once to 16 bits measures −97.5 dB at 1 kHz on the +same analysis, and the node's input and its output are each 16-bit, so two +roundings alone come to about −94.5 dB. The interpolation adds little to +what 16-bit audio already has. + +`tests/test_cpython_drift.py` holds these as traits, with bars a little +looser than the figures, and `tests/parity/drift_probe.py` checks that +MicroPython, CircuitPython and CPython render the node byte for byte. diff --git a/docs/how-it-works/index.md b/docs/how-it-works/index.md index cfd6d711..88ea7f8c 100644 --- a/docs/how-it-works/index.md +++ b/docs/how-it-works/index.md @@ -137,5 +137,6 @@ push-PCM transports. - [Differences from upstream, with measurements](upstream-diff.md): the reasons behind each one. - [Upstream reports](../upstream-reports/README.md): the fixes sent to CircuitPython. - [The graph runner](audiograph.md): how `audiograph.Graph` renders a patch. +- [The drift resampler](drift-resampler.md): how `audiopump.DriftResampler` holds a Ring filled on another clock. - [Writing a pump port](pump-ports.md): how a platform plugs in a clock, a thread and a sink. - [Wheels and platforms](wheels-and-platforms.md): what the release builds, and what it skips. diff --git a/docs/reference/audiopump.md b/docs/reference/audiopump.md index c6fcecf5..df832efb 100644 --- a/docs/reference/audiopump.md +++ b/docs/reference/audiopump.md @@ -30,6 +30,33 @@ the loop hands the thread back on a full ring and the caller's next `park_spin` — a driver that is not finished, which now says so rather than quietly losing blocks. +## `DriftResampler` + +```python +drift = audiopump.DriftResampler(ring, *, sample_rate=None, target=None, + max_ppm=1000, bandwidth=0.01) +``` + +Plays `ring`, which something else fills on its own clock, and steers its +ratio so the Ring holds `target` frames. Play the DriftResampler, not the +Ring: it is the Ring's only reader. On every runtime, CircuitPython included. + +| Argument | Meaning | +| --- | --- | +| `sample_rate` | The output's rate: the Ring's own (the default) up to twice it. | +| `target` | The fill to hold, in frames. Half the Ring by default; at least a block, and at most the Ring less 16 frames, the filter's look-ahead. Settable while playing; the fill moves there at the loop's pace. | +| `max_ppm` | How far the ratio may move from nominal, 0 to 10000. | +| `bandwidth` | How fast the loop follows the clocks once locked, 0.001 to 0.1 Hz. It starts eight times faster and narrows over about a minute. | + +The Ring's block must be 16 to 1024 frames. `ppm` (the ratio's offset from +nominal), `fill` (the smoothed fill it is holding) and `rate` (Ring frames per +output frame) say what the loop is doing; `underruns` counts blocks of +silence played because the Ring ran dry, after which it waits for `target` +again. `reset()` starts over and forgets the rate. Why it is built this way, +and what it measured: [the drift resampler](../how-it-works/drift-resampler.md). + +## What feeds the pump + Three kinds of source feed it and they are all the same kind of thing to it: a graph; `audiopump.Ring`, an audiosample node Python writes PCM into, so a pushed stream can sit behind a Mixer with effects on it like anything else; @@ -62,10 +89,10 @@ while playing: # from your app's timer, say every 20 ms their own tick. What CPython does not have is `pull()` under a heap lock, because it has no heap lock, and audiocore's pull-depth cap: a graph that leads back into itself is caught at Python's recursion limit instead, and -reported as the same fault. **CircuitPython carries only `Ring`**: its playback -layer already pulls natively from its own audio thread, and a second pull loop -would be a second owner of the same graph, so there a board's audio output -pulls the `Ring` itself X. It still takes the lock and the +reported as the same fault. **CircuitPython carries only `Ring` and +`DriftResampler`**: its playback layer already pulls natively from its own +audio thread, and a second pull loop would be a second owner of the same +graph, so there a board's audio output pulls the `Ring` itself. It still takes the lock and the hook table, on default hooks, which cost a load and a branch. The hardware — the channel a board writes into, and a live microphone as a diff --git a/micropython.cmake b/micropython.cmake index 2b3290cd..3600171c 100644 --- a/micropython.cmake +++ b/micropython.cmake @@ -30,6 +30,8 @@ target_sources(usermod_mpaudio INTERFACE ${MPAUDIO_SRC_DIR}/audiopump/audiopump_ring.c ${MPAUDIO_SRC_DIR}/audiopump/audiopump_events.c ${MPAUDIO_SRC_DIR}/audiopump/audiopump_tap.c + ${MPAUDIO_SRC_DIR}/audiopump/audiopump_drift.c + ${MPAUDIO_SRC_DIR}/shared/audiodsp_drift.c ${MPAUDIO_SRC_DIR}/shared/audiodsp_rawsample.c ${MPAUDIO_SRC_DIR}/shared/audiodsp_synth_dsp.c ${MPAUDIO_SRC_DIR}/shared/audiodsp_envelope.c diff --git a/micropython.mk b/micropython.mk index 763a31f3..9835b5ee 100644 --- a/micropython.mk +++ b/micropython.mk @@ -111,7 +111,9 @@ SRC_USERMOD_C += \ $(MPAUDIO_SRC_DIR)/audiopump/audiopump.c \ $(MPAUDIO_SRC_DIR)/audiopump/audiopump_ring.c \ $(MPAUDIO_SRC_DIR)/audiopump/audiopump_events.c \ - $(MPAUDIO_SRC_DIR)/audiopump/audiopump_tap.c + $(MPAUDIO_SRC_DIR)/audiopump/audiopump_tap.c \ + $(MPAUDIO_SRC_DIR)/audiopump/audiopump_drift.c \ + $(MPAUDIO_SRC_DIR)/shared/audiodsp_drift.c # --- module skeleton: empty `import ` targets for every top-level # module this tree will grow into (see docs/porting-plan.md). Tiers add diff --git a/setup.py b/setup.py index 00a4d0ef..962e306f 100644 --- a/setup.py +++ b/setup.py @@ -61,6 +61,7 @@ "src/shared/audiodsp_filter_f32.c", "src/shared/audiodsp_shaper.c", "src/shared/audiodsp_samplehold.c", + "src/shared/audiodsp_drift.c", "src/shared/audiodsp_ladder.c", "src/shared/audiodsp_oscillator.c", "src/shared/audiodsp_trig.c", diff --git a/src/audiopump/audiopump.c b/src/audiopump/audiopump.c index ca84e42d..23b806c0 100644 --- a/src/audiopump/audiopump.c +++ b/src/audiopump/audiopump.c @@ -45,6 +45,7 @@ #include "audiopump/audiopump_events.h" #include "audiopump/audiopump_ring.h" #include "audiopump/audiopump_tap.h" +#include "audiopump/audiopump_drift.h" // The same acquire/release pair the push ring uses, and here for the same // reason: `volatile` orders the compiler, not the machine, and a 64-bit store @@ -2381,6 +2382,7 @@ static const mp_rom_map_elem_t audiopump_globals_table[] = { // The push side. Every port, not just esp32: the ring is RAM and a // memcpy, so the unix build is where its correctness is settled. { MP_ROM_QSTR(MP_QSTR_Ring), MP_ROM_PTR(&audiopump_ring_type) }, + { MP_ROM_QSTR(MP_QSTR_DriftResampler), MP_ROM_PTR(&audiopump_drift_type) }, // Sequenced music, and the way out to a meter. Same reasoning: RAM, // arithmetic and a memcpy, so the unix build settles both. { MP_ROM_QSTR(MP_QSTR_Events), MP_ROM_PTR(&audiopump_events_type) }, diff --git a/src/audiopump/audiopump_drift.c b/src/audiopump/audiopump_drift.c new file mode 100644 index 00000000..227ed0dd --- /dev/null +++ b/src/audiopump/audiopump_drift.c @@ -0,0 +1,295 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: Copyright (c) 2026 PyDevices +// +// audiopump.DriftResampler -- plays a Ring filled on one clock at another +// clock's rate, steering its ratio to hold the Ring at a target fill. +// +// The DSP and the control loop are shared/audiodsp_drift.c, the same on every +// runtime. This file is the MicroPython half: it reads the Ring directly, at +// whatever pace the ratio asks for, through the Ring's own consumer calls +// (audiopump_ring.h), so it must be the Ring's only consumer. +// +// It lives in audiopump rather than beside audiospeed.Resampler because +// CircuitPython has its own audiospeed, and a module of ours may not shadow +// one of CircuitPython's; audiopump is ours on all three runtimes. + +#include + +#include "py/obj.h" +#include "py/runtime.h" + +#include "audiocore/__init__.h" +#include "cp_compat/context_manager_helpers.h" +#include "cp_compat/objproperty.h" +#include "shared/audiodsp_drift.h" +#include "shared/audiodsp_hot.h" +#include "shared/audiodsp_pump_lock.h" + +#include "audiopump/audiopump_drift.h" +#include "audiopump/audiopump_ring.h" + +typedef struct { + audiosample_base_t base; + mp_obj_t ring; + audiodsp_drift_t drift; + float *work; + int16_t *scratch; // the frames one block reads from the ring + uint32_t scratch_frames; + int16_t *out[2]; + uint8_t which; +} audiopump_drift_obj_t; + +// --- the pull: this runs on the pump thread -------------------------------- + +static void audiopump_drift_reset_buffer(audiopump_drift_obj_t *self, + bool single_channel_output, uint8_t channel) { + (void)single_channel_output; + (void)channel; + // Keeps the loop's state and what is queued, as the Ring does: play() + // calls this, and a rate estimate is worth keeping across it. + self->which = 0; +} + +static audioio_get_buffer_result_t AUDIODSP_HOT audiopump_drift_get_buffer( + audiopump_drift_obj_t *self, bool single_channel_output, uint8_t channel, + uint8_t **buffer, uint32_t *buffer_length) { + (void)single_channel_output; + (void)channel; + void *ring = MP_OBJ_TO_PTR(self->ring); + const uint32_t frame = self->drift.channels * 2u; + const uint32_t block = self->drift.frames * frame; + int16_t *dst = self->out[self->which]; + self->which ^= 1; + + const uint32_t underruns = self->drift.underruns; + const uint32_t level = audiopump_ring_readable(ring) / frame; + const uint32_t need = audiodsp_drift_plan(&self->drift, level); + if (need == 0u) { + memset(dst, 0, block); + if (self->drift.underruns != underruns) { + audiopump_ring_starve(ring, block); + } + } else { + // In pieces the size of the scratch buffer: the block that starts + // playback also drops what the ring holds beyond the target. + uint32_t left = need; + while (left != 0u) { + const uint32_t n = left < self->scratch_frames ? left : self->scratch_frames; + audiopump_ring_read(ring, (uint8_t *)self->scratch, n * frame); + audiodsp_drift_push(&self->drift, self->scratch, n); + left -= n; + } + audiodsp_drift_render(&self->drift, dst); + } + *buffer = (uint8_t *)dst; + *buffer_length = block; + return GET_BUFFER_MORE_DATA; +} + +static mp_obj_t audiopump_drift_sources(mp_obj_t self_in, uint8_t index) { + audiopump_drift_obj_t *self = MP_OBJ_TO_PTR(self_in); + return index == 0 ? self->ring : MP_OBJ_NULL; +} + +// --- Python ---------------------------------------------------------------- + +static mp_obj_t audiopump_drift_make_new(const mp_obj_type_t *type, + size_t n_args, size_t n_kw, const mp_obj_t *all_args) { + enum { ARG_ring, ARG_rate, ARG_target, ARG_max_ppm, ARG_bandwidth }; + static const mp_arg_t allowed[] = { + { MP_QSTR_ring, MP_ARG_REQUIRED | MP_ARG_OBJ, {} }, + { MP_QSTR_sample_rate, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE} }, + { MP_QSTR_target, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE} }, + { MP_QSTR_max_ppm, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 1000} }, + { MP_QSTR_bandwidth, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE} }, + }; + mp_arg_val_t args[MP_ARRAY_SIZE(allowed)]; + mp_arg_parse_all_kw_array(n_args, n_kw, all_args, MP_ARRAY_SIZE(allowed), + allowed, args); + + mp_obj_t ring_in = args[ARG_ring].u_obj; + if (!mp_obj_is_type(ring_in, &audiopump_ring_type)) { + mp_raise_TypeError(MP_ERROR_TEXT("ring must be an audiopump.Ring")); + } + audiosample_base_t *ring = MP_OBJ_TO_PTR(ring_in); + audiosample_check_for_deinit(ring); + const uint32_t channels = ring->channel_count; + const uint32_t frame = channels * 2u; + const uint32_t frames = ring->max_buffer_length / frame; + const uint32_t capacity_frames = audiopump_ring_capacity(MP_OBJ_TO_PTR(ring_in)) / frame; + const uint32_t in_rate = ring->sample_rate; + const uint32_t out_rate = args[ARG_rate].u_obj == mp_const_none + ? in_rate : (uint32_t)mp_obj_get_int(args[ARG_rate].u_obj); + const mp_int_t max_ppm = args[ARG_max_ppm].u_int; + const float bandwidth = args[ARG_bandwidth].u_obj == mp_const_none + ? 0.01f : (float)mp_obj_get_float(args[ARG_bandwidth].u_obj); + + if (frames < 16u || frames > AUDIODSP_DRIFT_MAX_FRAMES) { + mp_raise_ValueError(MP_ERROR_TEXT("the ring's frames must be 16 to 1024")); + } + if (!audiodsp_drift_rates_ok(in_rate, out_rate)) { + mp_raise_ValueError(MP_ERROR_TEXT( + "sample_rate must be from the ring's rate to twice it")); + } + if (max_ppm < 0 || max_ppm > (mp_int_t)AUDIODSP_DRIFT_MAX_PPM) { + mp_raise_ValueError(MP_ERROR_TEXT("max_ppm must be from 0 to 10000")); + } + if (!(bandwidth >= AUDIODSP_DRIFT_MIN_BANDWIDTH && bandwidth <= AUDIODSP_DRIFT_MAX_BANDWIDTH)) { + mp_raise_ValueError(MP_ERROR_TEXT("bandwidth must be from 0.001 to 0.1 Hz")); + } + const uint32_t least = audiodsp_drift_min_target(frames, in_rate, out_rate, + (uint32_t)max_ppm); + mp_int_t target = args[ARG_target].u_obj == mp_const_none + ? (mp_int_t)audiodsp_drift_default_target(least, capacity_frames) + : mp_obj_get_int(args[ARG_target].u_obj); + if (target < (mp_int_t)least || target > (mp_int_t)audiodsp_drift_max_target(capacity_frames)) { + mp_raise_ValueError(MP_ERROR_TEXT( + "target must cover a block and fit in the ring")); + } + + audiopump_drift_obj_t *self = mp_obj_malloc(audiopump_drift_obj_t, + (const mp_obj_type_t *)&audiopump_drift_type); + self->ring = ring_in; + const uint32_t work_frames = audiodsp_drift_capacity(frames, in_rate, + out_rate, (uint32_t)max_ppm); + self->work = m_new(float, work_frames * channels); + self->scratch_frames = least; + self->scratch = m_new(int16_t, least * channels); + self->out[0] = m_new(int16_t, frames * channels); + self->out[1] = m_new(int16_t, frames * channels); + memset(self->out[0], 0, frames * frame); + memset(self->out[1], 0, frames * frame); + self->which = 0; + audiodsp_drift_init(&self->drift, self->work, channels, frames, in_rate, + out_rate, (uint32_t)target, (uint32_t)max_ppm, bandwidth); + + self->base.sample_rate = out_rate; + self->base.max_buffer_length = frames * frame; + self->base.bits_per_sample = 16; + self->base.channel_count = (uint8_t)channels; + self->base.samples_signed = true; + self->base.single_buffer = false; + return MP_OBJ_FROM_PTR(self); +} + +static mp_obj_t audiopump_drift_deinit(mp_obj_t self_in) { + audiopump_drift_obj_t *self = MP_OBJ_TO_PTR(self_in); + audiodsp_pump_lock_acquire(); + audiosample_mark_deinit(&self->base); + audiodsp_pump_lock_release(); + return mp_const_none; +} +static MP_DEFINE_CONST_FUN_OBJ_1(audiopump_drift_deinit_obj, + audiopump_drift_deinit); + +static mp_obj_t audiopump_drift_reset(mp_obj_t self_in) { + audiopump_drift_obj_t *self = MP_OBJ_TO_PTR(self_in); + audiosample_check_for_deinit(&self->base); + audiodsp_pump_lock_acquire(); + audiodsp_drift_reset(&self->drift); + audiodsp_pump_lock_release(); + return mp_const_none; +} +static MP_DEFINE_CONST_FUN_OBJ_1(audiopump_drift_reset_obj, + audiopump_drift_reset); + +static mp_obj_t audiopump_drift_get_ppm(mp_obj_t self_in) { + audiopump_drift_obj_t *self = MP_OBJ_TO_PTR(self_in); + audiosample_check_for_deinit(&self->base); + return mp_obj_new_float((mp_float_t)audiodsp_drift_ppm(&self->drift)); +} +static MP_DEFINE_CONST_FUN_OBJ_1(audiopump_drift_get_ppm_obj, + audiopump_drift_get_ppm); +static MP_PROPERTY_GETTER(audiopump_drift_ppm_obj, + (mp_obj_t)&audiopump_drift_get_ppm_obj); + +static mp_obj_t audiopump_drift_get_fill(mp_obj_t self_in) { + audiopump_drift_obj_t *self = MP_OBJ_TO_PTR(self_in); + audiosample_check_for_deinit(&self->base); + return mp_obj_new_float((mp_float_t)audiodsp_drift_fill(&self->drift)); +} +static MP_DEFINE_CONST_FUN_OBJ_1(audiopump_drift_get_fill_obj, + audiopump_drift_get_fill); +static MP_PROPERTY_GETTER(audiopump_drift_fill_obj, + (mp_obj_t)&audiopump_drift_get_fill_obj); + +static mp_obj_t audiopump_drift_get_rate(mp_obj_t self_in) { + audiopump_drift_obj_t *self = MP_OBJ_TO_PTR(self_in); + audiosample_check_for_deinit(&self->base); + return mp_obj_new_float((mp_float_t)audiodsp_drift_rate(&self->drift)); +} +static MP_DEFINE_CONST_FUN_OBJ_1(audiopump_drift_get_rate_obj, + audiopump_drift_get_rate); +static MP_PROPERTY_GETTER(audiopump_drift_rate_obj, + (mp_obj_t)&audiopump_drift_get_rate_obj); + +static mp_obj_t audiopump_drift_get_underruns(mp_obj_t self_in) { + audiopump_drift_obj_t *self = MP_OBJ_TO_PTR(self_in); + audiosample_check_for_deinit(&self->base); + return mp_obj_new_int_from_uint(self->drift.underruns); +} +static MP_DEFINE_CONST_FUN_OBJ_1(audiopump_drift_get_underruns_obj, + audiopump_drift_get_underruns); +static MP_PROPERTY_GETTER(audiopump_drift_underruns_obj, + (mp_obj_t)&audiopump_drift_get_underruns_obj); + +static mp_obj_t audiopump_drift_get_target(mp_obj_t self_in) { + audiopump_drift_obj_t *self = MP_OBJ_TO_PTR(self_in); + audiosample_check_for_deinit(&self->base); + return mp_obj_new_int_from_uint(self->drift.target); +} +static MP_DEFINE_CONST_FUN_OBJ_1(audiopump_drift_get_target_obj, + audiopump_drift_get_target); + +static mp_obj_t audiopump_drift_set_target(mp_obj_t self_in, mp_obj_t value) { + audiopump_drift_obj_t *self = MP_OBJ_TO_PTR(self_in); + audiosample_check_for_deinit(&self->base); + const uint32_t frame = self->drift.channels * 2u; + const mp_int_t target = mp_obj_get_int(value); + const uint32_t least = audiodsp_drift_min_target(self->drift.frames, + self->drift.in_rate, self->drift.out_rate, self->drift.max_ppm); + const uint32_t most = audiopump_ring_capacity(MP_OBJ_TO_PTR(self->ring)) / frame; + if (target < (mp_int_t)least || target > (mp_int_t)audiodsp_drift_max_target(most)) { + mp_raise_ValueError(MP_ERROR_TEXT( + "target must cover a block and fit in the ring")); + } + self->drift.target = (uint32_t)target; + return mp_const_none; +} +static MP_DEFINE_CONST_FUN_OBJ_2(audiopump_drift_set_target_obj, + audiopump_drift_set_target); +static MP_PROPERTY_GETSET(audiopump_drift_target_obj, + (mp_obj_t)&audiopump_drift_get_target_obj, + (mp_obj_t)&audiopump_drift_set_target_obj); + +static const mp_rom_map_elem_t audiopump_driftresampler_locals_table[] = { + { MP_ROM_QSTR(MP_QSTR_deinit), MP_ROM_PTR(&audiopump_drift_deinit_obj) }, + { MP_ROM_QSTR(MP_QSTR___enter__), MP_ROM_PTR(&default___enter___obj) }, + { MP_ROM_QSTR(MP_QSTR___exit__), MP_ROM_PTR(&default___exit___obj) }, + { MP_ROM_QSTR(MP_QSTR_reset), MP_ROM_PTR(&audiopump_drift_reset_obj) }, + { MP_ROM_QSTR(MP_QSTR_ppm), MP_ROM_PTR(&audiopump_drift_ppm_obj) }, + { MP_ROM_QSTR(MP_QSTR_fill), MP_ROM_PTR(&audiopump_drift_fill_obj) }, + { MP_ROM_QSTR(MP_QSTR_rate), MP_ROM_PTR(&audiopump_drift_rate_obj) }, + { MP_ROM_QSTR(MP_QSTR_underruns), MP_ROM_PTR(&audiopump_drift_underruns_obj) }, + { MP_ROM_QSTR(MP_QSTR_target), MP_ROM_PTR(&audiopump_drift_target_obj) }, + AUDIOSAMPLE_FIELDS, +}; +static MP_DEFINE_CONST_DICT(audiopump_drift_locals, + audiopump_driftresampler_locals_table); + +static const audiosample_p_t audiopump_drift_proto = { + MP_PROTO_IMPLEMENT(MP_QSTR_protocol_audiosample) + .reset_buffer = (audiosample_reset_buffer_fun)audiopump_drift_reset_buffer, + .get_buffer = (audiosample_get_buffer_fun)audiopump_drift_get_buffer, + .sources = audiopump_drift_sources, +}; + +MP_DEFINE_CONST_OBJ_TYPE( + audiopump_drift_type, + MP_QSTR_DriftResampler, + MP_TYPE_FLAG_HAS_SPECIAL_ACCESSORS, + make_new, audiopump_drift_make_new, + attr, cp_compat_attr, + locals_dict, &audiopump_drift_locals, + protocol, &audiopump_drift_proto + ); diff --git a/src/audiopump/audiopump_drift.h b/src/audiopump/audiopump_drift.h new file mode 100644 index 00000000..e558e01d --- /dev/null +++ b/src/audiopump/audiopump_drift.h @@ -0,0 +1,10 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: Copyright (c) 2026 PyDevices +// +// audiopump.DriftResampler. See audiopump_drift.c. + +#pragma once + +#include "py/obj.h" + +extern const mp_obj_type_t audiopump_drift_type; diff --git a/src/audiopump/audiopump_ring.c b/src/audiopump/audiopump_ring.c index bb394124..7c02e188 100644 --- a/src/audiopump/audiopump_ring.c +++ b/src/audiopump/audiopump_ring.c @@ -196,6 +196,46 @@ static audioio_get_buffer_result_t AUDIODSP_HOT audiopump_ring_get_buffer( return GET_BUFFER_MORE_DATA; } +// The same consumer, for a node that reads at its own pace. See +// audiopump_ring.h. +uint32_t audiopump_ring_capacity(void *ctx) { + return ((audiopump_ring_obj_t *)ctx)->cap; +} + +uint32_t audiopump_ring_readable(void *ctx) { + audiopump_ring_obj_t *self = ctx; + return AUDIOPUMP_RING_LOAD_ACQ(&self->w) - self->r; +} + +void AUDIODSP_HOT audiopump_ring_read(void *ctx, uint8_t *dst, uint32_t len) { + audiopump_ring_obj_t *self = ctx; + const uint32_t r = self->r; + const uint32_t at = self->rpos; + uint32_t first = self->cap - at; + if (first > len) { + first = len; + } + memcpy(dst, self->ring + at, first); + if (len > first) { + memcpy(dst + first, self->ring, len - first); + } + uint64_t digest = self->out_digest; + for (uint32_t i = 0; i < len; i++) { + digest ^= dst[i]; + digest *= AUDIOPUMP_RING_FNV_PRIME; + } + self->out_digest = digest; + self->read += len; + self->rpos = (at + len) % self->cap; + AUDIOPUMP_RING_STORE_REL(&self->r, r + len); +} + +void audiopump_ring_starve(void *ctx, uint32_t len) { + audiopump_ring_obj_t *self = ctx; + self->underruns++; + self->starved += len; +} + // --- the producer side: this runs on the interpreter thread --------------- // The producer side again, for a C producer on another task: each packet of diff --git a/src/audiopump/audiopump_ring.h b/src/audiopump/audiopump_ring.h index fcfa2a8d..68ee121a 100644 --- a/src/audiopump/audiopump_ring.h +++ b/src/audiopump/audiopump_ring.h @@ -24,3 +24,18 @@ extern const mp_obj_type_t audiopump_ring_type; // at a time: while a C producer feeds the Ring, don't also write() to it. #define AUDIOPUMP_RING_HAS_SINK_WRITE (1) uint32_t audiopump_ring_sink_write(void *ctx, const uint8_t *src, uint32_t len); + +// The consumer side for a node that reads a Ring at its own pace instead of +// a block per pull (audiopump.DriftResampler). Call these only from the +// Ring's one consumer, never alongside a pull of the Ring itself. +// +// `audiopump_ring_capacity` is the ring's size in bytes and +// `audiopump_ring_readable` the bytes queued. `audiopump_ring_read` copies +// `len` of them out, whole frames and no more than readable, and gives the +// room back; the digest and the byte count see them as they would a pull. +// `audiopump_ring_starve` counts a block of silence the reader played +// because the ring ran short, as an underrun of the ring. +uint32_t audiopump_ring_capacity(void *ctx); +uint32_t audiopump_ring_readable(void *ctx); +void audiopump_ring_read(void *ctx, uint8_t *dst, uint32_t len); +void audiopump_ring_starve(void *ctx, uint32_t len); diff --git a/src/circuitpython_spike/shared-bindings/audiopump/DriftResampler.c b/src/circuitpython_spike/shared-bindings/audiopump/DriftResampler.c new file mode 100644 index 00000000..98b9ec10 --- /dev/null +++ b/src/circuitpython_spike/shared-bindings/audiopump/DriftResampler.c @@ -0,0 +1,254 @@ +// audiopump.DriftResampler bindings for CircuitPython. +// +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: Copyright (c) 2026 PyDevices + +#include + +#include "shared-bindings/audiopump/DriftResampler.h" +#include "shared-bindings/audiopump/Ring.h" +#include "shared-bindings/audiocore/__init__.h" + +#include "py/objproperty.h" +#include "py/runtime.h" +#include "shared/runtime/context_manager_helpers.h" + +//| class DriftResampler: +//| """Plays a `Ring` that is filled on another clock, holding it at a +//| target fill. +//| +//| Two devices that both run at 48 kHz differ by tens of parts per +//| million, so a Ring filled from one and played by the other slowly +//| empties or overflows. A DriftResampler reads the Ring at a ratio it +//| steers, in steps far finer than a part per million, so the Ring stays +//| at ``target`` frames and the output never hears a gap or a skip. +//| +//| Play the DriftResampler, not the Ring: it is the Ring's only reader:: +//| +//| ring = Ring(sample_rate=48000, channel_count=2, frames=256, capacity=8) +//| drift = DriftResampler(ring) +//| out.play(drift) +//| while True: +//| ring.write(read_from_the_other_clock()) +//| +//| It plays silence until the Ring first reaches ``target``. If the Ring +//| runs dry it plays a block of silence, counts an underrun on itself and +//| on the Ring, and waits for ``target`` again, keeping its rate.""" +//| +//| def __init__( +//| self, +//| ring: Ring, +//| *, +//| sample_rate: Optional[int] = None, +//| target: Optional[int] = None, +//| max_ppm: int = 1000, +//| bandwidth: float = 0.01, +//| ) -> None: +//| """Read ``ring`` and play it at ``sample_rate``. +//| +//| ``sample_rate`` is the output's rate: the Ring's own (the default) +//| up to twice it. ``target`` is the fill to hold, in frames, half the +//| Ring's capacity by default. ``max_ppm`` bounds how far the ratio +//| moves from nominal, up to 10000. ``bandwidth`` is how fast the loop +//| follows the clocks once settled, in hertz, from 0.001 to 0.1; it +//| starts eight times faster and narrows over about a minute. +//| The Ring's block must be 16 to 1024 frames.""" +//| ... +static mp_obj_t audiopump_driftresampler_make_new(const mp_obj_type_t *type, + size_t n_args, size_t n_kw, const mp_obj_t *all_args) { + enum { ARG_ring, ARG_rate, ARG_target, ARG_max_ppm, ARG_bandwidth }; + static const mp_arg_t allowed[] = { + { MP_QSTR_ring, MP_ARG_REQUIRED | MP_ARG_OBJ, {} }, + { MP_QSTR_sample_rate, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE} }, + { MP_QSTR_target, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE} }, + { MP_QSTR_max_ppm, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 1000} }, + { MP_QSTR_bandwidth, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE} }, + }; + mp_arg_val_t args[MP_ARRAY_SIZE(allowed)]; + mp_arg_parse_all_kw_array(n_args, n_kw, all_args, MP_ARRAY_SIZE(allowed), + allowed, args); + + mp_obj_t ring_in = args[ARG_ring].u_obj; + if (!mp_obj_is_type(ring_in, &audiopump_ring_type)) { + mp_raise_TypeError(MP_ERROR_TEXT("ring must be an audiopump.Ring")); + } + audiopump_ring_obj_t *ring = MP_OBJ_TO_PTR(ring_in); + audiosample_check_for_deinit(&ring->base); + const uint32_t frame = ring->frame; + const uint32_t frames = ring->block / frame; + const uint32_t capacity_frames = ring->cap / frame; + const uint32_t in_rate = ring->base.sample_rate; + const uint32_t out_rate = args[ARG_rate].u_obj == mp_const_none + ? in_rate : (uint32_t)mp_obj_get_int(args[ARG_rate].u_obj); + const mp_int_t max_ppm = args[ARG_max_ppm].u_int; + const float bandwidth = args[ARG_bandwidth].u_obj == mp_const_none + ? 0.01f : (float)mp_obj_get_float(args[ARG_bandwidth].u_obj); + + if (frames < 16u || frames > AUDIODSP_DRIFT_MAX_FRAMES) { + mp_raise_ValueError(MP_ERROR_TEXT("the ring's frames must be 16 to 1024")); + } + if (!audiodsp_drift_rates_ok(in_rate, out_rate)) { + mp_raise_ValueError(MP_ERROR_TEXT( + "sample_rate must be from the ring's rate to twice it")); + } + if (max_ppm < 0 || max_ppm > (mp_int_t)AUDIODSP_DRIFT_MAX_PPM) { + mp_raise_ValueError(MP_ERROR_TEXT("max_ppm must be from 0 to 10000")); + } + if (!(bandwidth >= AUDIODSP_DRIFT_MIN_BANDWIDTH && bandwidth <= AUDIODSP_DRIFT_MAX_BANDWIDTH)) { + mp_raise_ValueError(MP_ERROR_TEXT("bandwidth must be from 0.001 to 0.1 Hz")); + } + const uint32_t least = audiodsp_drift_min_target(frames, in_rate, out_rate, + (uint32_t)max_ppm); + const mp_int_t target = args[ARG_target].u_obj == mp_const_none + ? (mp_int_t)audiodsp_drift_default_target(least, capacity_frames) + : mp_obj_get_int(args[ARG_target].u_obj); + if (target < (mp_int_t)least || target > (mp_int_t)audiodsp_drift_max_target(capacity_frames)) { + mp_raise_ValueError(MP_ERROR_TEXT( + "target must cover a block and fit in the ring")); + } + + audiopump_driftresampler_obj_t *self = mp_obj_malloc( + audiopump_driftresampler_obj_t, &audiopump_driftresampler_type); + audiopump_driftresampler_construct(self, ring, out_rate, (uint32_t)target, + (uint32_t)max_ppm, bandwidth); + return MP_OBJ_FROM_PTR(self); +} + +//| def deinit(self) -> None: +//| """Mark it deinitialised. An output still playing it stops at its +//| next pull. The Ring is left as it is.""" +//| ... +static mp_obj_t audiopump_driftresampler_obj_deinit(mp_obj_t self_in) { + audiopump_driftresampler_obj_t *self = MP_OBJ_TO_PTR(self_in); + audiopump_driftresampler_deinit(self); + return mp_const_none; +} +static MP_DEFINE_CONST_FUN_OBJ_1(audiopump_driftresampler_deinit_obj, + audiopump_driftresampler_obj_deinit); + +//| def reset(self) -> None: +//| """Start over: silence until the Ring reaches ``target``, and the +//| rate estimate forgotten. For when the clock feeding the Ring +//| changes.""" +//| ... +static mp_obj_t audiopump_driftresampler_obj_reset(mp_obj_t self_in) { + audiopump_driftresampler_obj_t *self = MP_OBJ_TO_PTR(self_in); + audiosample_check_for_deinit(&self->base); + audiodsp_drift_reset(&self->drift); + return mp_const_none; +} +static MP_DEFINE_CONST_FUN_OBJ_1(audiopump_driftresampler_reset_obj, + audiopump_driftresampler_obj_reset); + +//| ppm: float +//| """How far the ratio is from nominal, in parts per million: positive +//| when the Ring's clock runs fast against the output's. Read-only.""" +static mp_obj_t audiopump_driftresampler_get_ppm(mp_obj_t self_in) { + audiopump_driftresampler_obj_t *self = MP_OBJ_TO_PTR(self_in); + audiosample_check_for_deinit(&self->base); + return mp_obj_new_float((mp_float_t)audiodsp_drift_ppm(&self->drift)); +} +MP_DEFINE_CONST_FUN_OBJ_1(audiopump_driftresampler_get_ppm_obj, + audiopump_driftresampler_get_ppm); +MP_PROPERTY_GETTER(audiopump_driftresampler_ppm_obj, + (mp_obj_t)&audiopump_driftresampler_get_ppm_obj); + +//| fill: float +//| """The smoothed fill the loop is holding at ``target``, in frames. +//| Read-only.""" +static mp_obj_t audiopump_driftresampler_get_fill(mp_obj_t self_in) { + audiopump_driftresampler_obj_t *self = MP_OBJ_TO_PTR(self_in); + audiosample_check_for_deinit(&self->base); + return mp_obj_new_float((mp_float_t)audiodsp_drift_fill(&self->drift)); +} +MP_DEFINE_CONST_FUN_OBJ_1(audiopump_driftresampler_get_fill_obj, + audiopump_driftresampler_get_fill); +MP_PROPERTY_GETTER(audiopump_driftresampler_fill_obj, + (mp_obj_t)&audiopump_driftresampler_get_fill_obj); + +//| rate: float +//| """The ratio in use: Ring frames read per frame played. Read-only.""" +static mp_obj_t audiopump_driftresampler_get_rate(mp_obj_t self_in) { + audiopump_driftresampler_obj_t *self = MP_OBJ_TO_PTR(self_in); + audiosample_check_for_deinit(&self->base); + return mp_obj_new_float((mp_float_t)audiodsp_drift_rate(&self->drift)); +} +MP_DEFINE_CONST_FUN_OBJ_1(audiopump_driftresampler_get_rate_obj, + audiopump_driftresampler_get_rate); +MP_PROPERTY_GETTER(audiopump_driftresampler_rate_obj, + (mp_obj_t)&audiopump_driftresampler_get_rate_obj); + +//| underruns: int +//| """Blocks of silence played because the Ring ran dry. Read-only.""" +static mp_obj_t audiopump_driftresampler_get_underruns(mp_obj_t self_in) { + audiopump_driftresampler_obj_t *self = MP_OBJ_TO_PTR(self_in); + audiosample_check_for_deinit(&self->base); + return mp_obj_new_int_from_uint(self->drift.underruns); +} +MP_DEFINE_CONST_FUN_OBJ_1(audiopump_driftresampler_get_underruns_obj, + audiopump_driftresampler_get_underruns); +MP_PROPERTY_GETTER(audiopump_driftresampler_underruns_obj, + (mp_obj_t)&audiopump_driftresampler_get_underruns_obj); + +//| target: int +//| """The fill to hold, in Ring frames. Setting it moves the fill there +//| slowly, at the loop's pace, without a jump.""" +//| +static mp_obj_t audiopump_driftresampler_get_target(mp_obj_t self_in) { + audiopump_driftresampler_obj_t *self = MP_OBJ_TO_PTR(self_in); + audiosample_check_for_deinit(&self->base); + return mp_obj_new_int_from_uint(self->drift.target); +} +MP_DEFINE_CONST_FUN_OBJ_1(audiopump_driftresampler_get_target_obj, + audiopump_driftresampler_get_target); + +static mp_obj_t audiopump_driftresampler_set_target(mp_obj_t self_in, + mp_obj_t value) { + audiopump_driftresampler_obj_t *self = MP_OBJ_TO_PTR(self_in); + audiosample_check_for_deinit(&self->base); + const mp_int_t target = mp_obj_get_int(value); + const uint32_t least = audiodsp_drift_min_target(self->drift.frames, + self->drift.in_rate, self->drift.out_rate, self->drift.max_ppm); + const uint32_t most = self->ring->cap / self->ring->frame; + if (target < (mp_int_t)least || target > (mp_int_t)audiodsp_drift_max_target(most)) { + mp_raise_ValueError(MP_ERROR_TEXT( + "target must cover a block and fit in the ring")); + } + self->drift.target = (uint32_t)target; + return mp_const_none; +} +MP_DEFINE_CONST_FUN_OBJ_2(audiopump_driftresampler_set_target_obj, + audiopump_driftresampler_set_target); +MP_PROPERTY_GETSET(audiopump_driftresampler_target_obj, + (mp_obj_t)&audiopump_driftresampler_get_target_obj, + (mp_obj_t)&audiopump_driftresampler_set_target_obj); + +static const mp_rom_map_elem_t audiopump_driftresampler_locals_dict_table[] = { + { MP_ROM_QSTR(MP_QSTR_deinit), MP_ROM_PTR(&audiopump_driftresampler_deinit_obj) }, + { MP_ROM_QSTR(MP_QSTR___enter__), MP_ROM_PTR(&default___enter___obj) }, + { MP_ROM_QSTR(MP_QSTR___exit__), MP_ROM_PTR(&default___exit___obj) }, + { MP_ROM_QSTR(MP_QSTR_reset), MP_ROM_PTR(&audiopump_driftresampler_reset_obj) }, + { MP_ROM_QSTR(MP_QSTR_ppm), MP_ROM_PTR(&audiopump_driftresampler_ppm_obj) }, + { MP_ROM_QSTR(MP_QSTR_fill), MP_ROM_PTR(&audiopump_driftresampler_fill_obj) }, + { MP_ROM_QSTR(MP_QSTR_rate), MP_ROM_PTR(&audiopump_driftresampler_rate_obj) }, + { MP_ROM_QSTR(MP_QSTR_underruns), MP_ROM_PTR(&audiopump_driftresampler_underruns_obj) }, + { MP_ROM_QSTR(MP_QSTR_target), MP_ROM_PTR(&audiopump_driftresampler_target_obj) }, + AUDIOSAMPLE_FIELDS, +}; +static MP_DEFINE_CONST_DICT(audiopump_driftresampler_locals_dict, + audiopump_driftresampler_locals_dict_table); + +static const audiosample_p_t audiopump_driftresampler_proto = { + MP_PROTO_IMPLEMENT(MP_QSTR_protocol_audiosample) + .reset_buffer = (audiosample_reset_buffer_fun)audiopump_driftresampler_reset_buffer, + .get_buffer = (audiosample_get_buffer_fun)audiopump_driftresampler_get_buffer, +}; + +MP_DEFINE_CONST_OBJ_TYPE( + audiopump_driftresampler_type, + MP_QSTR_DriftResampler, + MP_TYPE_FLAG_HAS_SPECIAL_ACCESSORS, + make_new, audiopump_driftresampler_make_new, + locals_dict, &audiopump_driftresampler_locals_dict, + protocol, &audiopump_driftresampler_proto + ); diff --git a/src/circuitpython_spike/shared-bindings/audiopump/DriftResampler.h b/src/circuitpython_spike/shared-bindings/audiopump/DriftResampler.h new file mode 100644 index 00000000..4a7965ec --- /dev/null +++ b/src/circuitpython_spike/shared-bindings/audiopump/DriftResampler.h @@ -0,0 +1,7 @@ +// SPDX-License-Identifier: MIT + +#pragma once + +#include "shared-module/audiopump/DriftResampler.h" + +extern const mp_obj_type_t audiopump_driftresampler_type; diff --git a/src/circuitpython_spike/shared-bindings/audiopump/__init__.c b/src/circuitpython_spike/shared-bindings/audiopump/__init__.c index 429456f1..34c70fa7 100644 --- a/src/circuitpython_spike/shared-bindings/audiopump/__init__.c +++ b/src/circuitpython_spike/shared-bindings/audiopump/__init__.c @@ -7,6 +7,7 @@ #include "py/runtime.h" #include "shared-bindings/audiopump/__init__.h" +#include "shared-bindings/audiopump/DriftResampler.h" #include "shared-bindings/audiopump/Ring.h" //| """A stream sample your code fills, played by CircuitPython's audio outputs @@ -20,11 +21,16 @@ //| On MicroPython, audiodsp's `audiopump` also runs the pull loop itself. //| CircuitPython's outputs already pull, so here the module carries only //| `Ring`, and code written against `audiopump.Ring` runs on both. +//| +//| `DriftResampler` plays a `Ring` filled on another clock (a USB or network +//| stream, a second codec) at the output's rate, steering its ratio so the +//| Ring neither runs dry nor overflows. //| """ static const mp_rom_map_elem_t audiopump_module_globals_table[] = { { MP_ROM_QSTR(MP_QSTR___name__), MP_ROM_QSTR(MP_QSTR_audiopump) }, { MP_ROM_QSTR(MP_QSTR_Ring), MP_ROM_PTR(&audiopump_ring_type) }, + { MP_ROM_QSTR(MP_QSTR_DriftResampler), MP_ROM_PTR(&audiopump_driftresampler_type) }, }; static MP_DEFINE_CONST_DICT(audiopump_module_globals, audiopump_module_globals_table); diff --git a/src/circuitpython_spike/shared-module/audiopump/DriftResampler.c b/src/circuitpython_spike/shared-module/audiopump/DriftResampler.c new file mode 100644 index 00000000..d2c789e8 --- /dev/null +++ b/src/circuitpython_spike/shared-module/audiopump/DriftResampler.c @@ -0,0 +1,123 @@ +// audiopump.DriftResampler for CircuitPython. See DriftResampler.h. +// +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: Copyright (c) 2026 PyDevices + +#include + +#include "py/runtime.h" + +#include "shared-module/audiopump/DriftResampler.h" +#include "shared-bindings/audiocore/__init__.h" + +void audiopump_driftresampler_construct(audiopump_driftresampler_obj_t *self, + audiopump_ring_obj_t *ring, uint32_t out_rate, uint32_t target, + uint32_t max_ppm, float bandwidth) { + const uint32_t channels = ring->base.channel_count; + const uint32_t frame = channels * 2u; + const uint32_t frames = ring->block / frame; + const uint32_t in_rate = ring->base.sample_rate; + const uint32_t least = audiodsp_drift_min_target(frames, in_rate, out_rate, + max_ppm); + const uint32_t work_frames = audiodsp_drift_capacity(frames, in_rate, + out_rate, max_ppm); + + self->ring = ring; + self->work = m_malloc(sizeof(float) * work_frames * channels); + self->scratch_frames = least; + self->scratch = m_malloc(sizeof(int16_t) * least * channels); + self->out[0] = m_malloc(frames * frame); + self->out[1] = m_malloc(frames * frame); + memset(self->out[0], 0, frames * frame); + memset(self->out[1], 0, frames * frame); + self->which = 0; + self->last = 0; + self->pulls = 0; + self->left_reads = 0; + self->right_reads = 0; + audiodsp_drift_init(&self->drift, self->work, channels, frames, in_rate, + out_rate, target, max_ppm, bandwidth); + + self->base.sample_rate = out_rate; + self->base.max_buffer_length = frames * frame; + self->base.bits_per_sample = 16; + self->base.channel_count = (uint8_t)channels; + self->base.samples_signed = true; + self->base.single_buffer = false; +} + +void audiopump_driftresampler_deinit(audiopump_driftresampler_obj_t *self) { + audiosample_mark_deinit(&self->base); +} + +void audiopump_driftresampler_reset_buffer(audiopump_driftresampler_obj_t *self, + bool single_channel_output, uint8_t channel) { + if (single_channel_output && channel == 1) { + return; + } + // Keeps the loop's state and what is queued, as the Ring does: play() + // calls this, and a rate estimate is worth keeping across it. + self->which = 0; + self->left_reads = self->pulls; + self->right_reads = self->pulls; +} + +static int16_t *drift_pull_block(audiopump_driftresampler_obj_t *self) { + audiopump_ring_obj_t *ring = self->ring; + const uint32_t frame = self->drift.channels * 2u; + const uint32_t block = self->drift.frames * frame; + const uint8_t index = self->which; + int16_t *dst = self->out[index]; + self->which ^= 1; + self->last = index; + + const uint32_t underruns = self->drift.underruns; + const uint32_t level = audiopump_ring_readable(ring) / frame; + const uint32_t need = audiodsp_drift_plan(&self->drift, level); + if (need == 0u) { + memset(dst, 0, block); + if (self->drift.underruns != underruns) { + audiopump_ring_starve(ring, block); + } + } else { + // In pieces the size of the scratch buffer: the block that starts + // playback also drops what the ring holds beyond the target. + uint32_t left = need; + while (left != 0u) { + const uint32_t n = left < self->scratch_frames ? left : self->scratch_frames; + audiopump_ring_read(ring, (uint8_t *)self->scratch, n * frame); + audiodsp_drift_push(&self->drift, self->scratch, n); + left -= n; + } + audiodsp_drift_render(&self->drift, dst); + } + self->pulls++; + return dst; +} + +audioio_get_buffer_result_t audiopump_driftresampler_get_buffer( + audiopump_driftresampler_obj_t *self, bool single_channel_output, + uint8_t channel, uint8_t **buffer, uint32_t *buffer_length) { + if (!single_channel_output) { + channel = 0; + } + const uint32_t reads = channel == 1 ? self->right_reads : self->left_reads; + uint8_t *dst; + if (reads == self->pulls) { + dst = (uint8_t *)drift_pull_block(self); + } else { + dst = (uint8_t *)self->out[self->last]; + } + if (channel == 1) { + self->right_reads++; + if (self->drift.channels == 2u) { + // The right channel of interleaved stereo starts one sample in. + dst += 2; + } + } else { + self->left_reads++; + } + *buffer = dst; + *buffer_length = self->drift.frames * self->drift.channels * 2u; + return GET_BUFFER_MORE_DATA; +} diff --git a/src/circuitpython_spike/shared-module/audiopump/DriftResampler.h b/src/circuitpython_spike/shared-module/audiopump/DriftResampler.h new file mode 100644 index 00000000..778be1e0 --- /dev/null +++ b/src/circuitpython_spike/shared-module/audiopump/DriftResampler.h @@ -0,0 +1,54 @@ +// audiopump.DriftResampler for CircuitPython: plays a Ring filled on one +// clock at another clock's rate, steering its ratio to hold the Ring at a +// target fill. +// +// The DSP and the control loop are shared/audiodsp_drift.c, byte for byte the +// code MicroPython and CPython run. This half reads the Ring through its +// consumer calls (Ring.h), at whatever pace the ratio asks for, so it must be +// the Ring's only consumer: play the DriftResampler, never the Ring itself. +// +// It lives in audiopump rather than beside audiospeed.Resampler because +// audiospeed is CircuitPython's own module, and a module of ours never +// shadows one of CircuitPython's. +// +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: Copyright (c) 2026 PyDevices + +#pragma once + +#include +#include + +#include "py/obj.h" +#include "shared-module/audiocore/__init__.h" +#include "shared-module/audiopump/Ring.h" +#include "shared/audiodsp_drift.h" + +typedef struct { + audiosample_base_t base; + audiopump_ring_obj_t *ring; + audiodsp_drift_t drift; + float *work; + int16_t *scratch; // the frames one block reads from the ring + uint32_t scratch_frames; + int16_t *out[2]; + uint8_t which; // the next of out[] to fill + uint8_t last; // the one handed out most recently + // An output that plays the two channels from separate pulls asks once + // per channel per block; the second is handed the same block (Ring.h). + uint32_t pulls; + uint32_t left_reads; + uint32_t right_reads; +} audiopump_driftresampler_obj_t; + +// The arguments are checked by the binding. +void audiopump_driftresampler_construct(audiopump_driftresampler_obj_t *self, + audiopump_ring_obj_t *ring, uint32_t out_rate, uint32_t target, + uint32_t max_ppm, float bandwidth); +void audiopump_driftresampler_deinit(audiopump_driftresampler_obj_t *self); + +void audiopump_driftresampler_reset_buffer(audiopump_driftresampler_obj_t *self, + bool single_channel_output, uint8_t channel); +audioio_get_buffer_result_t audiopump_driftresampler_get_buffer( + audiopump_driftresampler_obj_t *self, bool single_channel_output, + uint8_t channel, uint8_t **buffer, uint32_t *buffer_length); diff --git a/src/circuitpython_spike/shared-module/audiopump/Ring.c b/src/circuitpython_spike/shared-module/audiopump/Ring.c index 5f047317..98b87034 100644 --- a/src/circuitpython_spike/shared-module/audiopump/Ring.c +++ b/src/circuitpython_spike/shared-module/audiopump/Ring.c @@ -155,6 +155,39 @@ static void ring_take_clear(audiopump_ring_obj_t *self) { RING_STORE_REL(&self->clear_done, seq); } +uint32_t audiopump_ring_readable(audiopump_ring_obj_t *self) { + ring_take_clear(self); + return RING_LOAD_ACQ(&self->w) - self->r; +} + +void audiopump_ring_read(audiopump_ring_obj_t *self, uint8_t *dst, + uint32_t len) { + const uint32_t r = self->r; + const uint32_t at = self->rpos; + uint32_t first = self->cap - at; + if (first > len) { + first = len; + } + memcpy(dst, self->ring + at, first); + if (len > first) { + memcpy(dst + first, self->ring, len - first); + } + uint64_t digest = self->out_digest; + for (uint32_t i = 0; i < len; i++) { + digest ^= dst[i]; + digest *= RING_FNV_PRIME; + } + self->out_digest = digest; + self->read += len; + self->rpos = (at + len) % self->cap; + RING_STORE_REL(&self->r, r + len); +} + +void audiopump_ring_starve(audiopump_ring_obj_t *self, uint32_t len) { + self->underruns++; + self->starved += len; +} + static uint8_t *ring_pull_block(audiopump_ring_obj_t *self) { ring_take_clear(self); diff --git a/src/circuitpython_spike/shared-module/audiopump/Ring.h b/src/circuitpython_spike/shared-module/audiopump/Ring.h index 223de391..026b7822 100644 --- a/src/circuitpython_spike/shared-module/audiopump/Ring.h +++ b/src/circuitpython_spike/shared-module/audiopump/Ring.h @@ -98,3 +98,13 @@ void audiopump_ring_reset_buffer(audiopump_ring_obj_t *self, audioio_get_buffer_result_t audiopump_ring_get_buffer( audiopump_ring_obj_t *self, bool single_channel_output, uint8_t channel, uint8_t **buffer, uint32_t *buffer_length); + +// The consumer side for a node that reads the ring at its own pace instead +// of a block per pull (audiopump.DriftResampler). Only from the ring's one +// consumer. `readable` carries out a pending clear() first, as a pull does. +uint32_t audiopump_ring_readable(audiopump_ring_obj_t *self); +// `len` bytes, whole frames, no more than readable. +void audiopump_ring_read(audiopump_ring_obj_t *self, uint8_t *dst, + uint32_t len); +// A block of silence the reader played because the ring ran short. +void audiopump_ring_starve(audiopump_ring_obj_t *self, uint32_t len); diff --git a/src/cpython/_audiodsp.c b/src/cpython/_audiodsp.c index 4478e7a7..08acf884 100644 --- a/src/cpython/_audiodsp.c +++ b/src/cpython/_audiodsp.c @@ -28,6 +28,7 @@ #include "shared/audiodsp_feedback_delay.h" #include "shared/audiodsp_shaper.h" #include "shared/audiodsp_samplehold.h" +#include "shared/audiodsp_drift.h" #include "shared/audiodsp_ladder.h" #include "shared/audiodsp_oscillator.h" #include "shared/audiodsp_tank.h" @@ -121,6 +122,7 @@ typedef struct { PyObject *feedback_delay_state_type; PyObject *waveshaper_state_type; PyObject *samplehold_state_type; + PyObject *drift_state_type; PyObject *ladder_state_type; PyObject *oscillator_state_type; PyObject *biquad_f32_state_type; @@ -1899,6 +1901,178 @@ static PyType_Spec samplehold_state_spec = { .slots = samplehold_state_slots, }; +// audiopump.DriftResampler's kernel: the interpolator's history and the +// control loop. The Python class reads the Ring and hands this the frames +// plan() asked for; everything that decides a sample or a ratio is here, in +// the same C the other two runtimes run. + +typedef struct { + PyObject_HEAD + audiodsp_drift_t drift; + float *work; +} audiodsp_drift_object_t; + +static int drift_state_init(audiodsp_drift_object_t *self, PyObject *args, + PyObject *kwargs) { + unsigned int channels, frames, in_rate, out_rate, target, max_ppm; + unsigned int capacity; + double bandwidth; + static char *keywords[] = {"channel_count", "frames", "in_rate", + "out_rate", "target", "max_ppm", "bandwidth", "capacity", NULL}; + if (!PyArg_ParseTupleAndKeywords(args, kwargs, "IIIIIIdI:DriftState", + keywords, &channels, &frames, &in_rate, &out_rate, &target, &max_ppm, + &bandwidth, &capacity)) return -1; + if (channels < 1 || channels > 2) { + PyErr_SetString(PyExc_ValueError, "channel_count must be 1 or 2"); + return -1; + } + if (frames < 16 || frames > AUDIODSP_DRIFT_MAX_FRAMES) { + PyErr_SetString(PyExc_ValueError, "the ring's frames must be 16 to 1024"); + return -1; + } + if (!audiodsp_drift_rates_ok(in_rate, out_rate)) { + PyErr_SetString(PyExc_ValueError, + "sample_rate must be from the ring's rate to twice it"); + return -1; + } + if (max_ppm > AUDIODSP_DRIFT_MAX_PPM) { + PyErr_SetString(PyExc_ValueError, "max_ppm must be from 0 to 10000"); + return -1; + } + const float bw = (float)bandwidth; + if (!(bw >= AUDIODSP_DRIFT_MIN_BANDWIDTH && bw <= AUDIODSP_DRIFT_MAX_BANDWIDTH)) { + PyErr_SetString(PyExc_ValueError, "bandwidth must be from 0.001 to 0.1 Hz"); + return -1; + } + const uint32_t least = audiodsp_drift_min_target(frames, in_rate, + out_rate, max_ppm); + if (target == 0) { + // None in Python. + target = audiodsp_drift_default_target(least, capacity); + } + if (target < least || target > audiodsp_drift_max_target(capacity)) { + PyErr_SetString(PyExc_ValueError, + "target must cover a block and fit in the ring"); + return -1; + } + const uint32_t work = audiodsp_drift_capacity(frames, in_rate, out_rate, + max_ppm); + PyMem_Free(self->work); + self->work = PyMem_Malloc(sizeof(float) * work * channels); + if (self->work == NULL) { + PyErr_NoMemory(); + return -1; + } + audiodsp_drift_init(&self->drift, self->work, channels, frames, in_rate, + out_rate, target, max_ppm, bw); + return 0; +} + +static void drift_state_dealloc(audiodsp_drift_object_t *self) { + PyTypeObject *type = Py_TYPE(self); + PyMem_Free(self->work); + type->tp_free((PyObject *)self); + Py_DECREF(type); +} + +static PyObject *drift_state_plan(audiodsp_drift_object_t *self, PyObject *arg) { + unsigned long level = PyLong_AsUnsignedLong(arg); + if (level == (unsigned long)-1 && PyErr_Occurred()) return NULL; + if (self->work == NULL) { + PyErr_SetString(PyExc_RuntimeError, "DriftState not initialised"); + return NULL; + } + return PyLong_FromUnsignedLong( + audiodsp_drift_plan(&self->drift, (uint32_t)level)); +} + +static PyObject *drift_state_render(audiodsp_drift_object_t *self, + PyObject *args) { + Py_buffer input = {0}; + if (!PyArg_ParseTuple(args, "y*:render", &input)) return NULL; + const Py_ssize_t frame = (Py_ssize_t)self->drift.channels * 2; + if (self->work == NULL || input.len != (Py_ssize_t)self->drift.pending * frame) { + PyBuffer_Release(&input); + PyErr_SetString(PyExc_ValueError, + "render takes exactly the frames plan() asked for"); + return NULL; + } + PyObject *result = PyBytes_FromStringAndSize(NULL, + (Py_ssize_t)self->drift.frames * frame); + if (result != NULL) { + audiodsp_drift_push(&self->drift, (const int16_t *)input.buf, + self->drift.pending); + audiodsp_drift_render(&self->drift, + (int16_t *)PyBytes_AS_STRING(result)); + } + PyBuffer_Release(&input); + return result; +} + +static PyObject *drift_state_reset(audiodsp_drift_object_t *self, + PyObject *unused) { + (void)unused; + if (self->work != NULL) audiodsp_drift_reset(&self->drift); + Py_RETURN_NONE; +} + +// (ppm, fill, rate, underruns, target) +static PyObject *drift_state_status(audiodsp_drift_object_t *self, + PyObject *unused) { + (void)unused; + return Py_BuildValue("(dddkk)", + (double)audiodsp_drift_ppm(&self->drift), + (double)audiodsp_drift_fill(&self->drift), + (double)audiodsp_drift_rate(&self->drift), + (unsigned long)self->drift.underruns, + (unsigned long)self->drift.target); +} + +static PyObject *drift_state_set_target(audiodsp_drift_object_t *self, + PyObject *args) { + unsigned int target, capacity; + if (!PyArg_ParseTuple(args, "II:set_target", &target, &capacity)) return NULL; + const audiodsp_drift_t *d = &self->drift; + if (target < audiodsp_drift_min_target(d->frames, d->in_rate, d->out_rate, + d->max_ppm) || target > audiodsp_drift_max_target(capacity)) { + PyErr_SetString(PyExc_ValueError, + "target must cover a block and fit in the ring"); + return NULL; + } + self->drift.target = target; + Py_RETURN_NONE; +} + +AUDIODSP_LOCKED_ARGS(drift_state_plan, audiodsp_drift_object_t) +AUDIODSP_LOCKED_ARGS(drift_state_render, audiodsp_drift_object_t) +AUDIODSP_LOCKED_ARGS(drift_state_reset, audiodsp_drift_object_t) +AUDIODSP_LOCKED_ARGS(drift_state_status, audiodsp_drift_object_t) +AUDIODSP_LOCKED_ARGS(drift_state_set_target, audiodsp_drift_object_t) +static PyMethodDef drift_state_methods[] = { + {"plan", (PyCFunction)drift_state_plan_locked, METH_O, NULL}, + {"render", (PyCFunction)drift_state_render_locked, METH_VARARGS, NULL}, + {"reset", (PyCFunction)drift_state_reset_locked, METH_NOARGS, NULL}, + {"status", (PyCFunction)drift_state_status_locked, METH_NOARGS, NULL}, + {"set_target", (PyCFunction)drift_state_set_target_locked, METH_VARARGS, NULL}, + {NULL, NULL, 0, NULL}, +}; + +AUDIODSP_LOCKED_INIT(drift_state_init, audiodsp_drift_object_t) +static PyType_Slot drift_state_slots[] = { + {Py_tp_new, PyType_GenericNew}, + {Py_tp_init, drift_state_init_locked}, + {Py_tp_dealloc, drift_state_dealloc}, + {Py_tp_methods, drift_state_methods}, + {0, NULL}, +}; + +static PyType_Spec drift_state_spec = { + .name = "_audiodsp.DriftState", + .basicsize = sizeof(audiodsp_drift_object_t), + .flags = Py_TPFLAGS_DEFAULT | Py_TPFLAGS_HEAPTYPE, + .slots = drift_state_slots, +}; + // audioladder.Ladder's integrators and solver history. Small enough to sit in // the object -- unlike FeedbackDelayState above there is no line to allocate, // four floats a stage and three more a channel is the whole of it -- but a @@ -5730,6 +5904,11 @@ static int audiodsp_exec(PyObject *module) { AUDIODSP_SHAPER_FRAMES) < 0) return -1; if (PyModule_AddIntConstant(module, "SHAPER_MAX_OVERSAMPLE", AUDIODSP_SHAPER_MAX_OVERSAMPLE) < 0) return -1; + state->drift_state_type = PyType_FromModuleAndSpec(module, + &drift_state_spec, NULL); + if (state->drift_state_type == NULL) return -1; + if (PyModule_AddObjectRef(module, "DriftState", + state->drift_state_type) < 0) return -1; state->samplehold_state_type = PyType_FromModuleAndSpec(module, &samplehold_state_spec, NULL); if (state->samplehold_state_type == NULL) return -1; @@ -5868,6 +6047,7 @@ static int audiodsp_traverse(PyObject *module, visitproc visit, void *arg) { Py_VISIT(state->feedback_delay_state_type); Py_VISIT(state->waveshaper_state_type); Py_VISIT(state->samplehold_state_type); + Py_VISIT(state->drift_state_type); Py_VISIT(state->ladder_state_type); Py_VISIT(state->oscillator_state_type); Py_VISIT(state->biquad_f32_state_type); @@ -5896,6 +6076,7 @@ static int audiodsp_clear(PyObject *module) { Py_CLEAR(state->feedback_delay_state_type); Py_CLEAR(state->waveshaper_state_type); Py_CLEAR(state->samplehold_state_type); + Py_CLEAR(state->drift_state_type); Py_CLEAR(state->ladder_state_type); Py_CLEAR(state->oscillator_state_type); Py_CLEAR(state->biquad_f32_state_type); diff --git a/src/cpython/audiopump.py b/src/cpython/audiopump.py index 6bc27da8..0fcdd42b 100644 --- a/src/cpython/audiopump.py +++ b/src/cpython/audiopump.py @@ -1138,8 +1138,146 @@ def clear(self): self._rpos = self._wpos + # The consumer side for a node that reads at its own pace instead of a + # block per pull (DriftResampler). The same bookkeeping a pull does. + def _readable(self): + return self._w - self._r + + def _take(self, n): + at = self._rpos + first = min(self._cap - at, n) + out = bytes(self._ring[at:at + first]) + if n > first: + out += bytes(self._ring[:n - first]) + self._out_digest = _fnv(out, self._out_digest) + self._read += n + self._rpos = (at + n) % self._cap + self._r += n + return out + + def _starve(self, n): + self._underruns += 1 + self._starved += n + + +class DriftResampler(_AudioSample): + """Plays a Ring filled on one clock at another clock's rate. + + ring = audiopump.Ring(frames=256, capacity=8) + drift = audiopump.DriftResampler(ring, target=1024) + + Two devices that both run at 48 kHz differ by tens of parts per + million, so a Ring between an input and an output on separate clocks + slowly fills or empties. This node reads the Ring between its samples + with a 32-tap windowed sinc, and a slow PI loop steers the ratio to hold + the Ring at `target` frames. It must be the Ring's only consumer. + + `sample_rate` is the output's rate, from the Ring's own up to twice it. + `target` defaults to half the Ring. `max_ppm` bounds the correction and + `bandwidth` (hertz, 0.001 to 0.1) is how fast the loop follows once it + has locked; it starts eight times faster and narrows over about a minute. + + Until the Ring first reaches `target` the node plays silence. If the + Ring runs short of a block, it plays a block of silence, counts an + `underrun`, and waits for the target again, keeping its rate. `ppm`, + `fill` and `rate` say what the loop is doing. + """ + + def __init__(self, ring, *, sample_rate=None, target=None, max_ppm=1000, + bandwidth=0.01): + if not isinstance(ring, Ring): + raise TypeError("ring must be an audiopump.Ring") + ring._check() + self._ring = ring + self.channel_count = ring.channel_count + self.bits_per_sample = 16 + self.samples_signed = True + self.single_buffer = False + frame = ring._frame + self._frame = frame + self._frames = ring._block // frame + self._capacity = ring._cap // frame + in_rate = ring.sample_rate + self.sample_rate = in_rate if sample_rate is None else int(sample_rate) + max_ppm = int(max_ppm) + if max_ppm < 0 or max_ppm > 10000: + raise ValueError("max_ppm must be from 0 to 10000") + # 0 asks the kernel for its default; a real target is at least 2. + target = 0 if target is None else max(int(target), 1) + self._state = _audiodsp.DriftState( + channel_count=self.channel_count, frames=self._frames, + in_rate=in_rate, out_rate=self.sample_rate, target=int(target), + max_ppm=max_ppm, bandwidth=float(bandwidth), + capacity=self._capacity) + self.max_buffer_length = self._frames * frame + self._silence = bytes(self.max_buffer_length) + self._out = (bytearray(self.max_buffer_length), + bytearray(self.max_buffer_length)) + self._which = 0 + self._deinited = False + + @property + def ppm(self): + """The ratio's offset from nominal, in parts per million.""" + self._check() + return self._state.status()[0] + + @property + def fill(self): + """The Ring's fill as the loop sees it, smoothed, in frames.""" + self._check() + return self._state.status()[1] + + @property + def rate(self): + """Input frames read per output frame.""" + self._check() + return self._state.status()[2] + + @property + def underruns(self): + self._check() + return self._state.status()[3] + + @property + def target(self): + self._check() + return self._state.status()[4] + + @target.setter + def target(self, value): + self._check() + self._state.set_target(int(value), self._capacity) + + def reset(self): + """Forget the rate estimate and wait for the target again.""" + self._check() + self._state.reset() + + def _release(self): + pass + + def _reset_buffer(self, single_channel_output=False, audio_channel=0): + self._which = 0 + + def _get_buffer(self, single_channel_output=False, audio_channel=0): + self._check() + ring = self._ring + frame = self._frame + before = self._state.status()[3] + need = self._state.plan(ring._readable() // frame) + dst = self._out[self._which] + self._which ^= 1 + if need == 0: + dst[:] = self._silence + if self._state.status()[3] != before: + ring._starve(len(dst)) + else: + dst[:] = self._state.render(ring._take(need * frame)) + return GET_BUFFER_MORE_DATA, memoryview(dst) + __all__ = ( - "Events", "Ring", "Tap", "backpressure", "drain", "driver", "events", + "DriftResampler", "Events", "Ring", "Tap", "backpressure", "drain", "driver", "events", "fades", "fault", "info", "join", "lock_reset", "lock_stats", "now", "park", "pull", "reset", "retarget", "running", "service", "shutdown", "spawn", "stop", "tap", "threaded", "unpark", diff --git a/src/shared/audiodsp_drift.c b/src/shared/audiodsp_drift.c new file mode 100644 index 00000000..4bea4895 --- /dev/null +++ b/src/shared/audiodsp_drift.c @@ -0,0 +1,270 @@ +// The drift resampler's kernel. See audiodsp_drift.h for what it is, and +// docs/how-it-works/drift-resampler.md for the numbers behind each choice. +// +// SPDX-License-Identifier: MIT + +#include "shared/audiodsp_fp_contract.h" + +#include "shared/audiodsp_drift.h" +#include "shared/audiodsp_drift_table.h" + +#include + +// Frames of history before the read position, and of look-ahead after it. +#define DRIFT_BEHIND (AUDIODSP_DRIFT_TAPS / 2u - 1u) +#define DRIFT_AHEAD (AUDIODSP_DRIFT_TAPS / 2u) + +// The loop's damping, and the factor its starting bandwidth is above its +// final one. The smoothers sit at four times the loop's natural frequency, +// far enough above it to cost little phase and far enough below the block +// rate to take the ring's packet-sized sawtooth out of the measurement. +#define DRIFT_ZETA 0.8f +#define DRIFT_GEAR 7.0f +#define DRIFT_GEAR_SECONDS 15.0f +#define DRIFT_SMOOTH 4.0f +#define DRIFT_TWO_PI 6.28318531f + +static uint64_t drift_nominal(uint32_t in_rate, uint32_t out_rate) { + return (((uint64_t)in_rate << 32) + out_rate / 2u) / out_rate; +} + +static uint64_t drift_max_step(uint64_t nominal, uint32_t max_ppm) { + return nominal + (nominal / 1000000u) * max_ppm + 1u; +} + +bool audiodsp_drift_rates_ok(uint32_t in_rate, uint32_t out_rate) { + if (in_rate == 0u || out_rate == 0u) { + return false; + } + // From 1/2 (doubling the rate) up to 1 (the same rate). + return in_rate <= out_rate && + (uint64_t)in_rate * AUDIODSP_DRIFT_MIN_RATIO_DEN >= + (uint64_t)out_rate * AUDIODSP_DRIFT_MIN_RATIO_NUM; +} + +uint32_t audiodsp_drift_min_target(uint32_t frames, uint32_t in_rate, + uint32_t out_rate, uint32_t max_ppm) { + uint64_t step = drift_max_step(drift_nominal(in_rate, out_rate), max_ppm); + return (uint32_t)((step * frames) >> 32) + 2u; +} + +uint32_t audiodsp_drift_capacity(uint32_t frames, uint32_t in_rate, + uint32_t out_rate, uint32_t max_ppm) { + return AUDIODSP_DRIFT_TAPS + + audiodsp_drift_min_target(frames, in_rate, out_rate, max_ppm) + 2u; +} + +uint32_t audiodsp_drift_max_target(uint32_t ring_frames) { + return ring_frames > DRIFT_AHEAD ? ring_frames - DRIFT_AHEAD : 0u; +} + +uint32_t audiodsp_drift_default_target(uint32_t least, uint32_t ring_frames) { + const uint32_t half = ring_frames / 2u; + return half < least ? least : half; +} + +void audiodsp_drift_reset(audiodsp_drift_t *d) { + d->gear = DRIFT_GEAR; + d->s1 = (float)d->target; + d->s2 = (float)d->target; + d->integ = 0.0f; + d->dev = 0.0f; + d->step = d->nominal; + d->running = false; + // Silent history behind the first frame, so the filter starts from + // silence rather than from whatever the buffer held. + memset(d->work, 0, sizeof(float) * d->capacity * d->channels); + d->have = DRIFT_BEHIND; + d->idx = DRIFT_BEHIND; + d->frac = 0u; + d->pending = 0u; + d->skip = 0u; +} + +void audiodsp_drift_init(audiodsp_drift_t *d, float *work, uint32_t channels, + uint32_t frames, uint32_t in_rate, uint32_t out_rate, uint32_t target, + uint32_t max_ppm, float bandwidth) { + d->channels = channels; + d->frames = frames; + d->in_rate = in_rate; + d->out_rate = out_rate; + d->max_ppm = max_ppm; + d->nominal = drift_nominal(in_rate, out_rate); + d->max_dev = (float)max_ppm * 1e-6f; + d->bandwidth = bandwidth; + d->dt = (float)frames / (float)out_rate; + // exp(-dt/T) to first order; dt/T is under 1e-3 for any block allowed. + d->gear_decay = 1.0f - d->dt / DRIFT_GEAR_SECONDS; + d->target = target; + d->work = work; + d->capacity = audiodsp_drift_capacity(frames, in_rate, out_rate, max_ppm); + d->underruns = 0u; + audiodsp_drift_reset(d); +} + +static float drift_clamp(float v, float limit) { + if (v > limit) { + return limit; + } + if (v < -limit) { + return -limit; + } + return v; +} + +// Frames buffered ahead of the read position: the ring's, plus what the +// kernel has already taken from it beyond the filter's own look-ahead. The +// fractional part keeps the measure smooth as the position creeps. +static float drift_fill(const audiodsp_drift_t *d, uint32_t level) { + return (float)level + (float)(d->have - d->idx) - + (float)DRIFT_AHEAD - (float)(d->frac >> 8) * (1.0f / 16777216.0f); +} + +static void drift_control(audiodsp_drift_t *d, float fill) { + const float w = DRIFT_TWO_PI * d->bandwidth * (1.0f + d->gear); + d->gear *= d->gear_decay; + float a = DRIFT_SMOOTH * w * d->dt; + if (a > 1.0f) { + a = 1.0f; + } + d->s1 += a * (fill - d->s1); + d->s2 += a * (d->s1 - d->s2); + const float e = d->s2 - (float)d->target; + // A fill error of e frames changes the fill at -e * kp * out_rate frames + // a second through the proportional path; with ki the integrator's rate, + // the loop is e'' + 2 zeta w e' + w^2 e = 0. + const float rate = (float)d->out_rate; + const float kp = 2.0f * DRIFT_ZETA * w / rate; + const float ki = w * w / rate; + d->integ = drift_clamp(d->integ + ki * e * d->dt, d->max_dev); + d->dev = drift_clamp(kp * e + d->integ, d->max_dev); + d->step = d->nominal + (uint64_t)(int64_t)((double)d->dev * (double)d->nominal); +} + +static uint32_t drift_need(const audiodsp_drift_t *d) { + const uint64_t last = (uint64_t)d->frac + (uint64_t)(d->frames - 1u) * d->step; + const uint32_t idx_last = d->idx + (uint32_t)(last >> 32); + return idx_last + DRIFT_AHEAD + 1u - d->have; +} + +uint32_t audiodsp_drift_plan(audiodsp_drift_t *d, uint32_t level) { + d->pending = 0u; + d->skip = 0u; + if (!d->running) { + // The fill counts what is beyond the filter's look-ahead, so the + // ring holds that much more than the target when the fill is there. + const uint32_t start = d->target + DRIFT_AHEAD; + if (level < start) { + return 0u; + } + // Start at the target: drop what is beyond it, and start the + // smoothers there so they do not see a step. + d->running = true; + d->skip = level - start; + level = start; + d->s1 = d->s2 = drift_fill(d, level); + } + drift_control(d, drift_fill(d, level)); + const uint32_t need = drift_need(d); + if (need > level) { + // A whole block of silence, and back to waiting for the target: the + // frames in the ring stay there, and the rate estimate is kept. + d->underruns++; + d->running = false; + d->skip = 0u; + return 0u; + } + d->pending = need + d->skip; + return d->pending; +} + +void audiodsp_drift_push(audiodsp_drift_t *d, const int16_t *in, + uint32_t frames) { + if (frames > d->pending) { + frames = d->pending; + } + d->pending -= frames; + if (d->skip != 0u) { + const uint32_t dropped = frames < d->skip ? frames : d->skip; + d->skip -= dropped; + frames -= dropped; + in += (size_t)dropped * d->channels; + } + float *dst = d->work + (size_t)d->have * d->channels; + const uint32_t n = frames * d->channels; + for (uint32_t i = 0; i < n; i++) { + dst[i] = (float)in[i]; + } + d->have += frames; +} + +static int16_t drift_to_int16(float v) { + v += v >= 0.0f ? 0.5f : -0.5f; + if (v >= 32767.0f) { + return 32767; + } + if (v <= -32768.0f) { + return -32768; + } + return (int16_t)(int32_t)v; +} + +void audiodsp_drift_render(audiodsp_drift_t *d, int16_t *out) { + const uint32_t channels = d->channels; + const uint64_t step = d->step; + uint32_t idx = d->idx; + uint32_t frac = d->frac; + float h[AUDIODSP_DRIFT_TAPS]; + + for (uint32_t o = 0; o < d->frames; o++) { + const uint32_t p = frac >> (32u - AUDIODSP_DRIFT_PHASE_BITS); + const float a = (float)((frac << AUDIODSP_DRIFT_PHASE_BITS) >> 8) * + (1.0f / 16777216.0f); + const float *h0 = audiodsp_drift_table[p]; + const float *h1 = audiodsp_drift_table[p + 1u]; + for (uint32_t k = 0; k < AUDIODSP_DRIFT_TAPS; k++) { + h[k] = h0[k] + a * (h1[k] - h0[k]); + } + const float *x = d->work + (size_t)(idx - DRIFT_BEHIND) * channels; + if (channels == 2u) { + float left = 0.0f; + float right = 0.0f; + for (uint32_t k = 0; k < AUDIODSP_DRIFT_TAPS; k++) { + left += h[k] * x[2u * k]; + right += h[k] * x[2u * k + 1u]; + } + out[2u * o] = drift_to_int16(left); + out[2u * o + 1u] = drift_to_int16(right); + } else { + float acc = 0.0f; + for (uint32_t k = 0; k < AUDIODSP_DRIFT_TAPS; k++) { + acc += h[k] * x[k]; + } + out[o] = drift_to_int16(acc); + } + const uint64_t next = (uint64_t)frac + step; + idx += (uint32_t)(next >> 32); + frac = (uint32_t)next; + } + + // Keep the history the next block's first frame needs, at the front. + const uint32_t drop = idx - DRIFT_BEHIND; + const uint32_t keep = d->have - drop; + memmove(d->work, d->work + (size_t)drop * channels, + sizeof(float) * keep * channels); + d->have = keep; + d->idx = idx - drop; + d->frac = frac; +} + +float audiodsp_drift_ppm(const audiodsp_drift_t *d) { + return d->dev * 1e6f; +} + +float audiodsp_drift_fill(const audiodsp_drift_t *d) { + return d->s2; +} + +float audiodsp_drift_rate(const audiodsp_drift_t *d) { + return (float)((double)d->step * (1.0 / 4294967296.0)); +} diff --git a/src/shared/audiodsp_drift.h b/src/shared/audiodsp_drift.h new file mode 100644 index 00000000..f0ca2677 --- /dev/null +++ b/src/shared/audiodsp_drift.h @@ -0,0 +1,145 @@ +// Runtime-neutral kernel of audiopump.DriftResampler: reads a stream that +// arrives on one clock and plays it on another, steering the conversion +// ratio so the buffer between the two clocks stays at a target fill. +// +// New code -- not a CircuitPython port. Two devices that both say 48 kHz +// differ by tens of parts per million, so a Ring between an input and an +// output on separate clocks slowly fills or empties. `audiospeed.Resampler` +// converts at a fixed ratio and `SpeedChanger`'s 16.16 rate moves in steps of +// about 15 ppm, which is the size of the error being corrected. Here the +// read position is 32.32 fixed point, a step of 0.0002 ppm. +// +// Two parts, both described in docs/how-it-works/drift-resampler.md: +// +// * Interpolation: a 32-tap Kaiser-windowed sinc, tabulated at 128 +// fractional positions with linear interpolation between neighbouring +// rows (audiodsp_drift_table.h, written by tools/design_drift_table.py). +// Per output frame that is 32 coefficient blends and 32 multiply-adds a +// channel, all single-precision float. +// +// * Control: once per block the fill (frames buffered ahead of the read +// position) goes through two one-pole smoothers into a PI loop whose +// output is the ratio's offset from nominal, clamped to +/-max_ppm. The +// loop starts at 8x its bandwidth and narrows to it over about a minute, +// so it locks quickly and then holds still. +// +// The pulling loop is not here: each runtime reaches its Ring differently. +// A block is plan(level) -> read that many frames -> push() -> render(). +// +// SPDX-License-Identifier: MIT + +#pragma once + +#include +#include + +#define AUDIODSP_DRIFT_TAPS 32u +#define AUDIODSP_DRIFT_PHASE_BITS 7u +#define AUDIODSP_DRIFT_PHASES (1u << AUDIODSP_DRIFT_PHASE_BITS) +//: The largest block, in output frames. +#define AUDIODSP_DRIFT_MAX_FRAMES 1024u +//: Bounds on max_ppm and on the nominal ratio (input rate over output rate). +//: The filter's cutoff is set for the input rate, so the node converts up or +//: across, never down: going down would fold everything between the two +//: Nyquist frequencies back into the band. +#define AUDIODSP_DRIFT_MAX_PPM 10000u +#define AUDIODSP_DRIFT_MIN_RATIO_NUM 1u +#define AUDIODSP_DRIFT_MIN_RATIO_DEN 2u +//: Bounds on the loop's settled bandwidth, in hertz. Above about 0.1 Hz the +//: loop starts following the producer's packet-sized steps in the fill, and +//: the ratio swings by hundreds of ppm from one block to the next. +#define AUDIODSP_DRIFT_MIN_BANDWIDTH 0.001f +#define AUDIODSP_DRIFT_MAX_BANDWIDTH 0.1f + +typedef struct { + // Configuration, fixed at init. + uint32_t channels; // 1 or 2 + uint32_t frames; // output frames per block + uint32_t in_rate; // Hz + uint32_t out_rate; // Hz + uint32_t max_ppm; + uint64_t nominal; // input frames per output frame, 32.32 + float max_dev; // ratio offset limit, as a fraction + float bandwidth; // the loop's final natural frequency, Hz + float dt; // seconds per block + float gear_decay; // per-block decay of the start-up gear + + // Settable while running. + uint32_t target; // fill to hold, in input frames + + // Control state. + float gear; // extra bandwidth multiple, 7 -> 0 + float s1; // fill, smoothed once + float s2; // fill, smoothed twice + float integ; // the loop's integrator: ratio offset + float dev; // the ratio offset in use this block + uint64_t step; // the ratio in use, 32.32 + bool running; // false: waiting for the ring to reach target + + // Resampling state. `work` is interleaved float frames; `idx` is the + // frame the read position is in, `frac` how far past it. + float *work; + uint32_t capacity; // frames work holds + uint32_t have; // frames work holds now + uint32_t idx; + uint32_t frac; + uint32_t pending; // frames plan() asked for and push() has not had + uint32_t skip; // of those, how many to drop first + + uint32_t underruns; +} audiodsp_drift_t; + +//: True if the pair of rates is one this node converts between. +bool audiodsp_drift_rates_ok(uint32_t in_rate, uint32_t out_rate); + +//: Frames of `work` a node of this shape needs; allocate that times +//: channels floats. Arguments already checked. +uint32_t audiodsp_drift_capacity(uint32_t frames, uint32_t in_rate, + uint32_t out_rate, uint32_t max_ppm); + +//: The input frames one block can consume at most, which a target must +//: cover for the node ever to start. +uint32_t audiodsp_drift_min_target(uint32_t frames, uint32_t in_rate, + uint32_t out_rate, uint32_t max_ppm); + +//: The largest target a ring of `ring_frames` frames allows: the fill is +//: counted beyond the filter's look-ahead, which the ring must also hold. +uint32_t audiodsp_drift_max_target(uint32_t ring_frames); + +//: The target a node gets when none is given: half the ring, or `least` +//: (audiodsp_drift_min_target) if half the ring is less. +uint32_t audiodsp_drift_default_target(uint32_t least, uint32_t ring_frames); + +//: `work` is `audiodsp_drift_capacity(...) * channels` floats. `bandwidth` +//: is in hertz and positive; the bindings check every argument. +void audiodsp_drift_init(audiodsp_drift_t *d, float *work, uint32_t channels, + uint32_t frames, uint32_t in_rate, uint32_t out_rate, uint32_t target, + uint32_t max_ppm, float bandwidth); + +//: Back to the state init left: waiting, history silent, loop at its +//: start-up gear. The rate estimate is forgotten too. +void audiodsp_drift_reset(audiodsp_drift_t *d); + +//: Called once a block with the frames the ring holds. Returns how many to +//: read and push() before render(), or 0 for a block of silence (waiting for +//: the ring to fill, or an underrun, which is counted and sends the node back +//: to waiting). On the block that starts playback, the count includes what +//: the ring holds beyond the target, which push() drops, so playback starts +//: at the target rather than the loop having to drain the excess. +uint32_t audiodsp_drift_plan(audiodsp_drift_t *d, uint32_t level); + +//: Interleaved 16-bit frames, in order; may be called more than once per +//: block (a ring's two segments) as long as the total is what plan() said. +void audiodsp_drift_push(audiodsp_drift_t *d, const int16_t *in, + uint32_t frames); + +//: Writes one block, `frames * channels` samples. Only after a plan() that +//: returned non-zero and the pushes it asked for. +void audiodsp_drift_render(audiodsp_drift_t *d, int16_t *out); + +//: The ratio's offset from nominal, in ppm, and the smoothed fill, in +//: frames, for display. +float audiodsp_drift_ppm(const audiodsp_drift_t *d); +float audiodsp_drift_fill(const audiodsp_drift_t *d); +//: The ratio in use, input frames per output frame. +float audiodsp_drift_rate(const audiodsp_drift_t *d); diff --git a/src/shared/audiodsp_drift_table.h b/src/shared/audiodsp_drift_table.h new file mode 100644 index 00000000..eb068b4a --- /dev/null +++ b/src/shared/audiodsp_drift_table.h @@ -0,0 +1,1305 @@ +// Generated by tools/design_drift_table.py. Do not edit by hand: +// change the script and run it with --write. +// +// 32-tap Kaiser-windowed sinc (beta 9, cutoff 0.45 of +// the input rate) at 128 + 1 fractional positions, each row +// scaled to sum to 1. +// +// SPDX-License-Identifier: MIT + +#pragma once + +#include "shared/audiodsp_drift.h" + +static const float audiodsp_drift_table[AUDIODSP_DRIFT_PHASES + 1][AUDIODSP_DRIFT_TAPS] = { + { + -0.000105453786f, 0.00030489912f, -0.000615486701f, 0.000918322883f, + -0.000898136932f, 5.53861733e-18f, 0.00253729708f, -0.00754460646f, + 0.0156906731f, -0.0272067897f, 0.0416538678f, -0.0578258894f, + 0.0738556832f, -0.0875297114f, 0.0967472121f, 0.90001893f, + 0.0967472121f, -0.0875297114f, 0.0738556832f, -0.0578258894f, + 0.0416538678f, -0.0272067897f, 0.0156906731f, -0.00754460646f, + 0.00253729708f, 5.53861733e-18f, -0.000898136932f, 0.000918322883f, + -0.000615486701f, 0.00030489912f, -0.000105453786f, 1.73018179e-05f, + }, + { + -0.000104354323f, 0.000300368381f, -0.000601771055f, 0.000885566929f, + -0.000833073107f, -0.000110545298f, 0.00269928761f, -0.00774616515f, + 0.0158890188f, -0.0273141582f, 0.0415222421f, -0.0572237112f, + 0.0724064484f, -0.0844730511f, 0.0894247666f, 0.899944603f, + 0.104140602f, -0.0905700624f, 0.0752780885f, -0.0584029555f, + 0.0417655073f, -0.0270851925f, 0.0154832257f, -0.00733789755f, + 0.00237283553f, 0.000111441834f, -0.000963330793f, 0.000950939604f, + -0.000629039947f, 0.000309326686f, -0.000106509382f, 1.75231889e-05f, + }, + { + -0.000103212587f, 0.000295738748f, -0.000587902847f, 0.000852690951f, + -0.000768172264f, -0.000220144575f, 0.00285874377f, -0.00794250984f, + 0.0160782225f, -0.0274073388f, 0.0413708314f, -0.0565968864f, + 0.0709313005f, -0.0814016834f, 0.0821749046f, 0.899722159f, + 0.111603335f, -0.0935925692f, 0.0766728222f, -0.0589544997f, + 0.0418570116f, -0.0269493554f, 0.0152667286f, -0.00712611247f, + 0.00220597023f, 0.000223729847f, -0.00102862192f, 0.000983398757f, + -0.000642421772f, 0.000313647179f, -0.000107519627f, 1.77341371e-05f, + }, + { + -0.000102030113f, 0.000291014294f, -0.000573891331f, 0.000819713518f, + -0.000703466532f, -0.000328749273f, 0.0030156034f, -0.00813357532f, + 0.0162582491f, -0.0274863727f, 0.0411998294f, -0.0559458807f, + 0.0694311261f, -0.0783171728f, 0.0749992207f, 0.899351537f, + 0.119133741f, -0.096595645f, 0.0780390128f, -0.0594800897f, + 0.0419282205f, -0.0267992616f, 0.0150412368f, -0.00690932432f, + 0.00203676918f, 0.000336812896f, -0.00109397725f, 0.00101568119f, + -0.000655622745f, 0.000317856495f, -0.000108482971f, 1.79341514e-05f, + }, + { + -0.000100808451f, 0.000286199123f, -0.000559745939f, 0.000786653196f, + -0.000638987869f, -0.000436311733f, 0.00316980551f, -0.00831930153f, + 0.0164290667f, -0.0275513008f, 0.0410094336f, -0.0552711524f, + 0.0679068193f, -0.0752210692f, 0.067899242f, 0.898832917f, + 0.126730129f, -0.0995777175f, 0.0793758109f, -0.0599793121f, + 0.0419789702f, -0.0266349018f, 0.0148068052f, -0.00668760994f, + 0.00186530151f, 0.000450638996f, -0.00115936319f, 0.00104776816f, + -0.000668633322f, 0.000321950531f, -0.000109397901f, 1.81227206e-05f, + }, + { + -9.95491428e-05f, 0.00028129734f, -0.000545475923f, 0.000753528206f, + -0.00057476759f, -0.000542785216f, 0.00332129095f, -0.0084996298f, + 0.0165906437f, -0.0276021753f, 0.0407998487f, -0.0545731783f, + 0.066359289f, -0.0721149296f, 0.0608764924f, 0.898166418f, + 0.134390771f, -0.102537185f, 0.0806823447f, -0.0604517609f, + 0.0420091227f, -0.0264562722f, 0.014563499f, -0.0064610485f, + 0.00169163849f, 0.000565155409f, -0.00122474553f, 0.00107964023f, + -0.000681444129f, 0.00032592527f, -0.000110262888f, 1.82993372e-05f, + }, + { + -9.82537313e-05f, 0.000276313018f, -0.000531090656f, 0.000720356766f, + -0.000510836719f, -0.000648123911f, 0.00347000291f, -0.0086745061f, + 0.0167429578f, -0.0276390538f, 0.040571291f, -0.0538524352f, + 0.0647894219f, -0.0690002814f, 0.0539324433f, 0.897352278f, + 0.142113924f, -0.105472468f, 0.0819577947f, -0.0608970374f, + 0.0420185439f, -0.0262633823f, 0.0143113835f, -0.00622972287f, + 0.00151585287f, 0.000680308614f, -0.00129009038f, 0.00111127854f, + -0.000694045739f, 0.000329776638f, -0.000111076421f, 1.84634864e-05f, + }, + { + -9.69237517e-05f, 0.000271250232f, -0.000516599219f, 0.000687156862f, + -0.000447225873f, -0.000752282969f, 0.00361588574f, -0.00884387735f, + 0.0168859903f, -0.0276620016f, 0.0403239802f, -0.0531094186f, + 0.0631981492f, -0.0658786595f, 0.0470685475f, 0.896390736f, + 0.149897799f, -0.108381964f, 0.0832013115f, -0.0613147616f, + 0.0420071036f, -0.0260562375f, 0.0140505321f, -0.00599371828f, + 0.00133801904f, 0.000796044304f, -0.00135536294f, 0.00114266353f, + -0.000706428662f, 0.000333500619f, -0.000111836984f, 1.8614659e-05f, + }, + { + -9.55607466e-05f, 0.000266113027f, -0.000502010866f, 0.000653946481f, + -0.000383965147f, -0.000855218444f, 0.00375888613f, -0.00900769606f, + 0.0170197263f, -0.0276710875f, 0.0400581546f, -0.0523446277f, + 0.0615863726f, -0.0627515838f, 0.0402862132f, 0.89528209f, + 0.157740578f, -0.111264087f, 0.0844120756f, -0.0617045611f, + 0.0419746861f, -0.0258348621f, 0.0137810213f, -0.00575312367f, + 0.00115821289f, 0.000912307412f, -0.00142052874f, 0.00117377599f, + -0.000718583527f, 0.000337093195f, -0.000112543086f, 1.87523437e-05f, + }, + { + -9.41662438e-05f, 0.000260905537f, -0.000487334793f, 0.000620743202f, + -0.000321084226f, -0.000956887554f, 0.00389895239f, -0.00916591566f, + 0.0171441585f, -0.0276663918f, 0.0397740491f, -0.0515585616f, + 0.0599550121f, -0.0596205555f, 0.0335868262f, 0.894026697f, + 0.165640444f, -0.114117227f, 0.0855892673f, -0.0620660707f, + 0.0419211872f, -0.0255992822f, 0.0135029331f, -0.00550802983f, + 0.000976512092f, 0.00102904218f, -0.00148555252f, 0.00120459648f, + -0.000730500906f, 0.00034055041f, -0.000113193237f, 1.88760332e-05f, + }, + { + -9.27417641e-05f, 0.000255631719f, -0.000472579966f, 0.000587564544f, + -0.000258612359f, -0.00105724821f, 0.00403603492f, -0.00931849517f, + 0.0172592811f, -0.0276479963f, 0.0394719094f, -0.050751742f, + 0.0583049878f, -0.0564870797f, 0.0269717295f, 0.892624855f, + 0.173595518f, -0.116939791f, 0.086732097f, -0.0623989403f, + 0.04184651f, -0.0253495276f, 0.013216354f, -0.00525853131f, + 0.000792995619f, 0.00114619231f, -0.00155039912f, 0.00123510568f, + -0.000742171484f, 0.000343868247f, -0.000113785965f, 1.89852162e-05f, + }, + { + -9.12888427e-05f, 0.000250295649f, -0.000457755465f, 0.000554427796f, + -0.000196578272f, -0.00115625968f, 0.00417008577f, -0.00946539361f, + 0.0173650961f, -0.0276159905f, 0.0391519926f, -0.0499246828f, + 0.0566372275f, -0.0533526428f, 0.0204422399f, 0.891077042f, + 0.181603938f, -0.119730175f, 0.0878397599f, -0.0627028272f, + 0.0417505652f, -0.0250856448f, 0.0129213771f, -0.00500472542f, + 0.000607743976f, 0.00126370054f, -0.00161503302f, 0.00126528414f, + -0.000753585773f, 0.000347042835f, -0.00011431981f, 1.90793908e-05f, + }, + { + -8.98089784e-05f, 0.000244901312f, -0.000442870194f, 0.000521350128f, + -0.000135010225f, -0.00125388207f, 0.00430105906f, -0.00960657652f, + 0.0174616054f, -0.027570473f, 0.0388145596f, -0.0490779169f, + 0.054952655f, -0.0502187237f, 0.0139996316f, 0.889383733f, + 0.189663753f, -0.122486785f, 0.0889114738f, -0.0629774183f, + 0.0416332819f, -0.0248076785f, 0.0126180965f, -0.00474671181f, + 0.000420839351f, 0.00138150924f, -0.00167941849f, 0.00129511254f, + -0.000764734636f, 0.000350070244f, -0.000114793322f, 1.91580493e-05f, + }, + { + -8.83036773e-05f, 0.000239452711f, -0.000427933031f, 0.000488348422f, + -7.39359893e-05f, -0.0013500764f, 0.00442891056f, -0.00974201038f, + 0.01754882f, -0.0275115464f, 0.0384598747f, -0.0482119694f, + 0.0532522015f, -0.0470867865f, 0.00764515065f, 0.887545407f, + 0.197773054f, -0.12520802f, 0.0899464786f, -0.0632223859f, + 0.0414945856f, -0.0245156847f, 0.0123066157f, -0.0044845934f, + 0.00023236507f, 0.00149956008f, -0.00174351991f, 0.0013245712f, + -0.0007756087f, 0.000352946634f, -0.000115205075f, 1.92206935e-05f, + }, + { + -8.67744384e-05f, 0.000233953775f, -0.000412952766f, 0.000455439382f, + -1.33828025e-05f, -0.00144480506f, 0.00455359789f, -0.00987166539f, + 0.0176267549f, -0.0274393167f, 0.0380882211f, -0.0473273918f, + 0.051536791f, -0.0439582765f, 0.00138000422f, 0.885562658f, + 0.205929875f, -0.127892286f, 0.090944007f, -0.0634374395f, + 0.0413344279f, -0.0242097266f, 0.0119870389f, -0.00421847589f, + 4.24060891e-05f, 0.00161779416f, -0.00180730107f, 0.00135364069f, + -0.000786198769f, 0.000355668133f, -0.000115553652f, 1.92668249e-05f, + }, + { + -8.52227458e-05f, 0.000228408462f, -0.000397938129f, 0.000422639481f, + 4.66225902e-05f, -0.00153803127f, 0.00467508193f, -0.00999551453f, + 0.0176954288f, -0.0273539033f, 0.0376998782f, -0.0464247204f, + 0.0498073585f, -0.0408346355f, -0.00479463488f, 0.883436024f, + 0.214132234f, -0.130537972f, 0.0919033214f, -0.0636222884f, + 0.0411527604f, -0.0238898769f, 0.0116594769f, -0.00394846732f, + -0.000148951382f, 0.00173615187f, -0.00187072577f, 0.00138230179f, + -0.000796495646f, 0.000358230987f, -0.000115837654f, 1.92959506e-05f, + }, + { + -8.36500767e-05f, 0.000222820672f, -0.000382897735f, 0.000389964989f, + 0.00010605398f, -0.00162971939f, 0.00479332311f, -0.0101135345f, + 0.0177548639f, -0.0272554215f, 0.0372951329f, -0.0455045104f, + 0.0480648354f, -0.0377172902f, -0.0108776297f, 0.88116622f, + 0.22237812f, -0.133143514f, 0.0928236842f, -0.0637766495f, + 0.0409495458f, -0.0235562101f, 0.0113240471f, -0.00367467897f, + -0.000341619743f, 0.00185457338f, -0.0019337578f, 0.00141053495f, + -0.00080649025f, 0.000360631413f, -0.000116055715f, 1.9307583e-05f, + }, + { + -8.20578862e-05f, 0.000217194276f, -0.00036784017f, 0.000357431883f, + 0.0001648857f, -0.00171983475f, 0.00490828604f, -0.0102257049f, + 0.0178050902f, -0.0271440018f, 0.0368742831f, -0.0445673168f, + 0.0463101529f, -0.0346076451f, -0.0168678798f, 0.878753841f, + 0.230665505f, -0.135707304f, 0.0937043875f, -0.0639002696f, + 0.0407247581f, -0.0232088137f, 0.0109808668f, -0.00339722424f, + -0.000535510131f, 0.00197299803f, -0.00199636072f, 0.00143832085f, + -0.00081617356f, 0.000362865743f, -0.000116206465f, 1.93012347e-05f, + }, + { + -8.04476222e-05f, 0.000211533115f, -0.000352773844f, 0.000325055968f, + 0.000223092618f, -0.00180834404f, 0.00501993531f, -0.0103320088f, + 0.0178461354f, -0.0270197745f, 0.0364376307f, -0.0436137021f, + 0.0445442423f, -0.0315071009f, -0.0227643196f, 0.876199663f, + 0.238992333f, -0.138227761f, 0.0945447162f, -0.0639928952f, + 0.040478386f, -0.022847781f, 0.0106300628f, -0.00311621954f, + -0.000730532338f, 0.0020913647f, -0.0020584981f, 0.00146564026f, + -0.000825536554f, 0.00036493034f, -0.000116288582f, 1.92764292e-05f, + }, + { + -7.88207326e-05f, 0.000205841017f, -0.000337707141f, 0.000292852754f, + 0.000280650158f, -0.00189521466f, 0.00512823975f, -0.0104324305f, + 0.0178780425f, -0.0268828776f, 0.035985481f, -0.0426442362f, + 0.0427680351f, -0.0284170341f, -0.0285659209f, 0.87350446f, + 0.247356534f, -0.140703291f, 0.0953439847f, -0.0640542954f, + 0.0402104259f, -0.0224732105f, 0.0102717644f, -0.00283178361f, + -0.000926595123f, 0.0022096124f, -0.00212013349f, 0.00149247376f, + -0.000834570383f, 0.000366821536f, -0.000116300747f, 1.92326916e-05f, + }, + { + -7.71786363e-05f, 0.000200121751f, -0.000322648295f, 0.00026083758f, + 0.000337534351f, -0.00198041555f, 0.00523316814f, -0.0105269607f, + 0.0179008469f, -0.0267334562f, 0.0355181508f, -0.0416594855f, + 0.0409824587f, -0.0253388118f, -0.0342716947f, 0.870669007f, + 0.255756021f, -0.143132344f, 0.096101515f, -0.0640842393f, + 0.0399208851f, -0.0220852084f, 0.00990610477f, -0.00254403776f, + -0.00112360588f, 0.00232767873f, -0.00218122988f, 0.00151880225f, + -0.000843266142f, 0.000368535839f, -0.000116241674f, 1.91695544e-05f, + }, + { + -7.55227302e-05f, 0.000194379041f, -0.000307605485f, 0.000229025405f, + 0.000393721712f, -0.00206391653f, 0.00533469208f, -0.010615591f, + 0.0179145988f, -0.0265716556f, 0.0350359604f, -0.0406600274f, + 0.0391884446f, -0.0222737808f, -0.0398806818f, 0.867694139f, + 0.264188707f, -0.145513341f, 0.096816659f, -0.0640825331f, + 0.0396097749f, -0.0216838922f, 0.00953322276f, -0.00225310563f, + -0.00132147083f, 0.00244550197f, -0.00224175071f, 0.00154460664f, + -0.000851615157f, 0.000370069727f, -0.00011611011f, 1.90865612e-05f, + }, + { + -7.38544186e-05f, 0.000188616585f, -0.000292586745f, 0.000197431029f, + 0.000449189392f, -0.00214568852f, 0.00543278456f, -0.0106983157f, + 0.017919343f, -0.0263976324f, 0.0345392302f, -0.0396464355f, + 0.0373869166f, -0.0192232747f, -0.0453919694f, 0.86458087f, + 0.272652417f, -0.147844717f, 0.0974887609f, -0.0640489832f, + 0.0392771326f, -0.021269381f, 0.00915326085f, -0.0019591141f, + -0.00152009528f, 0.00256301905f, -0.00230165967f, 0.00156986807f, + -0.000859608874f, 0.000371419766f, -0.000115904812f, 1.89832554e-05f, + }, + { + -7.21750694e-05f, 0.000182838063f, -0.000277600047f, 0.000166068945f, + 0.000503915129f, -0.00222570379f, 0.00552742137f, -0.010775133f, + 0.0179151371f, -0.0262115467f, 0.0340282954f, -0.0386192985f, + 0.0355787985f, -0.0161886066f, -0.0508046746f, 0.861329973f, + 0.281145036f, -0.150124922f, 0.0981171876f, -0.0639833957f, + 0.0389229916f, -0.0208418071f, 0.0087663671f, -0.00166219112f, + -0.00171938329f, 0.00268016732f, -0.00236091972f, 0.00159456744f, + -0.000867238618f, 0.00037258261f, -0.000115624578f, 1.88591948e-05f, + }, + { + -7.04860431e-05f, 0.000177047084f, -0.000262653193f, 0.000134953341f, + 0.00055787724f, -0.00230393535f, 0.00561858015f, -0.0108460439f, + 0.017902039f, -0.0260135606f, 0.0335034914f, -0.0375792012f, + 0.0337650143f, -0.0131710758f, -0.0561179556f, 0.857942581f, + 0.289664388f, -0.152352408f, 0.0987013429f, -0.0638856217f, + 0.0385474041f, -0.0204013046f, 0.00837269332f, -0.00136246812f, + -0.00191923801f, 0.00279688323f, -0.00241949433f, 0.00161868602f, + -0.000874496007f, 0.000373554911f, -0.000115268231f, 1.87139394e-05f, + }, + { + -6.87886786e-05f, 0.0001712472f, -0.000247753895f, 0.000104098181f, + 0.000611054653f, -0.00238035759f, 0.00570623949f, -0.0109110521f, + 0.0178801101f, -0.0258038435f, 0.0329651535f, -0.0365267284f, + 0.0319464765f, -0.0101719601f, -0.0613310076f, 0.854419589f, + 0.298208266f, -0.154525638f, 0.0992406234f, -0.0637555122f, + 0.0381504335f, -0.0199480206f, 0.00797239598f, -0.00106007815f, + -0.00211956142f, 0.00291310367f, -0.00247734692f, 0.00164220529f, + -0.000881372776f, 0.000374333496f, -0.00011483462f, 1.85470653e-05f, + }, + { + -6.70843074e-05f, 0.000165441947f, -0.000232909777f, 7.35171197e-05f, + 0.000663426879f, -0.00245494628f, 0.00579038123f, -0.010970165f, + 0.0178494155f, -0.0255825706f, 0.0324136317f, -0.0354624726f, + 0.0301241018f, -0.00719252322f, -0.0664430633f, 0.850762069f, + 0.306774527f, -0.156643078f, 0.0997344479f, -0.0635929182f, + 0.0377321467f, -0.0194821022f, 0.00756563526f, -0.000755156798f, + -0.00232025515f, 0.00302876462f, -0.00253444095f, 0.00166510663f, + -0.000887860719f, 0.000374915136f, -0.00011432264f, 1.83581524e-05f, + }, + { + -6.53742172e-05f, 0.000159634801f, -0.000218128276f, 4.32235101e-05f, + 0.00071497407f, -0.00252767792f, 0.00587098813f, -0.0110233929f, + 0.017810028f, -0.0253499206f, 0.0318492763f, -0.0343870297f, + 0.0282988027f, -0.0042340057f, -0.0714533925f, 0.846971154f, + 0.315360934f, -0.158703223f, 0.100182258f, -0.0633977354f, + 0.0372926332f, -0.0190037116f, 0.00715257553f, -0.00044784174f, + -0.00252121896f, 0.00314380229f, -0.00259073987f, 0.00168737164f, + -0.000893951743f, 0.000375296775f, -0.000113731199f, 1.81467913e-05f, + }, + { + -6.36597033e-05f, 0.000153829198f, -0.000203416756f, 1.32304358e-05f, + 0.000765676901f, -0.00259853038f, 0.00594804529f, -0.0110707479f, + 0.0177620202f, -0.025106078f, 0.0312724374f, -0.0333009884f, + 0.0264714807f, -0.00129763153f, -0.0763613135f, 0.843047917f, + 0.323965251f, -0.160704538f, 0.100583509f, -0.0631698593f, + 0.0368319787f, -0.0185130127f, 0.00673338631f, -0.000138272706f, + -0.00272235274f, 0.0032581524f, -0.00264620734f, 0.00170898228f, + -0.000899637875f, 0.000375475356f, -0.000113059257f, 1.79125836e-05f, + }, + { + -6.19420243e-05f, 0.000148028514f, -0.000188782447f, -1.64493304e-05f, + 0.000815516803f, -0.00266748248f, 0.00602153968f, -0.0111122467f, + 0.0177054685f, -0.0248512272f, 0.0306834765f, -0.0322049521f, + 0.0246430375f, 0.00161539461f, -0.0811661631f, 0.838993609f, + 0.332585275f, -0.162645563f, 0.100937687f, -0.0629091859f, + 0.0363502987f, -0.0180101767f, 0.00630824082f, 0.000173408436f, + -0.00292355497f, 0.0033717507f, -0.00270080729f, 0.00172992051f, + -0.000904911314f, 0.000375447999f, -0.000112305803f, 1.76551457e-05f, + }, + { + -6.02224209e-05f, 0.000142236051f, -0.000174232424f, -4.58033137e-05f, + 0.000864475791f, -0.00273451465f, 0.00609146059f, -0.0111479061f, + 0.0176404584f, -0.0245855656f, 0.0300827548f, -0.031099515f, + 0.0228143707f, 0.00450388901f, -0.085867323f, 0.834809422f, + 0.34121868f, -0.164524779f, 0.101244271f, -0.0626156628f, + 0.0358476974f, -0.0174953863f, 0.00587731553f, 0.00048705787f, + -0.00312472368f, 0.00348453294f, -0.0027545034f, 0.00175016839f, + -0.000909764494f, 0.000375211792f, -0.000111469853f, 1.73740973e-05f, + }, + { + -5.85021153e-05f, 0.000136455114f, -0.000159773655f, -7.48193415e-05f, + 0.00091253646f, -0.00279960805f, 0.00615779823f, -0.0111777494f, + 0.0175670739f, -0.0243092868f, 0.0294706356f, -0.0299852788f, + 0.0209863689f, 0.00736668706f, -0.0904642269f, 0.83049655f, + 0.349863261f, -0.166340724f, 0.101502784f, -0.0622892156f, + 0.0353243127f, -0.0169688277f, 0.00544079253f, 0.000802529801f, + -0.00332575594f, 0.0035964346f, -0.00280726003f, 0.00176970835f, + -0.000914189732f, 0.000374764f, -0.000110550478f, 1.70690764e-05f, + }, + { + -5.67823117e-05f, 0.000130688917f, -0.000145412967f, -0.00010348556f, + 0.000959681987f, -0.0028627452f, 0.00622054609f, -0.0112017998f, + 0.0174854062f, -0.0240225941f, 0.0288474914f, -0.0288628414f, + 0.0191599168f, 0.0102026463f, -0.0949563384f, 0.826056361f, + 0.358516723f, -0.168091953f, 0.101712748f, -0.0619298108f, + 0.0347802751f, -0.0164306946f, 0.00499885576f, 0.00111967651f, + -0.00352654839f, 0.00370739121f, -0.00285904133f, 0.00178852293f, + -0.000918179809f, 0.000374101917f, -0.000109546781f, 1.67397338e-05f, + }, + { + -5.50641853e-05f, 0.000124940605f, -0.000131157023f, -0.00013179041f, + 0.00100589637f, -0.00292390981f, 0.00627969857f, -0.0112200854f, + 0.0173955448f, -0.0237256885f, 0.0282136947f, -0.0277328063f, + 0.0173358899f, 0.0130106425f, -0.0993431434f, 0.821490109f, + 0.367176801f, -0.169777006f, 0.101873733f, -0.0615374185f, + 0.0342157409f, -0.0158811864f, 0.00455169473f, 0.00143834867f, + -0.00372699695f, 0.00381733826f, -0.00290981168f, 0.00180659513f, + -0.000921727449f, 0.000373222952f, -0.000108457898f, 1.63857276e-05f, + }, + { + -5.33488856e-05f, 0.000119213299f, -0.000117012394f, -0.000159722666f, + 0.00105116412f, -0.00298308628f, 0.00633525196f, -0.011232635f, + 0.0172975883f, -0.0234187823f, 0.0275696218f, -0.0265957713f, + 0.0155151617f, 0.0157895759f, -0.103624195f, 0.816799164f, + 0.375841141f, -0.171394452f, 0.101985291f, -0.0611120313f, + 0.0336308628f, -0.0153205106f, 0.00409950269f, 0.001758395f, + -0.00392699707f, 0.00392621104f, -0.00295953569f, 0.00182390772f, + -0.000924825668f, 0.000372124632f, -0.000107283013f, 1.60067339e-05f, + }, + { + -5.16375476e-05f, 0.000113510039f, -0.000102985476f, -0.00018727142f, + 0.00109547027f, -0.0030402611f, 0.00638720486f, -0.0112394821f, + 0.0171916354f, -0.023102086f, 0.026915649f, -0.02545234f, + 0.0136985937f, 0.018538367f, -0.107799068f, 0.811984956f, + 0.384507507f, -0.172942877f, 0.102047034f, -0.060653653f, + 0.0330258198f, -0.0147488844f, 0.003642475f, 0.0020796624f, + -0.00412644353f, 0.0040339455f, -0.00300817867f, 0.00184044405f, + -0.000927467539f, 0.000370804511f, -0.000106021362f, 1.56024416e-05f, + }, + { + -4.9931281e-05f, 0.000107833817f, -8.90825177e-05f, -0.000214426065f, + 0.00113880076f, -0.00309542078f, 0.00643555773f, -0.0112406602f, + 0.0170777924f, -0.0227758158f, 0.0262521598f, -0.0243031103f, + 0.0118870437f, 0.021255957f, -0.111867368f, 0.807048857f, + 0.393173516f, -0.174420878f, 0.102058552f, -0.0601623058f, + 0.0324007869f, -0.0141665265f, 0.00318081258f, 0.00240199687f, + -0.00432523061f, 0.00414047716f, -0.00305570546f, 0.00185618771f, + -0.000929646427f, 0.00036926032f, -0.00010467221f, 1.51725544e-05f, + }, + { + -4.82311661e-05f, 0.000102187572f, -7.53096392e-05f, -0.000241176327f, + 0.00118114206f, -0.00314855389f, 0.00648031244f, -0.0112362094f, + 0.0169561617f, -0.0224401932f, 0.0255795382f, -0.023148682f, + 0.0100813583f, 0.0239413101f, -0.11582876f, 0.801992416f, + 0.401836902f, -0.175827071f, 0.102019496f, -0.0596380271f, + 0.0317559615f, -0.0135736652f, 0.00271471869f, 0.00272524171f, + -0.00452325307f, 0.00424574222f, -0.00310208183f, 0.00187112251f, + -0.000931355753f, 0.000367489847f, -0.000103234866f, 1.47167921e-05f, + }, + { + -4.65382655e-05f, 9.65741565e-05f, -6.16728139e-05f, -0.000267512281f, + 0.00122248125f, -0.00319964997f, 0.00652147317f, -0.0112261679f, + 0.0168268587f, -0.02209544f, 0.0248981696f, -0.0219896547f, + 0.00828237925f, 0.0265934132f, -0.119682938f, 0.796817124f, + 0.410495281f, -0.177160099f, 0.101929516f, -0.059080869f, + 0.0310915448f, -0.0129705351f, 0.00224440033f, 0.00304923998f, + -0.00472040381f, 0.00434967596f, -0.00314727332f, 0.00188523252f, + -0.000932589232f, 0.000365490967f, -0.000101708705f, 1.42348863e-05f, + }, + { + -4.48536084e-05f, 9.09963856e-05f, -4.81778588e-05f, -0.000293424295f, + 0.00126280601f, -0.00324869901f, 0.00655904552f, -0.0112105804f, + 0.0166899916f, -0.0217417832f, 0.0242084414f, -0.0208266266f, + 0.00649093743f, 0.0292112771f, -0.123429619f, 0.79152441f, + 0.419146359f, -0.178418607f, 0.101788282f, -0.0584909096f, + 0.0304077566f, -0.0123573784f, 0.00177006831f, 0.00337383244f, + -0.00491657667f, 0.00445221504f, -0.00319124642f, 0.00189850212f, + -0.000933340809f, 0.000363261701f, -0.000100093122f, 1.37265888e-05f, + }, + { + -4.31782064e-05f, 8.54570026e-05f, -3.48304311e-05f, -0.000318903098f, + 0.00130210491f, -0.00329569331f, 0.00659303693f, -0.0111894915f, + 0.0165456794f, -0.0213794503f, 0.0235107485f, -0.0196601916f, + 0.00470785564f, 0.0317939334f, -0.127068579f, 0.786115944f, + 0.427787781f, -0.179601267f, 0.101595506f, -0.0578682311f, + 0.0297048185f, -0.01173444f, 0.00129193638f, 0.0036988589f, + -0.00511166407f, 0.00455329521f, -0.00323396735f, 0.00191091583f, + -0.000933604431f, 0.000360800157f, -9.83875871e-05f, 1.31916649e-05f, + }, + { + -4.15130453e-05f, 7.99586851e-05f, -2.16360477e-05f, -0.000343939726f, + 0.00134036678f, -0.00334062544f, 0.00662345719f, -0.0111629497f, + 0.0163940396f, -0.0210086759f, 0.0228054784f, -0.0184909441f, + 0.00293394784f, 0.0343404375f, -0.130599633f, 0.780593336f, + 0.436417222f, -0.180706754f, 0.101350896f, -0.0572129339f, + 0.0289829671f, -0.011101976f, 0.000810221303f, 0.00402415777f, + -0.00530556031f, 0.00465285406f, -0.00327540329f, 0.00192245888f, + -0.000933374511f, 0.00035810456f, -9.6591597e-05f, 1.26298974e-05f, + }, + { + -3.98590819e-05f, 7.45040525e-05f, -8.60007094e-06f, -0.000368525536f, + 0.00137758139f, -0.00338348909f, 0.00665031699f, -0.0111310054f, + 0.016235197f, -0.0206296947f, 0.022093026f, -0.0173194781f, + 0.00117001776f, 0.0368498713f, -0.134022623f, 0.774958193f, + 0.445032269f, -0.181733787f, 0.101054192f, -0.0565251447f, + 0.0282424539f, -0.0104602445f, 0.000325143425f, 0.00434956606f, + -0.00549815688f, 0.00475082779f, -0.00331552117f, 0.00193311635f, + -0.000932645402f, 0.000355173252f, -9.47047156e-05f, 1.20410887e-05f, + }, + { + -3.82172475e-05f, 6.90956513e-05f, 4.27229452e-06f, -0.000392652291f, + 0.00141373917f, -0.00342427962f, 0.00667362893f, -0.0110937115f, + 0.0160692763f, -0.0202427432f, 0.021373786f, -0.0161463823f, + -0.000583140645f, 0.0393213332f, -0.137337416f, 0.769212127f, + 0.453630656f, -0.182681099f, 0.100705177f, -0.055804994f, + 0.0274835322f, -0.00980951265f, -0.000163074146f, 0.00467492081f, + -0.00568934716f, 0.00484715449f, -0.00335428887f, 0.00194287405f, + -0.000931411923f, 0.000352004659f, -9.27265501e-05f, 1.14250561e-05f, + }, + { + -3.65884516e-05f, 6.37359626e-05f, 1.69760006e-05f, -0.000416312017f, + 0.00144883106f, -0.00346299331f, 0.00669340836f, -0.011051123f, + 0.0158964023f, -0.0198480636f, 0.0206481554f, -0.0149722425f, + -0.0023247432f, 0.0417539589f, -0.140543953f, 0.763356924f, + 0.46220997f, -0.183547407f, 0.10030362f, -0.0550526381f, + 0.0267064739f, -0.00915005151f, -0.000654205272f, 0.00500005623f, + -0.00587902358f, 0.00494177127f, -0.00339167449f, 0.00195171777f, + -0.000929669011f, 0.000348597416f, -9.06567584e-05f, 1.07816368e-05f, + }, + { + -3.4973571e-05f, 5.84274021e-05f, 2.95061564e-05f, -0.000439497118f, + 0.00148284878f, -0.00349962711f, 0.00670967111f, -0.0110032978f, + 0.0157167092f, -0.0194458999f, 0.0199165307f, -0.0137976436f, + -0.00405401643f, 0.0441468917f, -0.143642187f, 0.757394195f, + 0.470767885f, -0.184331506f, 0.0998493433f, -0.0542682484f, + 0.0259115566f, -0.00848214049f, -0.00114802097f, 0.00532480748f, + -0.00606707903f, 0.00503461622f, -0.00342764659f, 0.00195963401f, + -0.000927411718f, 0.000344950153f, -8.84950641e-05f, 1.01106862e-05f, + }, + { + -3.33734642e-05f, 5.31723308e-05f, 4.18580275e-05f, -0.000462200318f, + 0.00151578442f, -0.00353417965f, 0.00672243536f, -0.0109502962f, + 0.0155303273f, -0.0190364942f, 0.0191793088f, -0.0126231657f, + -0.0057701976f, 0.0464993082f, -0.146632105f, 0.751325846f, + 0.479302019f, -0.185032189f, 0.0993421748f, -0.0534520075f, + 0.0250990707f, -0.00780606316f, -0.0016442897f, 0.00564900693f, + -0.00625340641f, 0.00512562785f, -0.00346217398f, 0.00196660962f, + -0.000924635504f, 0.000341061736f, -8.62412271e-05f, 9.41208054e-06f, + }, + { + -3.17889499e-05f, 4.79730297e-05f, 5.40270412e-05f, -0.000484414632f, + 0.00154763111f, -0.00356665044f, 0.00673172064f, -0.0108921789f, + 0.0153373946f, -0.0186200943f, 0.0184368864f, -0.0114493882f, + -0.00747253513f, 0.048810415f, -0.149513751f, 0.745153546f, + 0.487810045f, -0.185648248f, 0.0987819806f, -0.0526041202f, + 0.0242693145f, -0.00712210825f, -0.00214277697f, 0.00597248739f, + -0.00643789815f, 0.00521474564f, -0.00349522661f, 0.00197263132f, + -0.000921335944f, 0.000336931058f, -8.389508e-05f, 8.6857126e-06f, + }, + { + -3.02208337e-05f, 4.28317071e-05f, 6.60087826e-05f, -0.00050613354f, + 0.00157838222f, -0.00359703996f, 0.00673754839f, -0.010829011f, + 0.0151380487f, -0.0181969516f, 0.0176896658f, -0.010276882f, + -0.00916028861f, 0.0510794297f, -0.152287185f, 0.738879204f, + 0.496289581f, -0.186178535f, 0.0981686339f, -0.0517248064f, + 0.0234226007f, -0.00643057236f, -0.00264324574f, 0.00629508123f, + -0.00662044808f, 0.00530190859f, -0.00352677447f, 0.00197768677f, + -0.000917508965f, 0.000332557218f, -8.14565065e-05f, 7.93149957e-06f, + }, + { + -2.86698851e-05f, 3.77505094e-05f, 7.77989917e-05f, -0.000527350698f, + 0.00160803227f, -0.00362535007f, 0.00673994143f, -0.0107608596f, + 0.0149324303f, -0.0177673176f, 0.01693804f, -0.00910621975f, + -0.0108327297f, 0.053305611f, -0.154952526f, 0.732504547f, + 0.504738271f, -0.186621904f, 0.0975020528f, -0.0508143008f, + 0.0225592479f, -0.0057317554f, -0.00314545725f, 0.00661661895f, + -0.00680094957f, 0.00538705708f, -0.00355678773f, 0.0019817641f, + -0.000913150609f, 0.00032793943f, -7.89254409e-05f, 7.14937505e-06f, + }, + { + -2.71368535e-05f, 3.27315211e-05f, 8.93935867e-05f, -0.000548060285f, + 0.00163657579f, -0.00365158333f, 0.00673892396f, -0.0106877917f, + 0.0147206821f, -0.0173314437f, 0.0161824115f, -0.00793796405f, + -0.012489141f, 0.0554882251f, -0.157509923f, 0.726031542f, + 0.513153791f, -0.186977252f, 0.0967821479f, -0.0498728603f, + 0.0216795895f, -0.00502596423f, -0.00364916935f, 0.00693693198f, + -0.00697929552f, 0.00547013199f, -0.00358523824f, 0.0019848512f, + -0.000908257323f, 0.000323076994f, -7.63018834e-05f, 6.33929676e-06f, + }, + { + -2.56224539e-05f, 2.77767467e-05f, 0.000100788631f, -0.000568256655f, + 0.00166400848f, -0.00367574347f, 0.00673452253f, -0.0106098782f, + 0.0145029482f, -0.0168895852f, 0.0154231768f, -0.00677267835f, + -0.0141288182f, 0.057626579f, -0.15995954f, 0.719462037f, + 0.521533847f, -0.187243462f, 0.0960088968f, -0.0489007458f, + 0.0207839627f, -0.00431350945f, -0.00415413827f, 0.00725584943f, + -0.00715538068f, 0.00555107417f, -0.00361209712f, 0.0019869369f, + -0.000902825675f, 0.000317969359f, -7.35858848e-05f, 5.5012415e-06f, + }, + { + -2.41273792e-05f, 2.28881254e-05f, 0.000111980356f, -0.00058793457f, + 0.00169032626f, -0.00369783561f, 0.00672676461f, -0.0105271926f, + 0.0142793776f, -0.016441999f, 0.0146607328f, -0.00561091769f, + -0.0157510713f, 0.0597199984f, -0.1623016f, 0.712797999f, + 0.529875994f, -0.187419474f, 0.0951822624f, -0.0478982553f, + 0.0198727194f, -0.00359470863f, -0.00466011837f, 0.00757320225f, + -0.00732909935f, 0.00562982587f, -0.00363733713f, 0.00198801048f, + -0.000896852522f, 0.000312616176f, -7.07775762e-05f, 4.63520973e-06f, + }, + { + -2.26522916e-05f, 1.80675252e-05f, 0.000122965153f, -0.00060708914f, + 0.0017155261f, -0.00371786556f, 0.00671567954f, -0.0104398085f, + 0.0140501177f, -0.0159889422f, 0.0138954772f, -0.00445323391f, + -0.0173552167f, 0.0617678314f, -0.164536357f, 0.706041336f, + 0.538177907f, -0.187504262f, 0.0943022519f, -0.0468656793f, + 0.0189462192f, -0.00286988285f, -0.00516686123f, 0.00788881909f, + -0.00750034675f, 0.00570632936f, -0.00366093102f, 0.00198806124f, + -0.000890335068f, 0.000307017181f, -6.78771248e-05f, 3.74122237e-06f, + }, + { + -2.11978295e-05f, 1.3316745e-05f, 0.000133739595f, -0.000625715766f, + 0.00173960521f, -0.00373583985f, 0.00670129713f, -0.0103478031f, + 0.013815321f, -0.0155306719f, 0.0131278047f, -0.00330017274f, + -0.0189405922f, 0.0637694597f, -0.166664109f, 0.699194014f, + 0.546437323f, -0.187496796f, 0.0933688954f, -0.0458033457f, + 0.0180048347f, -0.00213935901f, -0.00567411771f, 0.00820252951f, + -0.00766901812f, 0.00578052877f, -0.00368285249f, 0.00198707962f, + -0.000883270695f, 0.0003011722f, -6.48847781e-05f, 2.81932489e-06f, + }, + { + -1.97646023e-05f, 8.63751575e-06f, 0.000144300357f, -0.000643810316f, + 0.00176256162f, -0.00375176687f, 0.00668365043f, -0.0102512529f, + 0.0135751404f, -0.0150674507f, 0.0123581104f, -0.00215227529f, + -0.0205065403f, 0.0657242835f, -0.168685183f, 0.69225806f, + 0.554651916f, -0.187396079f, 0.0923822522f, -0.0447115861f, + 0.0170489419f, -0.00140346796f, -0.00618163543f, 0.00851416215f, + -0.00783500914f, 0.00585236726f, -0.00370307569f, 0.00198505563f, + -0.000875657075f, 0.000295081292f, -6.18008271e-05f, 1.8695838e-06f, + }, + { + -1.83531902e-05f, 4.03149534e-06f, 0.000154644367f, -0.000661368773f, + 0.0017843938f, -0.00376565545f, 0.0066627725f, -0.0101502389f, + 0.0133297294f, -0.0145995375f, 0.0115867881f, -0.00101007626f, + -0.0220524259f, 0.0676317289f, -0.170599923f, 0.685235441f, + 0.562819302f, -0.187201172f, 0.0913423896f, -0.0435907543f, + 0.0160789303f, -0.000662545965f, -0.00668916153f, 0.00882354751f, + -0.00799821783f, 0.00592179084f, -0.00372157572f, 0.00198198063f, + -0.000867492286f, 0.000288744603f, -5.86256538e-05f, 8.920901e-07f, + }, + { + -1.6964148e-05f, -4.99727662e-07f, 0.000164768629f, -0.000678387587f, + 0.00180510117f, -0.00377751514f, 0.00663869828f, -0.010044843f, + 0.0130792446f, -0.014127194f, 0.010814229f, 0.000125895196f, + -0.0235776193f, 0.0694912523f, -0.17240873f, 0.678128302f, + 0.570937216f, -0.186911136f, 0.0902494267f, -0.0424412154f, + 0.0150951995f, 8.30668287e-05f, -0.00719644083f, 0.0091305133f, + -0.00815854035f, 0.00598874455f, -0.00373832788f, 0.00197784579f, + -0.000858774467f, 0.00028216248f, -5.53596765e-05f, -1.13042134e-07f, + }, + { + -1.55980051e-05f, -4.95463428e-06f, 0.000174670378f, -0.000694863615f, + 0.00182468351f, -0.00378735736f, 0.00661146455f, -0.00993514806f, + 0.0128238453f, -0.0136506837f, 0.0100408243f, 0.00125511608f, + -0.0250815097f, 0.071302332f, -0.174112037f, 0.670938611f, + 0.579003394f, -0.186525062f, 0.0891034752f, -0.0412633531f, + 0.0140981553f, 0.000833025668f, -0.00770321768f, 0.00943489f, + -0.00831587613f, 0.00605317531f, -0.00375330891f, 0.00197264319f, + -0.000849502161f, 0.000275335304f, -5.20033936e-05f, -1.14567467e-06f, + }, + { + -1.42552608e-05f, -9.33177671e-06f, 0.000184346965f, -0.000710793887f, + 0.00184314104f, -0.00379519351f, 0.00658110855f, -0.00982123986f, + 0.0125636896f, -0.0131702675f, 0.0092669623f, 0.00237706956f, + -0.0265634973f, 0.0730644837f, -0.175710306f, 0.663668513f, + 0.58701545f, -0.18604207f, 0.087904714f, -0.0400575735f, + 0.0130882142f, 0.0015869817f, -0.0082092341f, 0.00973650534f, + -0.0084701227f, 0.0061150305f, -0.0037664955f, 0.00196636515f, + -0.000839674263f, 0.000268263742f, -4.85573619e-05f, -2.20564493e-06f, + }, + { + -1.2936388e-05f, -1.36297758e-05f, 0.000193795902f, -0.000726175786f, + 0.00186047493f, -0.00380103663f, 0.00654766941f, -0.0097032031f, + 0.0122989388f, -0.0126862098f, 0.00849302951f, 0.00349124591f, + -0.0280230008f, 0.0747772381f, -0.177204043f, 0.656320035f, + 0.59497118f, -0.185461313f, 0.0866532996f, -0.0388242863f, + 0.0120657999f, 0.002344582f, -0.00871423259f, 0.0100351917f, + -0.00862118136f, 0.00617425889f, -0.00377786579f, 0.00195900444f, + -0.000829289842f, 0.000260948553f, -4.50222105e-05f, -3.292766e-06f, + }, + { + -1.16418332e-05f, -1.78473238e-05f, 0.000203014904f, -0.000741007214f, + 0.00187668682f, -0.00380490022f, 0.00651118765f, -0.0095811272f, + 0.0120297549f, -0.0121987741f, 0.00771940965f, 0.00459714187f, + -0.0294594467f, 0.0764401406f, -0.17859377f, 0.648895383f, + 0.602868319f, -0.184781969f, 0.0853494555f, -0.0375639275f, + 0.0110313483f, 0.00310546998f, -0.00921795145f, 0.0103307767f, + -0.00876895245f, 0.00623080926f, -0.00378739834f, 0.0019505549f, + -0.000818348373f, 0.000253390579f, -4.13986272e-05f, -4.40682606e-06f, + }, + { + -1.03720167e-05f, -2.19831818e-05f, 0.000212001803f, -0.000755286077f, + 0.00189177867f, -0.00380679942f, 0.00647170423f, -0.00945510156f, + 0.0117563009f, -0.0117082223f, 0.00694648549f, 0.00569426082f, + -0.0308722835f, 0.078052789f, -0.179880053f, 0.641396642f, + 0.610704601f, -0.184003264f, 0.0839934126f, -0.036276944f, + 0.00998529885f, 0.00386928511f, -0.00972013175f, 0.0106230937f, + -0.00891333818f, 0.00628463319f, -0.0037950722f, 0.00194101024f, + -0.000806849566f, 0.000245590898f, -3.76873722e-05f, -5.54758753e-06f, + }, + { + -9.12733049e-06f, -2.60361812e-05f, 0.000220754606f, -0.00076901092f, + 0.00190575339f, -0.00380674913f, 0.00642926292f, -0.00932521746f, + 0.0114787426f, -0.0112148197f, 0.00617463561f, 0.00678211404f, + -0.0322609693f, 0.0796147883f, -0.181063503f, 0.633826017f, + 0.618477821f, -0.183124423f, 0.0825854242f, -0.0349638015f, + 0.00892810337f, 0.00463566324f, -0.0102205109f, 0.0109119713f, + -0.0090542417f, 0.0063356813f, -0.00380086782f, 0.00193036499f, + -0.000794793421f, 0.000237550761f, -3.38892714e-05f, -6.71478892e-06f, + }, + { + -7.90814101e-06f, -3.00052216e-05f, 0.000229271478f, -0.000782180345f, + 0.00191861414f, -0.00380476634f, 0.00638390658f, -0.00919156615f, + 0.0111972438f, -0.0107188281f, 0.00540423673f, 0.00786021911f, + -0.0336249769f, 0.0811257809f, -0.182144746f, 0.626185715f, + 0.626185715f, -0.182144746f, 0.0811257809f, -0.0336249769f, + 0.00786021911f, 0.00540423673f, -0.0107188281f, 0.0111972438f, + -0.00919156615f, 0.00638390658f, -0.00380476634f, 0.00191861414f, + -0.000782180345f, 0.000229271478f, -3.00052216e-05f, -7.90814101e-06f, + }, + { + -6.71478892e-06f, -3.38892714e-05f, 0.000237550761f, -0.000794793421f, + 0.00193036499f, -0.00380086782f, 0.0063356813f, -0.0090542417f, + 0.0109119713f, -0.0102205109f, 0.00463566324f, 0.00892810337f, + -0.0349638015f, 0.0825854242f, -0.183124423f, 0.618477821f, + 0.633826017f, -0.181063503f, 0.0796147883f, -0.0322609693f, + 0.00678211404f, 0.00617463561f, -0.0112148197f, 0.0114787426f, + -0.00932521746f, 0.00642926292f, -0.00380674913f, 0.00190575339f, + -0.00076901092f, 0.000220754606f, -2.60361812e-05f, -9.12733049e-06f, + }, + { + -5.54758753e-06f, -3.76873722e-05f, 0.000245590898f, -0.000806849566f, + 0.00194101024f, -0.0037950722f, 0.00628463319f, -0.00891333818f, + 0.0106230937f, -0.00972013175f, 0.00386928511f, 0.00998529885f, + -0.036276944f, 0.0839934126f, -0.184003264f, 0.610704601f, + 0.641396642f, -0.179880053f, 0.078052789f, -0.0308722835f, + 0.00569426082f, 0.00694648549f, -0.0117082223f, 0.0117563009f, + -0.00945510156f, 0.00647170423f, -0.00380679942f, 0.00189177867f, + -0.000755286077f, 0.000212001803f, -2.19831818e-05f, -1.03720167e-05f, + }, + { + -4.40682606e-06f, -4.13986272e-05f, 0.000253390579f, -0.000818348373f, + 0.0019505549f, -0.00378739834f, 0.00623080926f, -0.00876895245f, + 0.0103307767f, -0.00921795145f, 0.00310546998f, 0.0110313483f, + -0.0375639275f, 0.0853494555f, -0.184781969f, 0.602868319f, + 0.648895383f, -0.17859377f, 0.0764401406f, -0.0294594467f, + 0.00459714187f, 0.00771940965f, -0.0121987741f, 0.0120297549f, + -0.0095811272f, 0.00651118765f, -0.00380490022f, 0.00187668682f, + -0.000741007214f, 0.000203014904f, -1.78473238e-05f, -1.16418332e-05f, + }, + { + -3.292766e-06f, -4.50222105e-05f, 0.000260948553f, -0.000829289842f, + 0.00195900444f, -0.00377786579f, 0.00617425889f, -0.00862118136f, + 0.0100351917f, -0.00871423259f, 0.002344582f, 0.0120657999f, + -0.0388242863f, 0.0866532996f, -0.185461313f, 0.59497118f, + 0.656320035f, -0.177204043f, 0.0747772381f, -0.0280230008f, + 0.00349124591f, 0.00849302951f, -0.0126862098f, 0.0122989388f, + -0.0097032031f, 0.00654766941f, -0.00380103663f, 0.00186047493f, + -0.000726175786f, 0.000193795902f, -1.36297758e-05f, -1.2936388e-05f, + }, + { + -2.20564493e-06f, -4.85573619e-05f, 0.000268263742f, -0.000839674263f, + 0.00196636515f, -0.0037664955f, 0.0061150305f, -0.0084701227f, + 0.00973650534f, -0.0082092341f, 0.0015869817f, 0.0130882142f, + -0.0400575735f, 0.087904714f, -0.18604207f, 0.58701545f, + 0.663668513f, -0.175710306f, 0.0730644837f, -0.0265634973f, + 0.00237706956f, 0.0092669623f, -0.0131702675f, 0.0125636896f, + -0.00982123986f, 0.00658110855f, -0.00379519351f, 0.00184314104f, + -0.000710793887f, 0.000184346965f, -9.33177671e-06f, -1.42552608e-05f, + }, + { + -1.14567467e-06f, -5.20033936e-05f, 0.000275335304f, -0.000849502161f, + 0.00197264319f, -0.00375330891f, 0.00605317531f, -0.00831587613f, + 0.00943489f, -0.00770321768f, 0.000833025668f, 0.0140981553f, + -0.0412633531f, 0.0891034752f, -0.186525062f, 0.579003394f, + 0.670938611f, -0.174112037f, 0.071302332f, -0.0250815097f, + 0.00125511608f, 0.0100408243f, -0.0136506837f, 0.0128238453f, + -0.00993514806f, 0.00661146455f, -0.00378735736f, 0.00182468351f, + -0.000694863615f, 0.000174670378f, -4.95463428e-06f, -1.55980051e-05f, + }, + { + -1.13042134e-07f, -5.53596765e-05f, 0.00028216248f, -0.000858774467f, + 0.00197784579f, -0.00373832788f, 0.00598874455f, -0.00815854035f, + 0.0091305133f, -0.00719644083f, 8.30668287e-05f, 0.0150951995f, + -0.0424412154f, 0.0902494267f, -0.186911136f, 0.570937216f, + 0.678128302f, -0.17240873f, 0.0694912523f, -0.0235776193f, + 0.000125895196f, 0.010814229f, -0.014127194f, 0.0130792446f, + -0.010044843f, 0.00663869828f, -0.00377751514f, 0.00180510117f, + -0.000678387587f, 0.000164768629f, -4.99727662e-07f, -1.6964148e-05f, + }, + { + 8.920901e-07f, -5.86256538e-05f, 0.000288744603f, -0.000867492286f, + 0.00198198063f, -0.00372157572f, 0.00592179084f, -0.00799821783f, + 0.00882354751f, -0.00668916153f, -0.000662545965f, 0.0160789303f, + -0.0435907543f, 0.0913423896f, -0.187201172f, 0.562819302f, + 0.685235441f, -0.170599923f, 0.0676317289f, -0.0220524259f, + -0.00101007626f, 0.0115867881f, -0.0145995375f, 0.0133297294f, + -0.0101502389f, 0.0066627725f, -0.00376565545f, 0.0017843938f, + -0.000661368773f, 0.000154644367f, 4.03149534e-06f, -1.83531902e-05f, + }, + { + 1.8695838e-06f, -6.18008271e-05f, 0.000295081292f, -0.000875657075f, + 0.00198505563f, -0.00370307569f, 0.00585236726f, -0.00783500914f, + 0.00851416215f, -0.00618163543f, -0.00140346796f, 0.0170489419f, + -0.0447115861f, 0.0923822522f, -0.187396079f, 0.554651916f, + 0.69225806f, -0.168685183f, 0.0657242835f, -0.0205065403f, + -0.00215227529f, 0.0123581104f, -0.0150674507f, 0.0135751404f, + -0.0102512529f, 0.00668365043f, -0.00375176687f, 0.00176256162f, + -0.000643810316f, 0.000144300357f, 8.63751575e-06f, -1.97646023e-05f, + }, + { + 2.81932489e-06f, -6.48847781e-05f, 0.0003011722f, -0.000883270695f, + 0.00198707962f, -0.00368285249f, 0.00578052877f, -0.00766901812f, + 0.00820252951f, -0.00567411771f, -0.00213935901f, 0.0180048347f, + -0.0458033457f, 0.0933688954f, -0.187496796f, 0.546437323f, + 0.699194014f, -0.166664109f, 0.0637694597f, -0.0189405922f, + -0.00330017274f, 0.0131278047f, -0.0155306719f, 0.013815321f, + -0.0103478031f, 0.00670129713f, -0.00373583985f, 0.00173960521f, + -0.000625715766f, 0.000133739595f, 1.3316745e-05f, -2.11978295e-05f, + }, + { + 3.74122237e-06f, -6.78771248e-05f, 0.000307017181f, -0.000890335068f, + 0.00198806124f, -0.00366093102f, 0.00570632936f, -0.00750034675f, + 0.00788881909f, -0.00516686123f, -0.00286988285f, 0.0189462192f, + -0.0468656793f, 0.0943022519f, -0.187504262f, 0.538177907f, + 0.706041336f, -0.164536357f, 0.0617678314f, -0.0173552167f, + -0.00445323391f, 0.0138954772f, -0.0159889422f, 0.0140501177f, + -0.0104398085f, 0.00671567954f, -0.00371786556f, 0.0017155261f, + -0.00060708914f, 0.000122965153f, 1.80675252e-05f, -2.26522916e-05f, + }, + { + 4.63520973e-06f, -7.07775762e-05f, 0.000312616176f, -0.000896852522f, + 0.00198801048f, -0.00363733713f, 0.00562982587f, -0.00732909935f, + 0.00757320225f, -0.00466011837f, -0.00359470863f, 0.0198727194f, + -0.0478982553f, 0.0951822624f, -0.187419474f, 0.529875994f, + 0.712797999f, -0.1623016f, 0.0597199984f, -0.0157510713f, + -0.00561091769f, 0.0146607328f, -0.016441999f, 0.0142793776f, + -0.0105271926f, 0.00672676461f, -0.00369783561f, 0.00169032626f, + -0.00058793457f, 0.000111980356f, 2.28881254e-05f, -2.41273792e-05f, + }, + { + 5.5012415e-06f, -7.35858848e-05f, 0.000317969359f, -0.000902825675f, + 0.0019869369f, -0.00361209712f, 0.00555107417f, -0.00715538068f, + 0.00725584943f, -0.00415413827f, -0.00431350945f, 0.0207839627f, + -0.0489007458f, 0.0960088968f, -0.187243462f, 0.521533847f, + 0.719462037f, -0.15995954f, 0.057626579f, -0.0141288182f, + -0.00677267835f, 0.0154231768f, -0.0168895852f, 0.0145029482f, + -0.0106098782f, 0.00673452253f, -0.00367574347f, 0.00166400848f, + -0.000568256655f, 0.000100788631f, 2.77767467e-05f, -2.56224539e-05f, + }, + { + 6.33929676e-06f, -7.63018834e-05f, 0.000323076994f, -0.000908257323f, + 0.0019848512f, -0.00358523824f, 0.00547013199f, -0.00697929552f, + 0.00693693198f, -0.00364916935f, -0.00502596423f, 0.0216795895f, + -0.0498728603f, 0.0967821479f, -0.186977252f, 0.513153791f, + 0.726031542f, -0.157509923f, 0.0554882251f, -0.012489141f, + -0.00793796405f, 0.0161824115f, -0.0173314437f, 0.0147206821f, + -0.0106877917f, 0.00673892396f, -0.00365158333f, 0.00163657579f, + -0.000548060285f, 8.93935867e-05f, 3.27315211e-05f, -2.71368535e-05f, + }, + { + 7.14937505e-06f, -7.89254409e-05f, 0.00032793943f, -0.000913150609f, + 0.0019817641f, -0.00355678773f, 0.00538705708f, -0.00680094957f, + 0.00661661895f, -0.00314545725f, -0.0057317554f, 0.0225592479f, + -0.0508143008f, 0.0975020528f, -0.186621904f, 0.504738271f, + 0.732504547f, -0.154952526f, 0.053305611f, -0.0108327297f, + -0.00910621975f, 0.01693804f, -0.0177673176f, 0.0149324303f, + -0.0107608596f, 0.00673994143f, -0.00362535007f, 0.00160803227f, + -0.000527350698f, 7.77989917e-05f, 3.77505094e-05f, -2.86698851e-05f, + }, + { + 7.93149957e-06f, -8.14565065e-05f, 0.000332557218f, -0.000917508965f, + 0.00197768677f, -0.00352677447f, 0.00530190859f, -0.00662044808f, + 0.00629508123f, -0.00264324574f, -0.00643057236f, 0.0234226007f, + -0.0517248064f, 0.0981686339f, -0.186178535f, 0.496289581f, + 0.738879204f, -0.152287185f, 0.0510794297f, -0.00916028861f, + -0.010276882f, 0.0176896658f, -0.0181969516f, 0.0151380487f, + -0.010829011f, 0.00673754839f, -0.00359703996f, 0.00157838222f, + -0.00050613354f, 6.60087826e-05f, 4.28317071e-05f, -3.02208337e-05f, + }, + { + 8.6857126e-06f, -8.389508e-05f, 0.000336931058f, -0.000921335944f, + 0.00197263132f, -0.00349522661f, 0.00521474564f, -0.00643789815f, + 0.00597248739f, -0.00214277697f, -0.00712210825f, 0.0242693145f, + -0.0526041202f, 0.0987819806f, -0.185648248f, 0.487810045f, + 0.745153546f, -0.149513751f, 0.048810415f, -0.00747253513f, + -0.0114493882f, 0.0184368864f, -0.0186200943f, 0.0153373946f, + -0.0108921789f, 0.00673172064f, -0.00356665044f, 0.00154763111f, + -0.000484414632f, 5.40270412e-05f, 4.79730297e-05f, -3.17889499e-05f, + }, + { + 9.41208054e-06f, -8.62412271e-05f, 0.000341061736f, -0.000924635504f, + 0.00196660962f, -0.00346217398f, 0.00512562785f, -0.00625340641f, + 0.00564900693f, -0.0016442897f, -0.00780606316f, 0.0250990707f, + -0.0534520075f, 0.0993421748f, -0.185032189f, 0.479302019f, + 0.751325846f, -0.146632105f, 0.0464993082f, -0.0057701976f, + -0.0126231657f, 0.0191793088f, -0.0190364942f, 0.0155303273f, + -0.0109502962f, 0.00672243536f, -0.00353417965f, 0.00151578442f, + -0.000462200318f, 4.18580275e-05f, 5.31723308e-05f, -3.33734642e-05f, + }, + { + 1.01106862e-05f, -8.84950641e-05f, 0.000344950153f, -0.000927411718f, + 0.00195963401f, -0.00342764659f, 0.00503461622f, -0.00606707903f, + 0.00532480748f, -0.00114802097f, -0.00848214049f, 0.0259115566f, + -0.0542682484f, 0.0998493433f, -0.184331506f, 0.470767885f, + 0.757394195f, -0.143642187f, 0.0441468917f, -0.00405401643f, + -0.0137976436f, 0.0199165307f, -0.0194458999f, 0.0157167092f, + -0.0110032978f, 0.00670967111f, -0.00349962711f, 0.00148284878f, + -0.000439497118f, 2.95061564e-05f, 5.84274021e-05f, -3.4973571e-05f, + }, + { + 1.07816368e-05f, -9.06567584e-05f, 0.000348597416f, -0.000929669011f, + 0.00195171777f, -0.00339167449f, 0.00494177127f, -0.00587902358f, + 0.00500005623f, -0.000654205272f, -0.00915005151f, 0.0267064739f, + -0.0550526381f, 0.10030362f, -0.183547407f, 0.46220997f, + 0.763356924f, -0.140543953f, 0.0417539589f, -0.0023247432f, + -0.0149722425f, 0.0206481554f, -0.0198480636f, 0.0158964023f, + -0.011051123f, 0.00669340836f, -0.00346299331f, 0.00144883106f, + -0.000416312017f, 1.69760006e-05f, 6.37359626e-05f, -3.65884516e-05f, + }, + { + 1.14250561e-05f, -9.27265501e-05f, 0.000352004659f, -0.000931411923f, + 0.00194287405f, -0.00335428887f, 0.00484715449f, -0.00568934716f, + 0.00467492081f, -0.000163074146f, -0.00980951265f, 0.0274835322f, + -0.055804994f, 0.100705177f, -0.182681099f, 0.453630656f, + 0.769212127f, -0.137337416f, 0.0393213332f, -0.000583140645f, + -0.0161463823f, 0.021373786f, -0.0202427432f, 0.0160692763f, + -0.0110937115f, 0.00667362893f, -0.00342427962f, 0.00141373917f, + -0.000392652291f, 4.27229452e-06f, 6.90956513e-05f, -3.82172475e-05f, + }, + { + 1.20410887e-05f, -9.47047156e-05f, 0.000355173252f, -0.000932645402f, + 0.00193311635f, -0.00331552117f, 0.00475082779f, -0.00549815688f, + 0.00434956606f, 0.000325143425f, -0.0104602445f, 0.0282424539f, + -0.0565251447f, 0.101054192f, -0.181733787f, 0.445032269f, + 0.774958193f, -0.134022623f, 0.0368498713f, 0.00117001776f, + -0.0173194781f, 0.022093026f, -0.0206296947f, 0.016235197f, + -0.0111310054f, 0.00665031699f, -0.00338348909f, 0.00137758139f, + -0.000368525536f, -8.60007094e-06f, 7.45040525e-05f, -3.98590819e-05f, + }, + { + 1.26298974e-05f, -9.6591597e-05f, 0.00035810456f, -0.000933374511f, + 0.00192245888f, -0.00327540329f, 0.00465285406f, -0.00530556031f, + 0.00402415777f, 0.000810221303f, -0.011101976f, 0.0289829671f, + -0.0572129339f, 0.101350896f, -0.180706754f, 0.436417222f, + 0.780593336f, -0.130599633f, 0.0343404375f, 0.00293394784f, + -0.0184909441f, 0.0228054784f, -0.0210086759f, 0.0163940396f, + -0.0111629497f, 0.00662345719f, -0.00334062544f, 0.00134036678f, + -0.000343939726f, -2.16360477e-05f, 7.99586851e-05f, -4.15130453e-05f, + }, + { + 1.31916649e-05f, -9.83875871e-05f, 0.000360800157f, -0.000933604431f, + 0.00191091583f, -0.00323396735f, 0.00455329521f, -0.00511166407f, + 0.0036988589f, 0.00129193638f, -0.01173444f, 0.0297048185f, + -0.0578682311f, 0.101595506f, -0.179601267f, 0.427787781f, + 0.786115944f, -0.127068579f, 0.0317939334f, 0.00470785564f, + -0.0196601916f, 0.0235107485f, -0.0213794503f, 0.0165456794f, + -0.0111894915f, 0.00659303693f, -0.00329569331f, 0.00130210491f, + -0.000318903098f, -3.48304311e-05f, 8.54570026e-05f, -4.31782064e-05f, + }, + { + 1.37265888e-05f, -0.000100093122f, 0.000363261701f, -0.000933340809f, + 0.00189850212f, -0.00319124642f, 0.00445221504f, -0.00491657667f, + 0.00337383244f, 0.00177006831f, -0.0123573784f, 0.0304077566f, + -0.0584909096f, 0.101788282f, -0.178418607f, 0.419146359f, + 0.79152441f, -0.123429619f, 0.0292112771f, 0.00649093743f, + -0.0208266266f, 0.0242084414f, -0.0217417832f, 0.0166899916f, + -0.0112105804f, 0.00655904552f, -0.00324869901f, 0.00126280601f, + -0.000293424295f, -4.81778588e-05f, 9.09963856e-05f, -4.48536084e-05f, + }, + { + 1.42348863e-05f, -0.000101708705f, 0.000365490967f, -0.000932589232f, + 0.00188523252f, -0.00314727332f, 0.00434967596f, -0.00472040381f, + 0.00304923998f, 0.00224440033f, -0.0129705351f, 0.0310915448f, + -0.059080869f, 0.101929516f, -0.177160099f, 0.410495281f, + 0.796817124f, -0.119682938f, 0.0265934132f, 0.00828237925f, + -0.0219896547f, 0.0248981696f, -0.02209544f, 0.0168268587f, + -0.0112261679f, 0.00652147317f, -0.00319964997f, 0.00122248125f, + -0.000267512281f, -6.16728139e-05f, 9.65741565e-05f, -4.65382655e-05f, + }, + { + 1.47167921e-05f, -0.000103234866f, 0.000367489847f, -0.000931355753f, + 0.00187112251f, -0.00310208183f, 0.00424574222f, -0.00452325307f, + 0.00272524171f, 0.00271471869f, -0.0135736652f, 0.0317559615f, + -0.0596380271f, 0.102019496f, -0.175827071f, 0.401836902f, + 0.801992416f, -0.11582876f, 0.0239413101f, 0.0100813583f, + -0.023148682f, 0.0255795382f, -0.0224401932f, 0.0169561617f, + -0.0112362094f, 0.00648031244f, -0.00314855389f, 0.00118114206f, + -0.000241176327f, -7.53096392e-05f, 0.000102187572f, -4.82311661e-05f, + }, + { + 1.51725544e-05f, -0.00010467221f, 0.00036926032f, -0.000929646427f, + 0.00185618771f, -0.00305570546f, 0.00414047716f, -0.00432523061f, + 0.00240199687f, 0.00318081258f, -0.0141665265f, 0.0324007869f, + -0.0601623058f, 0.102058552f, -0.174420878f, 0.393173516f, + 0.807048857f, -0.111867368f, 0.021255957f, 0.0118870437f, + -0.0243031103f, 0.0262521598f, -0.0227758158f, 0.0170777924f, + -0.0112406602f, 0.00643555773f, -0.00309542078f, 0.00113880076f, + -0.000214426065f, -8.90825177e-05f, 0.000107833817f, -4.9931281e-05f, + }, + { + 1.56024416e-05f, -0.000106021362f, 0.000370804511f, -0.000927467539f, + 0.00184044405f, -0.00300817867f, 0.0040339455f, -0.00412644353f, + 0.0020796624f, 0.003642475f, -0.0147488844f, 0.0330258198f, + -0.060653653f, 0.102047034f, -0.172942877f, 0.384507507f, + 0.811984956f, -0.107799068f, 0.018538367f, 0.0136985937f, + -0.02545234f, 0.026915649f, -0.023102086f, 0.0171916354f, + -0.0112394821f, 0.00638720486f, -0.0030402611f, 0.00109547027f, + -0.00018727142f, -0.000102985476f, 0.000113510039f, -5.16375476e-05f, + }, + { + 1.60067339e-05f, -0.000107283013f, 0.000372124632f, -0.000924825668f, + 0.00182390772f, -0.00295953569f, 0.00392621104f, -0.00392699707f, + 0.001758395f, 0.00409950269f, -0.0153205106f, 0.0336308628f, + -0.0611120313f, 0.101985291f, -0.171394452f, 0.375841141f, + 0.816799164f, -0.103624195f, 0.0157895759f, 0.0155151617f, + -0.0265957713f, 0.0275696218f, -0.0234187823f, 0.0172975883f, + -0.011232635f, 0.00633525196f, -0.00298308628f, 0.00105116412f, + -0.000159722666f, -0.000117012394f, 0.000119213299f, -5.33488856e-05f, + }, + { + 1.63857276e-05f, -0.000108457898f, 0.000373222952f, -0.000921727449f, + 0.00180659513f, -0.00290981168f, 0.00381733826f, -0.00372699695f, + 0.00143834867f, 0.00455169473f, -0.0158811864f, 0.0342157409f, + -0.0615374185f, 0.101873733f, -0.169777006f, 0.367176801f, + 0.821490109f, -0.0993431434f, 0.0130106425f, 0.0173358899f, + -0.0277328063f, 0.0282136947f, -0.0237256885f, 0.0173955448f, + -0.0112200854f, 0.00627969857f, -0.00292390981f, 0.00100589637f, + -0.00013179041f, -0.000131157023f, 0.000124940605f, -5.50641853e-05f, + }, + { + 1.67397338e-05f, -0.000109546781f, 0.000374101917f, -0.000918179809f, + 0.00178852293f, -0.00285904133f, 0.00370739121f, -0.00352654839f, + 0.00111967651f, 0.00499885576f, -0.0164306946f, 0.0347802751f, + -0.0619298108f, 0.101712748f, -0.168091953f, 0.358516723f, + 0.826056361f, -0.0949563384f, 0.0102026463f, 0.0191599168f, + -0.0288628414f, 0.0288474914f, -0.0240225941f, 0.0174854062f, + -0.0112017998f, 0.00622054609f, -0.0028627452f, 0.000959681987f, + -0.00010348556f, -0.000145412967f, 0.000130688917f, -5.67823117e-05f, + }, + { + 1.70690764e-05f, -0.000110550478f, 0.000374764f, -0.000914189732f, + 0.00176970835f, -0.00280726003f, 0.0035964346f, -0.00332575594f, + 0.000802529801f, 0.00544079253f, -0.0169688277f, 0.0353243127f, + -0.0622892156f, 0.101502784f, -0.166340724f, 0.349863261f, + 0.83049655f, -0.0904642269f, 0.00736668706f, 0.0209863689f, + -0.0299852788f, 0.0294706356f, -0.0243092868f, 0.0175670739f, + -0.0111777494f, 0.00615779823f, -0.00279960805f, 0.00091253646f, + -7.48193415e-05f, -0.000159773655f, 0.000136455114f, -5.85021153e-05f, + }, + { + 1.73740973e-05f, -0.000111469853f, 0.000375211792f, -0.000909764494f, + 0.00175016839f, -0.0027545034f, 0.00348453294f, -0.00312472368f, + 0.00048705787f, 0.00587731553f, -0.0174953863f, 0.0358476974f, + -0.0626156628f, 0.101244271f, -0.164524779f, 0.34121868f, + 0.834809422f, -0.085867323f, 0.00450388901f, 0.0228143707f, + -0.031099515f, 0.0300827548f, -0.0245855656f, 0.0176404584f, + -0.0111479061f, 0.00609146059f, -0.00273451465f, 0.000864475791f, + -4.58033137e-05f, -0.000174232424f, 0.000142236051f, -6.02224209e-05f, + }, + { + 1.76551457e-05f, -0.000112305803f, 0.000375447999f, -0.000904911314f, + 0.00172992051f, -0.00270080729f, 0.0033717507f, -0.00292355497f, + 0.000173408436f, 0.00630824082f, -0.0180101767f, 0.0363502987f, + -0.0629091859f, 0.100937687f, -0.162645563f, 0.332585275f, + 0.838993609f, -0.0811661631f, 0.00161539461f, 0.0246430375f, + -0.0322049521f, 0.0306834765f, -0.0248512272f, 0.0177054685f, + -0.0111122467f, 0.00602153968f, -0.00266748248f, 0.000815516803f, + -1.64493304e-05f, -0.000188782447f, 0.000148028514f, -6.19420243e-05f, + }, + { + 1.79125836e-05f, -0.000113059257f, 0.000375475356f, -0.000899637875f, + 0.00170898228f, -0.00264620734f, 0.0032581524f, -0.00272235274f, + -0.000138272706f, 0.00673338631f, -0.0185130127f, 0.0368319787f, + -0.0631698593f, 0.100583509f, -0.160704538f, 0.323965251f, + 0.843047917f, -0.0763613135f, -0.00129763153f, 0.0264714807f, + -0.0333009884f, 0.0312724374f, -0.025106078f, 0.0177620202f, + -0.0110707479f, 0.00594804529f, -0.00259853038f, 0.000765676901f, + 1.32304358e-05f, -0.000203416756f, 0.000153829198f, -6.36597033e-05f, + }, + { + 1.81467913e-05f, -0.000113731199f, 0.000375296775f, -0.000893951743f, + 0.00168737164f, -0.00259073987f, 0.00314380229f, -0.00252121896f, + -0.00044784174f, 0.00715257553f, -0.0190037116f, 0.0372926332f, + -0.0633977354f, 0.100182258f, -0.158703223f, 0.315360934f, + 0.846971154f, -0.0714533925f, -0.0042340057f, 0.0282988027f, + -0.0343870297f, 0.0318492763f, -0.0253499206f, 0.017810028f, + -0.0110233929f, 0.00587098813f, -0.00252767792f, 0.00071497407f, + 4.32235101e-05f, -0.000218128276f, 0.000159634801f, -6.53742172e-05f, + }, + { + 1.83581524e-05f, -0.00011432264f, 0.000374915136f, -0.000887860719f, + 0.00166510663f, -0.00253444095f, 0.00302876462f, -0.00232025515f, + -0.000755156798f, 0.00756563526f, -0.0194821022f, 0.0377321467f, + -0.0635929182f, 0.0997344479f, -0.156643078f, 0.306774527f, + 0.850762069f, -0.0664430633f, -0.00719252322f, 0.0301241018f, + -0.0354624726f, 0.0324136317f, -0.0255825706f, 0.0178494155f, + -0.010970165f, 0.00579038123f, -0.00245494628f, 0.000663426879f, + 7.35171197e-05f, -0.000232909777f, 0.000165441947f, -6.70843074e-05f, + }, + { + 1.85470653e-05f, -0.00011483462f, 0.000374333496f, -0.000881372776f, + 0.00164220529f, -0.00247734692f, 0.00291310367f, -0.00211956142f, + -0.00106007815f, 0.00797239598f, -0.0199480206f, 0.0381504335f, + -0.0637555122f, 0.0992406234f, -0.154525638f, 0.298208266f, + 0.854419589f, -0.0613310076f, -0.0101719601f, 0.0319464765f, + -0.0365267284f, 0.0329651535f, -0.0258038435f, 0.0178801101f, + -0.0109110521f, 0.00570623949f, -0.00238035759f, 0.000611054653f, + 0.000104098181f, -0.000247753895f, 0.0001712472f, -6.87886786e-05f, + }, + { + 1.87139394e-05f, -0.000115268231f, 0.000373554911f, -0.000874496007f, + 0.00161868602f, -0.00241949433f, 0.00279688323f, -0.00191923801f, + -0.00136246812f, 0.00837269332f, -0.0204013046f, 0.0385474041f, + -0.0638856217f, 0.0987013429f, -0.152352408f, 0.289664388f, + 0.857942581f, -0.0561179556f, -0.0131710758f, 0.0337650143f, + -0.0375792012f, 0.0335034914f, -0.0260135606f, 0.017902039f, + -0.0108460439f, 0.00561858015f, -0.00230393535f, 0.00055787724f, + 0.000134953341f, -0.000262653193f, 0.000177047084f, -7.04860431e-05f, + }, + { + 1.88591948e-05f, -0.000115624578f, 0.00037258261f, -0.000867238618f, + 0.00159456744f, -0.00236091972f, 0.00268016732f, -0.00171938329f, + -0.00166219112f, 0.0087663671f, -0.0208418071f, 0.0389229916f, + -0.0639833957f, 0.0981171876f, -0.150124922f, 0.281145036f, + 0.861329973f, -0.0508046746f, -0.0161886066f, 0.0355787985f, + -0.0386192985f, 0.0340282954f, -0.0262115467f, 0.0179151371f, + -0.010775133f, 0.00552742137f, -0.00222570379f, 0.000503915129f, + 0.000166068945f, -0.000277600047f, 0.000182838063f, -7.21750694e-05f, + }, + { + 1.89832554e-05f, -0.000115904812f, 0.000371419766f, -0.000859608874f, + 0.00156986807f, -0.00230165967f, 0.00256301905f, -0.00152009528f, + -0.0019591141f, 0.00915326085f, -0.021269381f, 0.0392771326f, + -0.0640489832f, 0.0974887609f, -0.147844717f, 0.272652417f, + 0.86458087f, -0.0453919694f, -0.0192232747f, 0.0373869166f, + -0.0396464355f, 0.0345392302f, -0.0263976324f, 0.017919343f, + -0.0106983157f, 0.00543278456f, -0.00214568852f, 0.000449189392f, + 0.000197431029f, -0.000292586745f, 0.000188616585f, -7.38544186e-05f, + }, + { + 1.90865612e-05f, -0.00011611011f, 0.000370069727f, -0.000851615157f, + 0.00154460664f, -0.00224175071f, 0.00244550197f, -0.00132147083f, + -0.00225310563f, 0.00953322276f, -0.0216838922f, 0.0396097749f, + -0.0640825331f, 0.096816659f, -0.145513341f, 0.264188707f, + 0.867694139f, -0.0398806818f, -0.0222737808f, 0.0391884446f, + -0.0406600274f, 0.0350359604f, -0.0265716556f, 0.0179145988f, + -0.010615591f, 0.00533469208f, -0.00206391653f, 0.000393721712f, + 0.000229025405f, -0.000307605485f, 0.000194379041f, -7.55227302e-05f, + }, + { + 1.91695544e-05f, -0.000116241674f, 0.000368535839f, -0.000843266142f, + 0.00151880225f, -0.00218122988f, 0.00232767873f, -0.00112360588f, + -0.00254403776f, 0.00990610477f, -0.0220852084f, 0.0399208851f, + -0.0640842393f, 0.096101515f, -0.143132344f, 0.255756021f, + 0.870669007f, -0.0342716947f, -0.0253388118f, 0.0409824587f, + -0.0416594855f, 0.0355181508f, -0.0267334562f, 0.0179008469f, + -0.0105269607f, 0.00523316814f, -0.00198041555f, 0.000337534351f, + 0.00026083758f, -0.000322648295f, 0.000200121751f, -7.71786363e-05f, + }, + { + 1.92326916e-05f, -0.000116300747f, 0.000366821536f, -0.000834570383f, + 0.00149247376f, -0.00212013349f, 0.0022096124f, -0.000926595123f, + -0.00283178361f, 0.0102717644f, -0.0224732105f, 0.0402104259f, + -0.0640542954f, 0.0953439847f, -0.140703291f, 0.247356534f, + 0.87350446f, -0.0285659209f, -0.0284170341f, 0.0427680351f, + -0.0426442362f, 0.035985481f, -0.0268828776f, 0.0178780425f, + -0.0104324305f, 0.00512823975f, -0.00189521466f, 0.000280650158f, + 0.000292852754f, -0.000337707141f, 0.000205841017f, -7.88207326e-05f, + }, + { + 1.92764292e-05f, -0.000116288582f, 0.00036493034f, -0.000825536554f, + 0.00146564026f, -0.0020584981f, 0.0020913647f, -0.000730532338f, + -0.00311621954f, 0.0106300628f, -0.022847781f, 0.040478386f, + -0.0639928952f, 0.0945447162f, -0.138227761f, 0.238992333f, + 0.876199663f, -0.0227643196f, -0.0315071009f, 0.0445442423f, + -0.0436137021f, 0.0364376307f, -0.0270197745f, 0.0178461354f, + -0.0103320088f, 0.00501993531f, -0.00180834404f, 0.000223092618f, + 0.000325055968f, -0.000352773844f, 0.000211533115f, -8.04476222e-05f, + }, + { + 1.93012347e-05f, -0.000116206465f, 0.000362865743f, -0.00081617356f, + 0.00143832085f, -0.00199636072f, 0.00197299803f, -0.000535510131f, + -0.00339722424f, 0.0109808668f, -0.0232088137f, 0.0407247581f, + -0.0639002696f, 0.0937043875f, -0.135707304f, 0.230665505f, + 0.878753841f, -0.0168678798f, -0.0346076451f, 0.0463101529f, + -0.0445673168f, 0.0368742831f, -0.0271440018f, 0.0178050902f, + -0.0102257049f, 0.00490828604f, -0.00171983475f, 0.0001648857f, + 0.000357431883f, -0.00036784017f, 0.000217194276f, -8.20578862e-05f, + }, + { + 1.9307583e-05f, -0.000116055715f, 0.000360631413f, -0.00080649025f, + 0.00141053495f, -0.0019337578f, 0.00185457338f, -0.000341619743f, + -0.00367467897f, 0.0113240471f, -0.0235562101f, 0.0409495458f, + -0.0637766495f, 0.0928236842f, -0.133143514f, 0.22237812f, + 0.88116622f, -0.0108776297f, -0.0377172902f, 0.0480648354f, + -0.0455045104f, 0.0372951329f, -0.0272554215f, 0.0177548639f, + -0.0101135345f, 0.00479332311f, -0.00162971939f, 0.00010605398f, + 0.000389964989f, -0.000382897735f, 0.000222820672f, -8.36500767e-05f, + }, + { + 1.92959506e-05f, -0.000115837654f, 0.000358230987f, -0.000796495646f, + 0.00138230179f, -0.00187072577f, 0.00173615187f, -0.000148951382f, + -0.00394846732f, 0.0116594769f, -0.0238898769f, 0.0411527604f, + -0.0636222884f, 0.0919033214f, -0.130537972f, 0.214132234f, + 0.883436024f, -0.00479463488f, -0.0408346355f, 0.0498073585f, + -0.0464247204f, 0.0376998782f, -0.0273539033f, 0.0176954288f, + -0.00999551453f, 0.00467508193f, -0.00153803127f, 4.66225902e-05f, + 0.000422639481f, -0.000397938129f, 0.000228408462f, -8.52227458e-05f, + }, + { + 1.92668249e-05f, -0.000115553652f, 0.000355668133f, -0.000786198769f, + 0.00135364069f, -0.00180730107f, 0.00161779416f, 4.24060891e-05f, + -0.00421847589f, 0.0119870389f, -0.0242097266f, 0.0413344279f, + -0.0634374395f, 0.090944007f, -0.127892286f, 0.205929875f, + 0.885562658f, 0.00138000422f, -0.0439582765f, 0.051536791f, + -0.0473273918f, 0.0380882211f, -0.0274393167f, 0.0176267549f, + -0.00987166539f, 0.00455359789f, -0.00144480506f, -1.33828025e-05f, + 0.000455439382f, -0.000412952766f, 0.000233953775f, -8.67744384e-05f, + }, + { + 1.92206935e-05f, -0.000115205075f, 0.000352946634f, -0.0007756087f, + 0.0013245712f, -0.00174351991f, 0.00149956008f, 0.00023236507f, + -0.0044845934f, 0.0123066157f, -0.0245156847f, 0.0414945856f, + -0.0632223859f, 0.0899464786f, -0.12520802f, 0.197773054f, + 0.887545407f, 0.00764515065f, -0.0470867865f, 0.0532522015f, + -0.0482119694f, 0.0384598747f, -0.0275115464f, 0.01754882f, + -0.00974201038f, 0.00442891056f, -0.0013500764f, -7.39359893e-05f, + 0.000488348422f, -0.000427933031f, 0.000239452711f, -8.83036773e-05f, + }, + { + 1.91580493e-05f, -0.000114793322f, 0.000350070244f, -0.000764734636f, + 0.00129511254f, -0.00167941849f, 0.00138150924f, 0.000420839351f, + -0.00474671181f, 0.0126180965f, -0.0248076785f, 0.0416332819f, + -0.0629774183f, 0.0889114738f, -0.122486785f, 0.189663753f, + 0.889383733f, 0.0139996316f, -0.0502187237f, 0.054952655f, + -0.0490779169f, 0.0388145596f, -0.027570473f, 0.0174616054f, + -0.00960657652f, 0.00430105906f, -0.00125388207f, -0.000135010225f, + 0.000521350128f, -0.000442870194f, 0.000244901312f, -8.98089784e-05f, + }, + { + 1.90793908e-05f, -0.00011431981f, 0.000347042835f, -0.000753585773f, + 0.00126528414f, -0.00161503302f, 0.00126370054f, 0.000607743976f, + -0.00500472542f, 0.0129213771f, -0.0250856448f, 0.0417505652f, + -0.0627028272f, 0.0878397599f, -0.119730175f, 0.181603938f, + 0.891077042f, 0.0204422399f, -0.0533526428f, 0.0566372275f, + -0.0499246828f, 0.0391519926f, -0.0276159905f, 0.0173650961f, + -0.00946539361f, 0.00417008577f, -0.00115625968f, -0.000196578272f, + 0.000554427796f, -0.000457755465f, 0.000250295649f, -9.12888427e-05f, + }, + { + 1.89852162e-05f, -0.000113785965f, 0.000343868247f, -0.000742171484f, + 0.00123510568f, -0.00155039912f, 0.00114619231f, 0.000792995619f, + -0.00525853131f, 0.013216354f, -0.0253495276f, 0.04184651f, + -0.0623989403f, 0.086732097f, -0.116939791f, 0.173595518f, + 0.892624855f, 0.0269717295f, -0.0564870797f, 0.0583049878f, + -0.050751742f, 0.0394719094f, -0.0276479963f, 0.0172592811f, + -0.00931849517f, 0.00403603492f, -0.00105724821f, -0.000258612359f, + 0.000587564544f, -0.000472579966f, 0.000255631719f, -9.27417641e-05f, + }, + { + 1.88760332e-05f, -0.000113193237f, 0.00034055041f, -0.000730500906f, + 0.00120459648f, -0.00148555252f, 0.00102904218f, 0.000976512092f, + -0.00550802983f, 0.0135029331f, -0.0255992822f, 0.0419211872f, + -0.0620660707f, 0.0855892673f, -0.114117227f, 0.165640444f, + 0.894026697f, 0.0335868262f, -0.0596205555f, 0.0599550121f, + -0.0515585616f, 0.0397740491f, -0.0276663918f, 0.0171441585f, + -0.00916591566f, 0.00389895239f, -0.000956887554f, -0.000321084226f, + 0.000620743202f, -0.000487334793f, 0.000260905537f, -9.41662438e-05f, + }, + { + 1.87523437e-05f, -0.000112543086f, 0.000337093195f, -0.000718583527f, + 0.00117377599f, -0.00142052874f, 0.000912307412f, 0.00115821289f, + -0.00575312367f, 0.0137810213f, -0.0258348621f, 0.0419746861f, + -0.0617045611f, 0.0844120756f, -0.111264087f, 0.157740578f, + 0.89528209f, 0.0402862132f, -0.0627515838f, 0.0615863726f, + -0.0523446277f, 0.0400581546f, -0.0276710875f, 0.0170197263f, + -0.00900769606f, 0.00375888613f, -0.000855218444f, -0.000383965147f, + 0.000653946481f, -0.000502010866f, 0.000266113027f, -9.55607466e-05f, + }, + { + 1.8614659e-05f, -0.000111836984f, 0.000333500619f, -0.000706428662f, + 0.00114266353f, -0.00135536294f, 0.000796044304f, 0.00133801904f, + -0.00599371828f, 0.0140505321f, -0.0260562375f, 0.0420071036f, + -0.0613147616f, 0.0832013115f, -0.108381964f, 0.149897799f, + 0.896390736f, 0.0470685475f, -0.0658786595f, 0.0631981492f, + -0.0531094186f, 0.0403239802f, -0.0276620016f, 0.0168859903f, + -0.00884387735f, 0.00361588574f, -0.000752282969f, -0.000447225873f, + 0.000687156862f, -0.000516599219f, 0.000271250232f, -9.69237517e-05f, + }, + { + 1.84634864e-05f, -0.000111076421f, 0.000329776638f, -0.000694045739f, + 0.00111127854f, -0.00129009038f, 0.000680308614f, 0.00151585287f, + -0.00622972287f, 0.0143113835f, -0.0262633823f, 0.0420185439f, + -0.0608970374f, 0.0819577947f, -0.105472468f, 0.142113924f, + 0.897352278f, 0.0539324433f, -0.0690002814f, 0.0647894219f, + -0.0538524352f, 0.040571291f, -0.0276390538f, 0.0167429578f, + -0.0086745061f, 0.00347000291f, -0.000648123911f, -0.000510836719f, + 0.000720356766f, -0.000531090656f, 0.000276313018f, -9.82537313e-05f, + }, + { + 1.82993372e-05f, -0.000110262888f, 0.00032592527f, -0.000681444129f, + 0.00107964023f, -0.00122474553f, 0.000565155409f, 0.00169163849f, + -0.0064610485f, 0.014563499f, -0.0264562722f, 0.0420091227f, + -0.0604517609f, 0.0806823447f, -0.102537185f, 0.134390771f, + 0.898166418f, 0.0608764924f, -0.0721149296f, 0.066359289f, + -0.0545731783f, 0.0407998487f, -0.0276021753f, 0.0165906437f, + -0.0084996298f, 0.00332129095f, -0.000542785216f, -0.00057476759f, + 0.000753528206f, -0.000545475923f, 0.00028129734f, -9.95491428e-05f, + }, + { + 1.81227206e-05f, -0.000109397901f, 0.000321950531f, -0.000668633322f, + 0.00104776816f, -0.00115936319f, 0.000450638996f, 0.00186530151f, + -0.00668760994f, 0.0148068052f, -0.0266349018f, 0.0419789702f, + -0.0599793121f, 0.0793758109f, -0.0995777175f, 0.126730129f, + 0.898832917f, 0.067899242f, -0.0752210692f, 0.0679068193f, + -0.0552711524f, 0.0410094336f, -0.0275513008f, 0.0164290667f, + -0.00831930153f, 0.00316980551f, -0.000436311733f, -0.000638987869f, + 0.000786653196f, -0.000559745939f, 0.000286199123f, -0.000100808451f, + }, + { + 1.79341514e-05f, -0.000108482971f, 0.000317856495f, -0.000655622745f, + 0.00101568119f, -0.00109397725f, 0.000336812896f, 0.00203676918f, + -0.00690932432f, 0.0150412368f, -0.0267992616f, 0.0419282205f, + -0.0594800897f, 0.0780390128f, -0.096595645f, 0.119133741f, + 0.899351537f, 0.0749992207f, -0.0783171728f, 0.0694311261f, + -0.0559458807f, 0.0411998294f, -0.0274863727f, 0.0162582491f, + -0.00813357532f, 0.0030156034f, -0.000328749273f, -0.000703466532f, + 0.000819713518f, -0.000573891331f, 0.000291014294f, -0.000102030113f, + }, + { + 1.77341371e-05f, -0.000107519627f, 0.000313647179f, -0.000642421772f, + 0.000983398757f, -0.00102862192f, 0.000223729847f, 0.00220597023f, + -0.00712611247f, 0.0152667286f, -0.0269493554f, 0.0418570116f, + -0.0589544997f, 0.0766728222f, -0.0935925692f, 0.111603335f, + 0.899722159f, 0.0821749046f, -0.0814016834f, 0.0709313005f, + -0.0565968864f, 0.0413708314f, -0.0274073388f, 0.0160782225f, + -0.00794250984f, 0.00285874377f, -0.000220144575f, -0.000768172264f, + 0.000852690951f, -0.000587902847f, 0.000295738748f, -0.000103212587f, + }, + { + 1.75231889e-05f, -0.000106509382f, 0.000309326686f, -0.000629039947f, + 0.000950939604f, -0.000963330793f, 0.000111441834f, 0.00237283553f, + -0.00733789755f, 0.0154832257f, -0.0270851925f, 0.0417655073f, + -0.0584029555f, 0.0752780885f, -0.0905700624f, 0.104140602f, + 0.899944603f, 0.0894247666f, -0.0844730511f, 0.0724064484f, + -0.0572237112f, 0.0415222421f, -0.0273141582f, 0.0158890188f, + -0.00774616515f, 0.00269928761f, -0.000110545298f, -0.000833073107f, + 0.000885566929f, -0.000601771055f, 0.000300368381f, -0.000104354323f, + }, + { + 1.73018179e-05f, -0.000105453786f, 0.00030489912f, -0.000615486701f, + 0.000918322883f, -0.000898136932f, 5.53861733e-18f, 0.00253729708f, + -0.00754460646f, 0.0156906731f, -0.0272067897f, 0.0416538678f, + -0.0578258894f, 0.0738556832f, -0.0875297114f, 0.0967472121f, + 0.90001893f, 0.0967472121f, -0.0875297114f, 0.0738556832f, + -0.0578258894f, 0.0416538678f, -0.0272067897f, 0.0156906731f, + -0.00754460646f, 0.00253729708f, 5.53861733e-18f, -0.000898136932f, + 0.000918322883f, -0.000615486701f, 0.00030489912f, -0.000105453786f, + }, +}; diff --git a/tests/parity/drift_probe.py b/tests/parity/drift_probe.py new file mode 100644 index 00000000..9d40bd66 --- /dev/null +++ b/tests/parity/drift_probe.py @@ -0,0 +1,93 @@ +"""DriftResampler following a producer on another clock. + + drift_probe.py audiopump + +A Ring is written on a schedule that runs fast or slow against the pulls, by +an integer accumulator, so every target sees the same writes. The node starts +at its target, steers to the producer's rate, rides out a stall (an +underrun) and a moved target, and converts 44.1 kHz up to 48 kHz. Each case +prints a checksum of everything played, the underrun counts and the ratio's +offset rounded to a whole ppm. The material is integer and the kernel is +shared C, so MicroPython, CircuitPython and CPython print the same lines. +""" + +import sys +from array import array + +import audiocore + +MODULE = sys.argv[1] if len(sys.argv) > 1 else "audiopump" +audiopump = __import__(MODULE) + + +def checksum(value, data): + for byte in data: + value = ((value ^ byte) * 16777619) & 0xffffffff + return value + + +class Source: + """An integer saw and its octave, `channels` wide, read in chunks.""" + + def __init__(self, channels, period=97): + self.channels = channels + self.period = period + self.at = 0 + + def take(self, frames): + out = array("h", [0] * (frames * self.channels)) + for i in range(frames): + step = (self.at + i) % self.period + saw = step * 20000 // (self.period - 1) - 10000 + fast = (2 * step) % self.period * 8000 // (self.period - 1) - 4000 + out[i * self.channels] = saw + fast + if self.channels == 2: + out[i * self.channels + 1] = saw - fast + self.at += frames + return out + + +def run(name, channels, in_rate, out_rate, frames, capacity, pulls, + ppm, chunk, stall=None, retarget=None, **options): + """`ppm` fast (positive) or slow against the pulls, in `chunk`-frame + writes. `stall` is a range of pulls with no writes; `retarget` a pull + and the target to set there.""" + ring = audiopump.Ring(sample_rate=in_rate, channel_count=channels, + frames=frames, capacity=capacity) + node = audiopump.DriftResampler(ring, sample_rate=out_rate, **options) + source = Source(channels) + # Input frames owed per pull, times out_rate * 10^6, so the accumulator + # stays an integer at any ppm. + scale = out_rate * 1000000 + per_pull = frames * in_rate * (1000000 + ppm) + owed = 0 + digest = 2166136261 + for pull in range(pulls): + if stall is None or not (stall[0] <= pull < stall[1]): + owed += per_pull + while owed >= chunk * scale: + owed -= chunk * scale + ring.write(source.take(chunk)) + if retarget is not None and pull == retarget[0]: + node.target = retarget[1] + digest = checksum(digest, bytes(audiocore.get_buffer(node)[1])) + if pull % 750 == 749: + print(name, pull + 1, "fill", int(node.fill + 0.5), + "ppm", int(node.ppm + (0.5 if node.ppm >= 0 else -0.5))) + stats = ring.stats() + print(name, "digest", digest, "underruns", node.underruns, stats[3], + "overruns", stats[4], "target", node.target) + + +# The bandwidth is 1/16, exact in binary: a decimal one would rest on each +# interpreter's float parser, and CircuitPython's parses to single precision. +run("stereo_fast", 2, 48000, 48000, 64, 16, 3000, 3000, 48, + max_ppm=5000, bandwidth=0.0625) +run("stereo_slow", 2, 48000, 48000, 64, 16, 3000, -3000, 48, + max_ppm=5000, bandwidth=0.0625) +run("mono_stall", 1, 48000, 48000, 32, 32, 6000, 1000, 40, + stall=(3000, 3060), max_ppm=2000, bandwidth=0.0625) +run("retarget", 2, 48000, 48000, 64, 16, 3000, 0, 64, + retarget=(1500, 300), max_ppm=10000, bandwidth=0.0625) +run("up_44k1", 2, 44100, 48000, 64, 16, 3000, 500, 49, + max_ppm=2000, bandwidth=0.0625) diff --git a/tests/parity/verify_dsp.py b/tests/parity/verify_dsp.py index b27bd804..87749352 100644 --- a/tests/parity/verify_dsp.py +++ b/tests/parity/verify_dsp.py @@ -351,6 +351,9 @@ "comes back empty, and ours return - a deliberate " "departure recorded in docs/how-it-works/upstream-diff.md"}, None), + # audiopump.DriftResampler following a producer on another clock: the + # control loop and the interpolation, three ways. + ("drift_probe.py", "audiopump", {}, None), ) diff --git a/tests/test_binding_parity.py b/tests/test_binding_parity.py index 1c83b99d..2be9f38b 100644 --- a/tests/test_binding_parity.py +++ b/tests/test_binding_parity.py @@ -166,12 +166,13 @@ def test_every_node_type_has_deinit_on_both_sides(self): class TheCircuitPythonRingIsThePumpsRing(unittest.TestCase): - """`audiopump.Ring` is bound twice too: `src/audiopump/audiopump_ring.c` + """`audiopump.Ring` and `audiopump.DriftResampler` are bound twice too: `src/audiopump/audiopump_ring.c` for MicroPython, where the pump pulls it, and `src/circuitpython_spike/shared-bindings/audiopump/Ring.c` for CircuitPython, where the board's audio output does. The point of the CircuitPython one is that a producer written against either runs on both, - so the two must expose the same names. + so the two must expose the same names. DriftResampler is the same: its + MicroPython binding is `src/audiopump/audiopump_drift.c`. It is held here by name rather than through `MODULES`, because the rest of MicroPython's `audiopump` (the pump itself, `Events`, `Tap`) has no @@ -182,24 +183,29 @@ class TheCircuitPythonRingIsThePumpsRing(unittest.TestCase): MICROPYTHON_ONLY = {"from_usb"} def test_the_two_rings_expose_the_same_names(self): - micropython = surfaces([str(ROOT / "src" / "audiopump" / "audiopump_ring.c")], "") + micropython = surfaces([str(ROOT / "src" / "audiopump" / "audiopump_ring.c"), + str(ROOT / "src" / "audiopump" / "audiopump_drift.c")], "") circuitpython = circuitpython_surfaces("audiopump") # The control: both tables are found, so agreement means something. - self.assertEqual(sorted(micropython), ["audiopump_ring"]) - self.assertEqual(sorted(circuitpython), ["audiopump_ring"]) + both = ["audiopump_driftresampler", "audiopump_ring"] + self.assertEqual(sorted(micropython), both) + self.assertEqual(sorted(circuitpython), both) self.assertIn("write", circuitpython["audiopump_ring"]) + self.assertIn("target", circuitpython["audiopump_driftresampler"]) self.assertEqual( micropython["audiopump_ring"] - self.MICROPYTHON_ONLY, circuitpython["audiopump_ring"]) + self.assertEqual(micropython["audiopump_driftresampler"], + circuitpython["audiopump_driftresampler"]) - def test_the_circuitpython_module_carries_only_the_ring(self): + def test_the_circuitpython_module_carries_only_the_ring_and_resampler(self): source = (ROOT / "src" / "circuitpython_spike" / "shared-bindings" / "audiopump" / "__init__.c").read_text() table = re.search(r"globals_table\[\]\s*=\s*\{(.*?)\n\};", source, re.DOTALL) self.assertIsNotNone(table) # The key of each entry, not its value: `__name__`'s value is a qstr too. keys = re.findall(r"\{\s*MP_ROM_QSTR\(MP_QSTR_(\w+)\)", table.group(1)) - self.assertEqual(sorted(keys), ["Ring", "__name__"]) + self.assertEqual(sorted(keys), ["DriftResampler", "Ring", "__name__"]) def test_the_circuitpython_build_compiles_it(self): glue = (ROOT / "circuitpython.mk").read_text() diff --git a/tests/test_cpython_drift.py b/tests/test_cpython_drift.py new file mode 100644 index 00000000..e8cbf015 --- /dev/null +++ b/tests/test_cpython_drift.py @@ -0,0 +1,269 @@ +"""`audiopump.DriftResampler`: lock, hold and sound, against arithmetic. + +A producer on one clock writes a Ring; the node plays it on another and +steers its ratio so the Ring holds its target fill. `drift_probe.py` proves +every target renders the node identically. These are the traits, measured +against the clocks the simulation sets rather than against a stored answer: + +| Trait | Bar | +| --- | --- | +| Locks at +100 ppm and -100 ppm, 48-frame packets | no underrun or overrun; smoothed fill within 8 frames of target after 20 s; mean ratio within 15 ppm of the offset | +| Locks with 480-frame packets (a USB or network burst) | no underrun or overrun; within 8 frames after 45 s | +| Follows a clock wandering +/-100 ppm over a minute | no underrun or overrun; within 80 frames after 20 s | +| THD+N of a 0.9 FS sine through it at +100 ppm, 1 kHz and 10 kHz | below -92 dB (two 16-bit quantizations alone are about -94.5) | +| Constant input | comes out exactly, whatever the fractional position | +| Starts at the target | silence until the Ring reaches it, then the fill is the target | +| A stall | one underrun, counted on the node and the Ring, then plays again | +| Arguments | every bound refused with ValueError, a non-Ring with TypeError | + +The THD+N bar is set by 16-bit audio, not by the filter: a sine quantized +once to 16 bits measures -97.5 dB here and quantized twice about -94.5, and +the resampler adds little to that (docs/how-it-works/drift-resampler.md). +""" + +import math +import unittest + +import numpy as np + +import audiocore +import audiopump + +RATE = 48000 +BLOCK = 256 + + +def blackman_harris(n): + k = np.arange(n) * (2.0 * math.pi / (n - 1)) + return (0.35875 - 0.48829 * np.cos(k) + 0.14128 * np.cos(2 * k) - + 0.01168 * np.cos(3 * k)) + + +def thd_n(y, notch_hz=10.0): + """Everything from 20 Hz to 20 kHz but the fundamental, in dB.""" + spectrum = np.abs(np.fft.rfft(y * blackman_harris(len(y)))) ** 2 + freqs = np.fft.rfftfreq(len(y), 1.0 / RATE) + band = (freqs >= 20) & (freqs <= 20000) + peak = freqs[np.argmax(np.where(band, spectrum, 0))] + fundamental = band & (np.abs(freqs - peak) <= notch_hz) + return 10 * math.log10(spectrum[band & ~fundamental].sum() / + spectrum[fundamental].sum()) + + +def sine_packet(freq, frames, channels=2, amp=0.9): + """`frames` frames of a sine whose period divides them, so every packet + is the same bytes and the stream is continuous.""" + n = np.arange(frames) + s = np.round(amp * 32767 * np.sin(2 * math.pi * freq * n / RATE)) + s = s.astype(" now: + break + due = after + self.ring.write(data) + _, out = audiocore.get_buffer(self.node) + if keep_audio: + audio.append(bytes(out)) + self.fill[block] = self.node.fill + self.ppm[block] = self.node.ppm + self.underruns[block] = self.node.underruns + self.overruns[block] = self.ring.stats()[4] + if keep_audio: + self.audio = np.frombuffer(b"".join(audio), "= len(packet) and ring.level() < 300 * 2: + ring.write(packet) + data = np.frombuffer(bytes(audiocore.get_buffer(node)[1]), " 10: + steady.append(data) + steady = np.concatenate(steady) + self.assertEqual(set(steady.tolist()), {12345}) + self.assertEqual(node.underruns, 0) + + +class TestBehaviour(unittest.TestCase): + + def test_starts_at_target(self): + ring = audiopump.Ring(sample_rate=RATE, channel_count=2, frames=64, + capacity=16) + node = audiopump.DriftResampler(ring, target=400) + packet = sine_packet(1000.0, 48) + silent = bytes(256) + while ring.level() // 4 < 300: + ring.write(packet) + self.assertEqual(bytes(audiocore.get_buffer(node)[1]), silent) + while ring.level() // 4 < 700: + ring.write(packet) + audiocore.get_buffer(node) + # It dropped what was beyond the target and played from there. + self.assertLess(abs(node.fill - 400), 2) + self.assertLess(abs(ring.level() // 4 - (400 - 64)), 20) + self.assertEqual(node.underruns, 0) + + def test_stall(self): + ring = audiopump.Ring(sample_rate=RATE, channel_count=2, frames=64, + capacity=16) + node = audiopump.DriftResampler(ring) + packet = sine_packet(1000.0, 64) + for _ in range(20): + ring.write(packet) + for _ in range(16): + audiocore.get_buffer(node) + self.assertEqual(node.underruns, 1) + self.assertEqual(ring.stats()[3], 1) + for _ in range(12): + ring.write(packet) + data = bytes(audiocore.get_buffer(node)[1]) + self.assertNotEqual(data, bytes(len(data))) + self.assertEqual(node.underruns, 1) + + def test_reset_waits_again(self): + ring = audiopump.Ring(sample_rate=RATE, channel_count=2, frames=64, + capacity=16) + node = audiopump.DriftResampler(ring) + packet = sine_packet(1000.0, 64) + for _ in range(12): + ring.write(packet) + self.assertNotEqual(bytes(audiocore.get_buffer(node)[1]), bytes(256)) + node.reset() + self.assertEqual(bytes(audiocore.get_buffer(node)[1]), bytes(256)) + self.assertEqual(node.ppm, 0.0) + + +class TestArguments(unittest.TestCase): + + def ring(self, **options): + settings = dict(sample_rate=RATE, channel_count=2, frames=256, + capacity=8) + settings.update(options) + return audiopump.Ring(**settings) + + def test_defaults(self): + node = audiopump.DriftResampler(self.ring()) + self.assertEqual(node.target, 1024) + self.assertEqual(node.sample_rate, RATE) + self.assertEqual(node.channel_count, 2) + self.assertEqual(node.underruns, 0) + self.assertAlmostEqual(node.rate, 1.0, places=6) + # Half a two-block Ring is less than a block plus the correction, so + # the default is the least target that works rather than an error. + small = audiopump.DriftResampler(self.ring(capacity=2)) + self.assertEqual(small.target, 258) + + def test_bounds(self): + ring = self.ring() + for options in (dict(sample_rate=44100), dict(sample_rate=96001), + dict(max_ppm=-1), dict(max_ppm=10001), + dict(bandwidth=0.0005), dict(bandwidth=0.2), + dict(target=100), dict(target=2033)): + with self.subTest(**options): + with self.assertRaises(ValueError): + audiopump.DriftResampler(ring, **options) + for frames in (8, 2048): + with self.subTest(frames=frames): + with self.assertRaises(ValueError): + audiopump.DriftResampler(self.ring(frames=frames, + capacity=2)) + with self.assertRaises(TypeError): + audiopump.DriftResampler(object()) + + def test_up_to_twice(self): + node = audiopump.DriftResampler(self.ring(sample_rate=24000), + sample_rate=RATE) + self.assertAlmostEqual(node.rate, 0.5, places=6) + + def test_target_setter(self): + node = audiopump.DriftResampler(self.ring()) + # The ring also holds the filter's 16 frames of look-ahead. + node.target = 2032 + self.assertEqual(node.target, 2032) + for value in (100, 2033): + with self.assertRaises(ValueError): + node.target = value + self.assertEqual(node.target, 2032) + + def test_deinit(self): + node = audiopump.DriftResampler(self.ring()) + node.deinit() + with self.assertRaises(ValueError): + node.ppm + + +if __name__ == "__main__": + unittest.main() diff --git a/tools/design_drift_table.py b/tools/design_drift_table.py new file mode 100644 index 00000000..43cb625e --- /dev/null +++ b/tools/design_drift_table.py @@ -0,0 +1,131 @@ +#!/usr/bin/env python3 +"""Write (or check) the DriftResampler's polyphase interpolation table. + + tools/design_drift_table.py --write regenerate src/shared/audiodsp_drift_table.h + tools/design_drift_table.py --verify check the shipped table is this script's + +`audiopump.DriftResampler` reads its input between samples with a 32-tap +Kaiser-windowed sinc. The kernel is tabulated at 128 fractional positions +(plus the 129th, the first one moved one whole sample, so row p + 1 always +exists) and the resampler interpolates linearly between two neighbouring +rows. Each row is scaled to sum to exactly 1, so a constant input comes out +unchanged at every fractional position. + +The table is a checked-in constant rather than something each runtime +computes, because CircuitPython's floats are single precision and its libm +is not glibc's: a table built at run time would be a different table on +each target, and every target has to render the node identically. + +Pure Python, no numpy: the numbers here are the ones the C ships. +""" + +import argparse +import math +import pathlib +import struct +import sys + +TAPS = 32 +PHASES = 128 +#: The cutoff as a fraction of the input sample rate. 0.45 keeps the passband +#: flat to 20 kHz at 48 kHz within about 1 dB and stops above 0.55. +CUTOFF = 0.45 +#: Kaiser's beta. 9 puts the stopband near -90 dB, below 16-bit audio's floor. +BETA = 9.0 + +ROOT = pathlib.Path(__file__).resolve().parent.parent +HEADER = ROOT / "src" / "shared" / "audiodsp_drift_table.h" + + +def bessel_i0(x): + """The modified Bessel function of order 0, by its power series.""" + total = term = 1.0 + k = 1 + while term > 1e-17 * total: + term *= (x / (2.0 * k)) ** 2 + total += term + k += 1 + return total + + +def kaiser(t): + x = t / (TAPS / 2) + if abs(x) > 1.0: + return 0.0 + return bessel_i0(BETA * math.sqrt(1.0 - x * x)) / bessel_i0(BETA) + + +def sinc(x): + return 1.0 if x == 0 else math.sin(math.pi * x) / (math.pi * x) + + +def f32(value): + return struct.unpack("