diff --git a/book/src/SUMMARY.md b/book/src/SUMMARY.md index 1923c5f..df515c7 100644 --- a/book/src/SUMMARY.md +++ b/book/src/SUMMARY.md @@ -49,3 +49,4 @@ [Appendix A: The C++ decision log](appendix/cpp-decisions.md) [Appendix B: Glossary](appendix/glossary.md) [Appendix C: Annotated bibliography](appendix/bibliography.md) +[Appendix D: The two budgets](appendix/two-budgets.md) diff --git a/book/src/appendix/two-budgets.md b/book/src/appendix/two-budgets.md new file mode 100644 index 0000000..3758270 --- /dev/null +++ b/book/src/appendix/two-budgets.md @@ -0,0 +1,187 @@ +# Appendix D: The two budgets + +> We shape our tools, and thereafter our tools shape us. +> +> — John M. Culkin (often attributed to Marshall McLuhan) + +*This appendix speaks in the first person: the surprise it describes was +mine. — T.* + +## The sentence that stopped me + +In the middle of the mailing-list exchange that became this book's +epilogue, Robert Bristow-Johnson mentioned, in passing, how he planned to +fix his filter-design script: "I can modify that by using an FFT (with a +really large size N = 2^20 = 1048576)." + +I stopped on that sentence. I have been writing real-time audio software +for thirty years, and my honest, immediate reaction was: *I have never +heard of anyone running an FFT that large.* A million points. It sounded +like something you'd book cluster time for. + +Here is what a million-point FFT actually costs, measured on the modest +shared virtual machine this book's figures are built on — not a good +machine; the point survives the handicap: + +| size | points | time | input | +|---|---|---|---| +| 2²⁰ | 1,048,576 | 165 ms | 8.4 MB | +| 2²¹ | 2,097,152 | 339 ms | 16.8 MB | +| 2²⁴ | 16,777,216 | 3.4 s | 134 MB | + +On a decent laptop with a tuned FFT library, divide by five or ten. The +"really large" FFT is a *sip of coffee*. And the embarrassing coda: this +project had been running 2²¹-point FFTs all along — every response +measurement in the verification notebook, every filter curve in the +book's figures. The huge FFT was in the room the whole time. I had +written the code that called it and never once registered its size as a +number that meant anything. + +That gap — between what the number *sounded like* and what it *cost* — +is worth an appendix, because it isn't ignorance. It's training. + +## Where the instinct comes from + +I grew up, professionally, inside the perform loop. In the Max/MSP world +the signal vector is 64 samples long; at 44.1 kHz that is 1.45 +milliseconds, and everything — every object, every patch, every product — +is priced in what it does to that deadline. You learn to feel a +per-sample multiply the way a carpenter feels the weight of a hammer. An +FFT, in that world, is not "an algorithm"; it is a *latency decision*: a +4096-point FFT means 93 milliseconds of buffering before the first +useful output, so you don't run one unless the musical idea demands it, +and you certainly never think about a size with seven digits, because +seven digits is twenty-two *seconds* of audio and no deadline that +matters is twenty-two seconds away. The thought "a million-point FFT" is +not rejected by that mindset. It simply never forms. + +Thirty years of that is not a bias you notice. It is a currency you +think in. Every computation gets priced in the only unit that ever +mattered — time per sample, paid at the deadline — including, silently, +the computations that never go anywhere near a deadline. + +Robert's vocabulary has the same shape with an older mint mark. He built +his instincts on SHArCs three decades ago, when a million points of +*memory* — eight megabytes of doubles — was the obstacle, never mind the +cycles. When he says "really large," he is speaking 1995 fluently. When +I hear "really large" and flinch, I am speaking perform-loop fluently. +We are both correct in a currency that no longer applies to the +transaction: his script runs once, offline, in MATLAB, on a laptop. + +## The two budgets + +So name the thing properly. Audio software lives under two budgets, and +they have almost nothing in common: + +**The sample budget.** Paid per sample, at a deadline, forever. Miss it +once and everyone hears the click. This budget is hard in both senses — +unforgiving and difficult — and it deserves every instinct the last +thirty years installed: count the multiplies, fear the allocation, +distrust the branch. In this library the sample budget governs exactly +two functions, `push()` and `pull()`, and the numbers it runs to are +a couple hundred nanoseconds per stereo frame. + +**The construction budget.** Paid once, off the deadline. Design a +filter, fit a compensation curve, solve a least-squares system, verify a +claim with a million-point FFT — none of it happens while anyone is +listening. Work here is not free, but it is priced in a different unit: +*is this fast enough to happen at startup, or in a design session, +without annoying a human?* Milliseconds are pennies. Whole seconds are +sometimes fine. + +The line between the two budgets is drawn straight through this +library's API, and always was: the constructor may allocate, may throw, +may spend milliseconds on Bessel functions; `pull()` is `noexcept`, +lock-free, allocation-free, and answers to the deadline. Half of this +book has been about defending that line. What the FFT moment taught me +is that I had internalized the line for *code* but not for *thinking* — +my design-time reasoning was still spending sample-budget currency, +thirty years after it stopped being the right money. + +## What "large" costs where + +The two budgets don't just price the same work differently — they +sometimes invert it, and this project has measured both directions. + +Robert's million-point FFT is the cheap direction. His design method +draws the ideal response on an N-point frequency grid and inverse-FFTs +it; N must be much larger than the filter so the grid is dense and the +circular wraparound has room to die out — his 2²⁰ over a 16,384-tap +filter is the standard 64× headroom, not extravagance. Cost: O(N log N), +about twenty million operations, a sip of coffee. Meanwhile, in the same +email, his `firls()` call "takes several seconds" and `firpm()` "gets +funky" at the same 16K taps — because dense least-squares design scales +like the *square* of the filter length and Remez exchange grows fragile +long before that. Offline is not one place. It has its own O(·) +geography, and the FFT lives in the flattest part of it. + +The inversion runs the other way on small machines, and the +instruction-count ledger caught this project relying on the cheap-side +assumption. The compensated filter design costs roughly seven million +double-precision operations — nothing, I would have said, it's the +construction budget. Then the Cortex-M33 leg of CI priced it: on a +target where every double multiply is a ~140-instruction library call, +"nothing" measured **1.9 billion instructions** — seconds of boot on a +part that might be inside a product with a power button. The design got +resized (one correction pass, fewer probes) to fit a budget I hadn't +known it was spending against. Same work, three targets, three prices: +16 million instructions on the M55, 131 million on Hexagon, 929 million +on the M33. "Construction-time" names *when* the money is spent, not +*how much there is*. + +And the cosine-series trick at the heart of the compensated design is, +seen this way, a currency exchange. Robert pays for his passband tilt +with a million-point FFT, because MATLAB makes FFTs free and his script +runs on a laptop. This library pays for the identical tilt with a +handful of shifted sinc evaluations and a 15-unknown least-squares +solve, because a dependency-free C++ header makes FFTs expensive and its +constructor sometimes runs on a microcontroller. Two implementations of +one filter, each shaped by which budget its environment discounts. + +## What the instinct is still for + +It would be the wrong reading of this appendix to conclude that the +old reflexes are obsolete. The sample budget is exactly as merciless as +it was in 1995; the perform loop still doesn't care about your feelings; +every instinct that flinches at an allocation inside `pull()` is earning +its keep tonight in somebody's live rig. The reflexes are not wrong. +They are *scoped*, and the scope is the deadline. + +The skill I apparently still get to learn, thirty years in, is noticing +which budget I'm standing in before I price the work. A million-point +FFT at the deadline is absurd. The same FFT in a design script is a +rounding error, and treating it as absurd there has real costs — it +makes you avoid the honest tool, under-verify the finished design, and +mistake an expert's era-stamped vocabulary ("really large") for a +warning when it was just an old accent, faithfully kept. + +The samples still have to arrive on time. Everything else has all the +time in the world. + +## Verify it yourself + +```sh +# The table above, on your machine (any Python with numpy): +python3 -c " +import numpy as np, time +for p in (20, 21, 24): + x = np.random.default_rng(1).standard_normal(1 << p) + t0 = time.perf_counter(); np.fft.rfft(x); t1 = time.perf_counter() + print(f'2^{p}: {1e3*(t1-t0):8.1f} ms {x.nbytes/1e6:7.1f} MB')" + +# The million-point-class FFTs this book was already running: every +# spec/response measurement in the verification notebook uses rfft at 2^21. +grep -n '1 << 21' notebooks/asrc_rbj_analysis.ipynb scripts/book_figures.py + +# The two budgets, drawn through the API: the constructor may throw and +# allocate; the audio path may not. The real-time contract is a grep away: +grep -n 'noexcept' include/srt/asrc.h | head + +# What the construction budget costs where it is NOT cheap: the M33 entry +# in the instruction-count ledger. +grep -n 'M33' docs/PERFORMANCE.md | head +``` + +Wall-clock caveat, as always in this book: the timings are one machine's +and the table's job is the order of magnitude, not the digits. That is +the entire point. diff --git a/book/src/epilogue/letter.md b/book/src/epilogue/letter.md index e28d206..6b1bdde 100644 --- a/book/src/epilogue/letter.md +++ b/book/src/epilogue/letter.md @@ -227,6 +227,37 @@ recorded in the notebook as an honest draw, with the multichannel case still settled by C1's measurements. Not every good question has a portable answer; the record is the deliverable. +## Postscript: the script arrives + +A week after the exchange, his MATLAB script did arrive — updated the same +day so the Kaiser-windowed method got the sinc² zeros too, and validated, +by his own report, at R=512 and N=32: the exact shape of `economy()`. With +it came one more claim, stated with his usual precision: *"for DC, this +should get you infinite S/N ratio... for every phase or fractional delay, +the FIR coefficients must add to 1."* + +Checked, and better than checked. In double precision the property falls +out of the rect construction *automatically*: the compensated designs' +branch sums are uniform to a spread of 1.8×10⁻¹⁵ — machine epsilon — +because zeros at k·fs and branch-DC uniformity are the same fact stated in +two domains (the plain design's 4.7×10⁻⁶ spread is its stopband leakage at +fs, visible from the other side). But his accompanying warning — that +16-bit coefficient quantization needs "tricky things" — turned out to be +pointing at a live defect: independent per-tap rounding was re-breaking +the property by several LSB, a bias this suite had measured and papered +over with a widened tolerance earlier in this very chapter's story. The +tricky thing is row-sum-preserving quantization: + +```cpp +{{#include ../../../include/srt/polyphase_filter.h:pw_row_sum}} +``` + +After it, the measured result is the one his claim demanded: a DC input +passes through the Q15 and Q31 converters **bit-exactly at every +fractional delay** — zero LSB of deviation over a 256-point μ sweep — and +`FixedPoint.RowSumsAreExact*` pins every row's sum to exactly one, in +coefficient units, forever. Infinite S/N at DC, delivered in 16 bits. + ## What this chapter is actually about Strip the DSP away and the shape of the episode is this: an expert @@ -261,6 +292,10 @@ ctest --test-dir build -R 'ProgramWeighted' --output-on-failure # The half-sample alignment sentinel that caught the rect-centering bug: ctest --test-dir build -R 'FractionalDelay' --output-on-failure + +# The postscript's claims: branch-DC uniformity at machine epsilon, exact +# fixed-point row sums, bit-exact DC at every fractional delay: +ctest --test-dir build -R 'BranchSums|RowSums|DcGain' --output-on-failure ``` And one experiment in the spirit of the thread: change `image_zeros` to diff --git a/examples/alsa_bridge.cpp b/examples/alsa_bridge.cpp index 8b36173..ff3b242 100644 --- a/examples/alsa_bridge.cpp +++ b/examples/alsa_bridge.cpp @@ -34,35 +34,35 @@ namespace { - std::atomic gStop{false}; + std::atomic g_stop{false}; - void onSigint(int) { - gStop.store(true, std::memory_order_relaxed); + void on_sigint(int) { + g_stop.store(true, std::memory_order_relaxed); } - const char* stateName(srt::converter_state s) { + const char* state_name(srt::converter_state s) { switch (s) { - case srt::converter_state::Filling: + case srt::converter_state::filling: return "Filling"; - case srt::converter_state::Acquiring: + case srt::converter_state::acquiring: return "Acquiring"; - case srt::converter_state::Locked: + case srt::converter_state::locked: return "Locked"; } return "?"; } - struct Args { - const char* inDev = "default"; - const char* outDev = "default"; - unsigned rate = 48000; - unsigned channels = 2; - snd_pcm_uframes_t period = 128; - std::size_t latency = 96; - const char* csvPath = nullptr; - const char* dumpPath = nullptr; - unsigned long seconds = 0; // 0 = run until SIGINT - double toneHz = 0.0; // 0 = pass captured audio through + struct cli_args { + const char* in_dev = "default"; + const char* out_dev = "default"; + unsigned rate = 48000; + unsigned channels = 2; + snd_pcm_uframes_t period = 128; + std::size_t latency = 96; + const char* csv_path = nullptr; + const char* dump_path = nullptr; + unsigned long seconds = 0; // 0 = run until SIGINT + double tone_hz = 0.0; // 0 = pass captured audio through }; void usage(const char* prog) { @@ -72,7 +72,7 @@ namespace { " --rate nominal sample rate of both devices (default 48000)\n" " --channels channel count (default 2)\n" " --period frames per ALSA period (default 128)\n" - " --latency converter targetLatencyFrames (default 96)\n" + " --latency converter target_latency_frames (default 96)\n" " --csv append per-second telemetry CSV\n" " --dump write post-ASRC interleaved float stream raw\n" " --seconds run time in seconds (default 0 = until SIGINT)\n" @@ -81,7 +81,7 @@ namespace { prog); } - bool parseArgs(int argc, char** argv, Args& a) { + bool parse_args(int argc, char** argv, cli_args& a) { const auto value = [&](int& i) -> const char* { if (i + 1 >= argc) { std::fprintf(stderr, "%s: missing value for %s\n", argv[0], argv[i]); @@ -131,21 +131,21 @@ namespace { return true; } - struct AlsaDevice { - snd_pcm_t* pcm = nullptr; - snd_pcm_format_t format = SND_PCM_FORMAT_UNKNOWN; - snd_pcm_uframes_t periodFrames = 0; + struct alsa_device { + snd_pcm_t* pcm = nullptr; + snd_pcm_format_t format = SND_PCM_FORMAT_UNKNOWN; + snd_pcm_uframes_t period_frames = 0; - AlsaDevice() = default; - AlsaDevice(const AlsaDevice&) = delete; - AlsaDevice& operator=(const AlsaDevice&) = delete; - ~AlsaDevice() { + alsa_device() = default; + alsa_device(const alsa_device&) = delete; + alsa_device& operator=(const alsa_device&) = delete; + ~alsa_device() { if (pcm != nullptr) snd_pcm_close(pcm); } }; - bool openDevice(AlsaDevice& dev, const char* name, snd_pcm_stream_t stream, const Args& a) { + bool open_device(alsa_device& dev, const char* name, snd_pcm_stream_t stream, const cli_args& a) { const char* dir = stream == SND_PCM_STREAM_CAPTURE ? "capture" : "playback"; int err = snd_pcm_open(&dev.pcm, name, stream, 0); if (err < 0) { @@ -179,16 +179,16 @@ namespace { } dev.period_frames = a.period; int sub = 0; - if ((err = snd_pcm_hw_params_set_period_size_near(dev.pcm, hw, &dev.periodFrames, &sub)) < 0) + if ((err = snd_pcm_hw_params_set_period_size_near(dev.pcm, hw, &dev.period_frames, &sub)) < 0) return fail("set period size", err); - snd_pcm_uframes_t bufFrames = 4 * dev.periodFrames; - if ((err = snd_pcm_hw_params_set_buffer_size_near(dev.pcm, hw, &bufFrames)) < 0) + snd_pcm_uframes_t buf_frames = 4 * dev.period_frames; + if ((err = snd_pcm_hw_params_set_buffer_size_near(dev.pcm, hw, &buf_frames)) < 0) return fail("set buffer size", err); if ((err = snd_pcm_hw_params(dev.pcm, hw)) < 0) return fail("hw_params commit", err); std::printf("%s '%s': %s, %u Hz, %u ch, period %lu, buffer %lu\n", dir, name, snd_pcm_format_name(dev.format), - rate, a.channels, static_cast(dev.periodFrames), - static_cast(bufFrames)); + rate, a.channels, static_cast(dev.period_frames), + static_cast(buf_frames)); return true; } @@ -207,8 +207,8 @@ namespace { } // namespace int main(int argc, char** argv) { - Args args; - if (!parseArgs(argc, argv, args)) { + cli_args args; + if (!parse_args(argc, argv, args)) { usage(argv[0]); return 2; } @@ -217,13 +217,13 @@ int main(int argc, char** argv) { return 2; } - AlsaDevice in; - AlsaDevice out; - if (!openDevice(in, args.inDev, SND_PCM_STREAM_CAPTURE, args) - || !openDevice(out, args.outDev, SND_PCM_STREAM_PLAYBACK, args)) + alsa_device in; + alsa_device out; + if (!open_device(in, args.in_dev, SND_PCM_STREAM_CAPTURE, args) + || !open_device(out, args.out_dev, SND_PCM_STREAM_PLAYBACK, args)) return 1; - srt::Config cfg; + srt::config cfg; cfg.sample_rate_hz = static_cast(args.rate); cfg.channels = args.channels; cfg.target_latency_frames = args.latency; @@ -231,46 +231,46 @@ int main(int argc, char** argv) { // comfortably above half the transfer block, or block-quantized // occupancy excursions can demote the servo stage spuriously. cfg.servo.unlock_threshold_frames = - std::max(cfg.servo.unlockThresholdFrames, 1.5 * static_cast(args.period)); + std::max(cfg.servo.unlock_threshold_frames, 1.5 * static_cast(args.period)); srt::async_sample_rate_converter asrc(cfg); - std::printf("designed latency: %.2f ms%s\n", asrc.designedLatencySeconds() * 1e3, - args.toneHz > 0.0 ? " (tone mode: captured samples discarded)" : ""); + std::printf("designed latency: %.2f ms%s\n", asrc.designed_latency_seconds() * 1e3, + args.tone_hz > 0.0 ? " (tone mode: captured samples discarded)" : ""); std::FILE* csv = nullptr; - if (args.csvPath != nullptr) { - csv = std::fopen(args.csvPath, "a"); + if (args.csv_path != nullptr) { + csv = std::fopen(args.csv_path, "a"); if (csv == nullptr) { - std::fprintf(stderr, "cannot open CSV file '%s'\n", args.csvPath); + std::fprintf(stderr, "cannot open CSV file '%s'\n", args.csv_path); return 1; } std::fprintf(csv, "time_s,state,ppm,fill,underruns,overruns,resyncs\n"); } std::FILE* dump = nullptr; - if (args.dumpPath != nullptr) { - dump = std::fopen(args.dumpPath, "wb"); + if (args.dump_path != nullptr) { + dump = std::fopen(args.dump_path, "wb"); if (dump == nullptr) { - std::fprintf(stderr, "cannot open dump file '%s'\n", args.dumpPath); + std::fprintf(stderr, "cannot open dump file '%s'\n", args.dump_path); return 1; } } - std::signal(SIGINT, onSigint); + std::signal(SIGINT, on_sigint); std::thread capture([&] { const std::size_t ch = args.channels; - const snd_pcm_uframes_t period = in.periodFrames; + const snd_pcm_uframes_t period = in.period_frames; std::vector raw(period * ch); std::vector buf(period * ch); double phase = 0.0; - const double dphi = 2.0 * std::numbers::pi * args.toneHz / static_cast(args.rate); - while (!gStop.load(std::memory_order_relaxed)) { + const double dphi = 2.0 * std::numbers::pi * args.tone_hz / static_cast(args.rate); + while (!g_stop.load(std::memory_order_relaxed)) { void* dst = in.format == SND_PCM_FORMAT_S16_LE ? static_cast(raw.data()) : static_cast(buf.data()); const snd_pcm_sframes_t n = snd_pcm_readi(in.pcm, dst, period); if (n < 0) { if (snd_pcm_recover(in.pcm, static_cast(n), 1) < 0) { std::fprintf(stderr, "capture failed: %s\n", snd_strerror(static_cast(n))); - gStop.store(true, std::memory_order_relaxed); + g_stop.store(true, std::memory_order_relaxed); return; } continue; @@ -278,7 +278,7 @@ int main(int argc, char** argv) { const auto frames = static_cast(n); if (frames == 0) continue; - if (args.toneHz > 0.0) { // tone mode: keep the device pacing, swap the payload + if (args.tone_hz > 0.0) { // tone mode: keep the device pacing, swap the payload for (std::size_t f = 0; f < frames; ++f) { const auto v = static_cast(0.5 * std::sin(phase)); phase += dphi; @@ -297,21 +297,21 @@ int main(int argc, char** argv) { std::thread playback([&] { const std::size_t ch = args.channels; - const snd_pcm_uframes_t period = out.periodFrames; + const snd_pcm_uframes_t period = out.period_frames; std::vector buf(period * ch); std::vector raw(period * ch); - bool dumpFailed = false; - while (!gStop.load(std::memory_order_relaxed)) { + bool dump_failed = false; + while (!g_stop.load(std::memory_order_relaxed)) { asrc.pull(buf.data(), period); // silence-pads while filling/underrun - if (dump != nullptr && !dumpFailed + if (dump != nullptr && !dump_failed && std::fwrite(buf.data(), sizeof(float), period * ch, dump) != period * ch) { std::fprintf(stderr, "dump write failed; disabling --dump\n"); - dumpFailed = true; + dump_failed = true; } if (out.format == SND_PCM_FORMAT_S16_LE) floatToS16(buf.data(), raw.data(), period * ch); snd_pcm_uframes_t done = 0; - while (done < period && !gStop.load(std::memory_order_relaxed)) { + while (done < period && !g_stop.load(std::memory_order_relaxed)) { const void* src = out.format == SND_PCM_FORMAT_S16_LE ? static_cast(raw.data() + done * ch) : static_cast(buf.data() + done * ch); @@ -319,7 +319,7 @@ int main(int argc, char** argv) { if (n < 0) { if (snd_pcm_recover(out.pcm, static_cast(n), 1) < 0) { std::fprintf(stderr, "playback failed: %s\n", snd_strerror(static_cast(n))); - gStop.store(true, std::memory_order_relaxed); + g_stop.store(true, std::memory_order_relaxed); return; } continue; @@ -332,26 +332,27 @@ int main(int argc, char** argv) { using clock = std::chrono::steady_clock; const auto t0 = clock::now(); - for (unsigned long sec = 1; !gStop.load(std::memory_order_relaxed); ++sec) { + for (unsigned long sec = 1; !g_stop.load(std::memory_order_relaxed); ++sec) { std::this_thread::sleep_until(t0 + std::chrono::seconds(sec)); - if (gStop.load(std::memory_order_relaxed)) + if (g_stop.load(std::memory_order_relaxed)) break; const auto st = asrc.status(); std::printf("t=%6lus state=%-9s ppm=%+8.2f fill=%8.1f under=%llu over=%llu " "resync=%llu\n", - sec, stateName(st.state), st.ppm, st.fifoFillFrames, static_cast(st.underruns), - static_cast(st.overruns), static_cast(st.resyncs)); + sec, state_name(st.state), st.ppm, st.fifo_fill_frames, + static_cast(st.underruns), static_cast(st.overruns), + static_cast(st.resyncs)); std::fflush(stdout); if (csv != nullptr) { - std::fprintf(csv, "%lu,%s,%.3f,%.2f,%llu,%llu,%llu\n", sec, stateName(st.state), st.ppm, st.fifoFillFrames, - static_cast(st.underruns), static_cast(st.overruns), - static_cast(st.resyncs)); + std::fprintf(csv, "%lu,%s,%.3f,%.2f,%llu,%llu,%llu\n", sec, state_name(st.state), st.ppm, + st.fifo_fill_frames, static_cast(st.underruns), + static_cast(st.overruns), static_cast(st.resyncs)); std::fflush(csv); } if (args.seconds != 0 && sec >= args.seconds) break; } - gStop.store(true, std::memory_order_relaxed); + g_stop.store(true, std::memory_order_relaxed); capture.join(); playback.join(); diff --git a/examples/drifting_clocks.cpp b/examples/drifting_clocks.cpp index 9973c4a..da640ab 100644 --- a/examples/drifting_clocks.cpp +++ b/examples/drifting_clocks.cpp @@ -29,10 +29,10 @@ namespace { - constexpr double k_k_fs = 48000.0; - constexpr double k_k_consumer_ppm = 500.0; // consumer clock runs 500 ppm fast - constexpr std::size_t k_k_chunk = 96; - constexpr double k_k_run_seconds = 20.0; + constexpr double k_fs = 48000.0; + constexpr double k_consumer_ppm = 500.0; // consumer clock runs 500 ppm fast + constexpr std::size_t k_chunk = 96; + constexpr double k_run_seconds = 20.0; const char* state_name(srt::converter_state s) { switch (s) { @@ -56,7 +56,7 @@ int main() { cfg.servo.unlock_threshold_frames = 96.0; srt::async_sample_rate_converter asrc(cfg); - std::printf("drifting_clocks: producer 48000.0 Hz, consumer %+.0f ppm, %g s\n", k_k_consumer_ppm, k_k_run_seconds); + std::printf("drifting_clocks: producer 48000.0 Hz, consumer %+.0f ppm, %g s\n", k_consumer_ppm, k_run_seconds); std::printf("designed latency: %.2f ms\n", asrc.designed_latency_seconds() * 1e3); std::atomic stop{false}; @@ -64,30 +64,30 @@ int main() { const auto t0 = clock::now(); std::thread producer([&] { - std::vector buf(k_k_chunk); - const double nu = 997.0 / k_k_fs; + std::vector buf(k_chunk); + const double nu = 997.0 / k_fs; std::uint64_t idx = 0; auto next = t0; const auto period = std::chrono::duration_cast( - std::chrono::duration(static_cast(k_k_chunk) / k_k_fs)); + std::chrono::duration(static_cast(k_chunk) / k_fs)); while (!stop.load(std::memory_order_relaxed)) { for (auto& v : buf) { v = static_cast(0.5 * std::sin(2.0 * std::numbers::pi * nu * static_cast(idx++))); } - asrc.push(buf.data(), k_k_chunk); + asrc.push(buf.data(), k_chunk); next += period; std::this_thread::sleep_until(next); } }); std::thread consumer([&] { - std::vector buf(k_k_chunk); - const double fs_out = k_k_fs * (1.0 + k_k_consumer_ppm * 1e-6); + std::vector buf(k_chunk); + const double fs_out = k_fs * (1.0 + k_consumer_ppm * 1e-6); auto next = t0; const auto period = std::chrono::duration_cast( - std::chrono::duration(static_cast(k_k_chunk) / fs_out)); + std::chrono::duration(static_cast(k_chunk) / fs_out)); while (!stop.load(std::memory_order_relaxed)) { - asrc.pull(buf.data(), k_k_chunk); + asrc.pull(buf.data(), k_chunk); next += period; std::this_thread::sleep_until(next); } @@ -95,7 +95,7 @@ int main() { double ppm_avg = 0.0; const double avg_alpha = 0.5 / 3.0; // 0.5 s prints, ~3 s averaging - for (double t = 0.5; t <= k_k_run_seconds; t += 0.5) { + for (double t = 0.5; t <= k_run_seconds; t += 0.5) { std::this_thread::sleep_for(std::chrono::milliseconds(500)); const auto st = asrc.status(); ppm_avg += avg_alpha * (st.ppm - ppm_avg); @@ -113,9 +113,8 @@ int main() { const auto st = asrc.status(); // The producer pushes at exactly 48000 Hz and the consumer pulls 500 ppm // fast, so the converter must consume input slower than unity: -500 ppm. - std::printf("\nfinal: state=%s ppm(avg)=%+.2f (expected ~%+.0f)\n", state_name(st.state), ppm_avg, - -k_k_consumer_ppm); - const bool ok = st.state == srt::converter_state::locked && std::abs(ppm_avg + k_k_consumer_ppm) < 150.0; + std::printf("\nfinal: state=%s ppm(avg)=%+.2f (expected ~%+.0f)\n", state_name(st.state), ppm_avg, -k_consumer_ppm); + const bool ok = st.state == srt::converter_state::locked && std::abs(ppm_avg + k_consumer_ppm) < 150.0; std::printf("%s\n", ok ? "OK" : "NOT CONVERGED (heavily loaded machine?)"); return ok ? 0 : 1; } diff --git a/include/srt/asrc.h b/include/srt/asrc.h index 0bff68c..ab24885 100644 --- a/include/srt/asrc.h +++ b/include/srt/asrc.h @@ -102,7 +102,7 @@ namespace srt { explicit basic_async_sample_rate_converter(const config& cfg) : m_cfg(validated(cfg)) , m_bank(m_cfg.filter, m_cfg.sample_rate_hz) - , m_resampler(m_bank, m_cfg.channels, k_k_pop_chunk_frames) + , m_resampler(m_bank, m_cfg.channels, k_pop_chunk_frames) , m_ring(ring_capacity_elems(m_cfg, m_bank.taps())) , m_servo(m_cfg.servo, m_cfg.sample_rate_hz, static_cast(m_cfg.target_latency_frames)) , m_target_frames(m_cfg.target_latency_frames) @@ -166,7 +166,7 @@ namespace srt { // the entry occupancy never grazes the pull size; configs that // already satisfy it (e.g. the 32-frame default transfer against // the 48-frame default setpoint) are left exactly as configured. - const std::size_t needed = frames + std::max(frames / 2, k_k_pop_chunk_frames); + const std::size_t needed = frames + std::max(frames / 2, k_pop_chunk_frames); const std::size_t new_target = std::clamp(needed, m_cfg.target_latency_frames, m_max_target_frames); if (new_target > m_target_frames) { m_target_frames = new_target; @@ -193,7 +193,7 @@ namespace srt { occ = backlog_frames(); m_servo.seed(occ); m_filling = false; - m_fade_frames_left = k_k_fade_frames; + m_fade_frames_left = k_fade_frames; } // ANCHOR_END: asrc_filling @@ -271,7 +271,7 @@ namespace srt { const polyphase_filter_bank& filter_bank() const noexcept { return m_bank; } private: - static constexpr std::size_t k_k_pop_chunk_frames = 16; + static constexpr std::size_t k_pop_chunk_frames = 16; static std::size_t ring_capacity_elems(const config& cfg, std::size_t taps) { const std::size_t fill_threshold = cfg.target_latency_frames + taps; @@ -310,9 +310,9 @@ namespace srt { /// even on FPU-less targets. void apply_fade_in(S* interleaved, std::size_t made_frames) noexcept { const std::size_t n = std::min(made_frames, m_fade_frames_left); - const std::size_t done = k_k_fade_frames - m_fade_frames_left; + const std::size_t done = k_fade_frames - m_fade_frames_left; for (std::size_t f = 0; f < n; ++f) { - const double g = static_cast(done + f + 1) / static_cast(k_k_fade_frames); + const double g = static_cast(done + f + 1) / static_cast(k_fade_frames); for (std::size_t c = 0; c < m_cfg.channels; ++c) { S& x = interleaved[f * m_cfg.channels + c]; x = scale_sample(x, g); @@ -370,7 +370,7 @@ namespace srt { return cfg; } - static constexpr std::size_t k_k_fade_frames = 64; + static constexpr std::size_t k_fade_frames = 64; config m_cfg; polyphase_filter_bank m_bank; diff --git a/include/srt/polyphase_filter.h b/include/srt/polyphase_filter.h index 10762a3..fe4f3c9 100644 --- a/include/srt/polyphase_filter.h +++ b/include/srt/polyphase_filter.h @@ -191,12 +191,56 @@ namespace srt { } m_table.resize((m_phases + 1) * m_taps); - for (std::size_t p = 0; p <= m_phases; ++p) { + std::vector remainder(m_taps); + for (std::size_t p = 0; p < m_phases; ++p) { + coeff* row = m_table.data() + p * m_taps; + double exact_sum = 0.0; + std::int64_t quant_sum = 0; for (std::size_t t = 0; t < m_taps; ++t) { - const std::size_t m = t * m_phases + p; // prototype index of (branch p, tap t) - const double v = (m < n) ? proto[m] : 0.0; - m_table[p * m_taps + (m_taps - 1 - t)] = sample_traits::make_coeff(v); + const std::size_t m = t * m_phases + p; // prototype index of (branch p, tap t) + const double v = (m < n) ? proto[m] : 0.0; + const double scaled = v * sample_traits::k_coeff_scale; + const coeff q = sample_traits::make_coeff(v); + row[m_taps - 1 - t] = q; + remainder[m_taps - 1 - t] = scaled - static_cast(q); + exact_sum += scaled; + quant_sum += static_cast(q); } + if constexpr (!std::is_floating_point_v) { + // ANCHOR: pw_row_sum + // Row-sum-preserving quantization ("the coefficients of every + // phase must add to one" — R. Bristow-Johnson, music-dsp). In + // double the image_zeros designs make every branch's DC sum + // identical to machine epsilon (the k*fs zeros ARE branch-DC + // uniformity, stated in frequency), but independent per-tap + // rounding rebroke it by several LSB. Distribute each row's + // total rounding residual to the taps that were rounded + // furthest from it (largest-remainder method): every row then + // sums to llround(exact * scale), so DC gain deviates by at + // most one coefficient LSB across all mu. + std::int64_t residual = static_cast(std::llround(exact_sum)) - quant_sum; + while (residual != 0) { + const double sgn = residual > 0 ? 1.0 : -1.0; + std::size_t best = 0; + for (std::size_t u = 1; u < m_taps; ++u) { + if (sgn * remainder[u] > sgn * remainder[best]) { + best = u; + } + } + row[best] = static_cast(row[best] + (residual > 0 ? 1 : -1)); + remainder[best] -= sgn; + residual -= residual > 0 ? 1 : -1; + } + // ANCHOR_END: pw_row_sum + } + } + // The extra row L is row 0 advanced one input sample; copy it from + // the QUANTIZED row 0 so the mu-wrap continuity invariant holds + // bit-exactly even after the row-sum correction above. + coeff* row_l = m_table.data() + m_phases * m_taps; + row_l[0] = coeff{0}; + for (std::size_t u = 1; u < m_taps; ++u) { + row_l[u] = m_table[u - 1]; } } // ANCHOR_END: bank_build @@ -426,7 +470,7 @@ namespace srt { public: /// Frame-major channel-parallel mode is compiled in only on CP targets /// and only for floating-point samples (see SRT_CHANNEL_PARALLEL). - static constexpr bool k_k_channel_parallel = SRT_CHANNEL_PARALLEL != 0 && std::is_floating_point_v; + static constexpr bool k_channel_parallel = SRT_CHANNEL_PARALLEL != 0 && std::is_floating_point_v; /// Allocates histories and the pop scratch buffer; setup time only. fractional_resampler(const polyphase_filter_bank& bank, std::size_t channels, std::size_t chunk_frames = 64) @@ -435,7 +479,7 @@ namespace srt { , m_chunk(chunk_frames) , m_hist_cap(bank.taps() + chunk_frames) , m_scratch(chunk_frames * channels) - , m_frame_major(k_k_channel_parallel && channels >= SRT_CP_MIN_CHANNELS) + , m_frame_major(k_channel_parallel && channels >= SRT_CP_MIN_CHANNELS) , m_hist(m_frame_major ? 1 : channels) , m_row(bank.taps()) { if (m_channels == 0 || m_chunk == 0) { @@ -531,7 +575,7 @@ namespace srt { if (m_channels == 1 && !k_prefer_dot_row) { // fused blend+mac; no scratch traffic out[n] = interpolate_phase(*m_bank, window(0), m); } - else if (k_k_channel_parallel && m_frame_major) { // constant-folds away off-host + else if (k_channel_parallel && m_frame_major) { // constant-folds away off-host // High channel counts: one blend, then all channels' dots in // a single channel-parallel pass over the frame-major window. blend_row_phase(*m_bank, m_row.data(), m); @@ -572,14 +616,14 @@ namespace srt { // Samples per frame slot; the gate is compile-time so non-CP // targets keep their previous codegen exactly (the runtime form // measured +6-8% on the M55 ratchet from hot-loop branch bloat). - const std::size_t w = (k_k_channel_parallel && m_frame_major) ? m_channels : 1; + const std::size_t w = (k_channel_parallel && m_frame_major) ? m_channels : 1; for (auto& h : m_hist) { std::memmove(h.data(), h.data() + (m_end - keep) * w, keep * w * sizeof(S)); } m_end = keep; } const S* frame = m_scratch.data() + m_scratch_pos * m_channels; - if (k_k_channel_parallel && m_frame_major) { // frames stay interleaved: one contiguous copy + if (k_channel_parallel && m_frame_major) { // frames stay interleaved: one contiguous copy std::memcpy(m_hist[0].data() + m_end * m_channels, frame, m_channels * sizeof(S)); } else { diff --git a/include/srt/sample_traits.h b/include/srt/sample_traits.h index a43f2ed..0347d0e 100644 --- a/include/srt/sample_traits.h +++ b/include/srt/sample_traits.h @@ -75,6 +75,9 @@ namespace srt { /// Convert a double-precision designed coefficient to storage form. static coeff make_coeff(double c) noexcept { return static_cast(c); } + /// Coefficient units per 1.0 (used by the bank's row-sum-preserving + /// quantization; unity for floating storage, where no correction runs). + static constexpr double k_coeff_scale = 1.0; /// Convert the intra-phase fraction (in [0,1)) once per output sample. static blend_factor make_blend_factor(double fr) noexcept { return static_cast(fr); } @@ -127,6 +130,7 @@ namespace srt { static coeff make_coeff(double c) noexcept { return detail::round_sat(c * 16384.0); // Q1.14 } + static constexpr double k_coeff_scale = 16384.0; // Q1.14 units per 1.0 // ANCHOR_END: st_q15_coeff static blend_factor make_blend_factor(double fr) noexcept { @@ -189,6 +193,7 @@ namespace srt { static coeff make_coeff(double c) noexcept { return detail::round_sat(c * 1073741824.0); // Q1.30 } + static constexpr double k_coeff_scale = 1073741824.0; // Q1.30 units per 1.0 static blend_factor make_blend_factor(double fr) noexcept { return static_cast(fr * 1048576.0); // Q20 diff --git a/include/srt/spsc_ring.h b/include/srt/spsc_ring.h index 5d971ba..d29b423 100644 --- a/include/srt/spsc_ring.h +++ b/include/srt/spsc_ring.h @@ -126,14 +126,14 @@ namespace srt { // deliberately avoided: it is ABI-fragile and warns on GCC). The // read-only members (buf_, mask_) lead so they share a line with each // other, not with either side's mutable state. - static constexpr std::size_t k_k_cache_line = 64; + static constexpr std::size_t k_cache_line = 64; std::vector m_buf; std::size_t m_mask; - alignas(k_k_cache_line) std::atomic m_head{0}; // written by producer - alignas(k_k_cache_line) std::size_t m_tail_cache{0}; // producer's view of tail - alignas(k_k_cache_line) std::atomic m_tail{0}; // written by consumer - alignas(k_k_cache_line) std::size_t m_head_cache{0}; // consumer's view of head + alignas(k_cache_line) std::atomic m_head{0}; // written by producer + alignas(k_cache_line) std::size_t m_tail_cache{0}; // producer's view of tail + alignas(k_cache_line) std::atomic m_tail{0}; // written by consumer + alignas(k_cache_line) std::size_t m_head_cache{0}; // consumer's view of head // ANCHOR_END: layout }; diff --git a/tests/test_asrc_lock.cpp b/tests/test_asrc_lock.cpp index 6619e06..ccfbca5 100644 --- a/tests/test_asrc_lock.cpp +++ b/tests/test_asrc_lock.cpp @@ -9,7 +9,7 @@ namespace { - constexpr double k_k_fs = 48000.0; + constexpr double k_fs = 48000.0; srt::config mono_config() { srt::config cfg; @@ -19,7 +19,7 @@ namespace { TEST(AsrcLock, LocksAndHoldsAtConstantOffset) { srt::async_sample_rate_converter asrc(mono_config()); - srt_test::two_clock_sim sim{.asrc = asrc, .fs_in = k_k_fs * (1.0 + 200e-6), .fs_out = k_k_fs, .channels = 1}; + srt_test::two_clock_sim sim{.asrc = asrc, .fs_in = k_fs * (1.0 + 200e-6), .fs_out = k_fs, .channels = 1}; bool locked_by2s = false; double ppm_sum = 0.0; double fill_sum = 0.0; @@ -53,7 +53,7 @@ namespace { TEST(AsrcLock, TracksDriftRampWithoutUnlocking) { srt::async_sample_rate_converter asrc(mono_config()); srt_test::two_clock_sim sim{ - .asrc = asrc, .fs_in = k_k_fs, .fs_out = k_k_fs, .channels = 1, .chunk_in = 1, .chunk_out = 1}; + .asrc = asrc, .fs_in = k_fs, .fs_out = k_fs, .channels = 1, .chunk_in = 1, .chunk_out = 1}; // Input clock drifts 0 -> +300 ppm over 30 s (10 ppm/s, far faster than // real oscillator wander), then holds for the loop to reconverge. sim.fs_in_scale = [](double t) { return 1.0 + 300e-6 * std::min(t, 30.0) / 30.0; }; @@ -81,15 +81,11 @@ namespace { // sine's second difference is bounded by A*omega^2; any window-shift // discontinuity would blow far past that bound. srt::async_sample_rate_converter asrc(mono_config()); - srt_test::two_clock_sim sim{.asrc = asrc, - .fs_in = k_k_fs * (1.0 + 500e-6), - .fs_out = k_k_fs, - .channels = 1, - .chunk_in = 1, - .chunk_out = 1}; - const double amp = 0.5; - const double nu = 1000.0 / k_k_fs; - sim.gen = [&](std::uint64_t i) { + srt_test::two_clock_sim sim{ + .asrc = asrc, .fs_in = k_fs * (1.0 + 500e-6), .fs_out = k_fs, .channels = 1, .chunk_in = 1, .chunk_out = 1}; + const double amp = 0.5; + const double nu = 1000.0 / k_fs; + sim.gen = [&](std::uint64_t i) { return static_cast(amp * std::sin(2.0 * std::numbers::pi * nu * static_cast(i))); }; std::vector tail; @@ -115,7 +111,7 @@ namespace { // Producer keeps pushing while the consumer stops pulling: occupancy blows // through the high watermark, the converter hard-resyncs, then relocks. srt::async_sample_rate_converter asrc(mono_config()); - srt_test::two_clock_sim sim{.asrc = asrc, .fs_in = k_k_fs * (1.0 + 100e-6), .fs_out = k_k_fs, .channels = 1}; + srt_test::two_clock_sim sim{.asrc = asrc, .fs_in = k_fs * (1.0 + 100e-6), .fs_out = k_fs, .channels = 1}; sim.run(10.0, [&](const float*, std::size_t, double) {}); ASSERT_EQ(asrc.status().state, srt::converter_state::locked); @@ -128,7 +124,7 @@ namespace { // Resume pulling; the converter must resync and relock without underruns // turning permanent. - srt_test::two_clock_sim resume{.asrc = asrc, .fs_in = k_k_fs * (1.0 + 100e-6), .fs_out = k_k_fs, .channels = 1}; + srt_test::two_clock_sim resume{.asrc = asrc, .fs_in = k_fs * (1.0 + 100e-6), .fs_out = k_fs, .channels = 1}; resume.run(10.0, [&](const float*, std::size_t, double) {}); const auto st = asrc.status(); EXPECT_EQ(st.state, srt::converter_state::locked); diff --git a/tests/test_asrc_program.cpp b/tests/test_asrc_program.cpp index 0d9f685..1de38ea 100644 --- a/tests/test_asrc_program.cpp +++ b/tests/test_asrc_program.cpp @@ -16,8 +16,8 @@ namespace { - constexpr double k_k_fs = 48000.0; - constexpr double k_k_eps = 200e-6; + constexpr double k_fs = 48000.0; + constexpr double k_eps = 200e-6; // ANCHOR: pw_measure // 24 pink-weighted tones, 60 Hz - 16 kHz, through a +200 ppm offset; the @@ -28,9 +28,9 @@ namespace { cfg.channels = 1; cfg.filter = spec; srt::async_sample_rate_converter asrc(cfg); - const double fs_in = k_k_fs * (1.0 + k_k_eps); + const double fs_in = k_fs * (1.0 + k_eps); srt_test::two_clock_sim sim{ - .asrc = asrc, .fs_in = fs_in, .fs_out = k_k_fs, .channels = 1, .chunk_in = 1, .chunk_out = 1}; + .asrc = asrc, .fs_in = fs_in, .fs_out = k_fs, .channels = 1, .chunk_in = 1, .chunk_out = 1}; const auto comb = srt_test::tone_comb::pink(24, 60.0, 16000.0, 0.9); sim.gen = [&](std::uint64_t i) { return static_cast(comb.sample_at(i, fs_in)); }; std::vector tail; @@ -43,7 +43,7 @@ namespace { }); EXPECT_EQ(asrc.status().underruns, 0u); EXPECT_EQ(asrc.status().state, srt::converter_state::locked); - const double snr = srt_test::program_weighted_snr_db(tail, comb, fs_in, k_k_fs); + const double snr = srt_test::program_weighted_snr_db(tail, comb, fs_in, k_fs); std::printf("[ measured ] program-weighted (24 pink tones), %zu phases x %zu taps: %.1f dB\n", spec.num_phases, spec.taps_per_phase, snr); return snr; @@ -57,14 +57,10 @@ namespace { cfg.channels = 1; cfg.filter = spec; srt::async_sample_rate_converter asrc(cfg); - srt_test::two_clock_sim sim{.asrc = asrc, - .fs_in = k_k_fs * (1.0 + k_k_eps), - .fs_out = k_k_fs, - .channels = 1, - .chunk_in = 1, - .chunk_out = 1}; - const double nu_in = freq_hz / k_k_fs; - sim.gen = [&](std::uint64_t i) { + srt_test::two_clock_sim sim{ + .asrc = asrc, .fs_in = k_fs * (1.0 + k_eps), .fs_out = k_fs, .channels = 1, .chunk_in = 1, .chunk_out = 1}; + const double nu_in = freq_hz / k_fs; + sim.gen = [&](std::uint64_t i) { return static_cast(0.5 * std::sin(2.0 * std::numbers::pi * nu_in * static_cast(i))); }; std::vector tail; @@ -74,7 +70,7 @@ namespace { tail.insert(tail.end(), x, x + frames); } }); - const auto fit = srt_test::fit_sine_tracked(tail, nu_in * (1.0 + k_k_eps)); + const auto fit = srt_test::fit_sine_tracked(tail, nu_in * (1.0 + k_eps)); const double snr = srt_test::snr_db(fit); std::printf("[ measured ] economy %5.0f Hz sine: %.1f dB\n", freq_hz, snr); return snr; @@ -91,12 +87,12 @@ namespace { double v = 0.0; for (std::size_t k = 0; k < comb.freq_hz.size(); ++k) { v += comb.amplitude[k] - * std::sin(2.0 * std::numbers::pi * comb.freq_hz[k] * rho / k_k_fs * static_cast(i) + * std::sin(2.0 * std::numbers::pi * comb.freq_hz[k] * rho / k_fs * static_cast(i) + comb.phase[k]); } tail[i] = static_cast(v); } - const double snr = srt_test::program_weighted_snr_db(tail, comb, k_k_fs * (1.0 + k_k_eps), k_k_fs); + const double snr = srt_test::program_weighted_snr_db(tail, comb, k_fs * (1.0 + k_eps), k_fs); std::printf("[ measured ] instrument floor (synthetic exact tones): %.1f dB\n", snr); // float storage of the tail quantizes at ~ -150 dB; the fit must reach // it (measured 151.9 dB). diff --git a/tests/test_asrc_quality.cpp b/tests/test_asrc_quality.cpp index 4881765..3adba2d 100644 --- a/tests/test_asrc_quality.cpp +++ b/tests/test_asrc_quality.cpp @@ -11,9 +11,9 @@ namespace { - constexpr double k_k_fs = 48000.0; - constexpr double k_k_eps = 200e-6; - constexpr double k_k_amp = 0.5; + constexpr double k_fs = 48000.0; + constexpr double k_eps = 200e-6; + constexpr double k_amp = 0.5; // Resamples a sine across a +200 ppm clock offset (sample-synchronous // transfer) and measures the residual after removing the fitted fundamental @@ -24,15 +24,11 @@ namespace { cfg.channels = 1; cfg.filter = spec; srt::async_sample_rate_converter asrc(cfg); - srt_test::two_clock_sim sim{.asrc = asrc, - .fs_in = k_k_fs * (1.0 + k_k_eps), - .fs_out = k_k_fs, - .channels = 1, - .chunk_in = 1, - .chunk_out = 1}; - const double nu_in = freq_hz / k_k_fs; - sim.gen = [&](std::uint64_t i) { - return static_cast(k_k_amp * std::sin(2.0 * std::numbers::pi * nu_in * static_cast(i))); + srt_test::two_clock_sim sim{ + .asrc = asrc, .fs_in = k_fs * (1.0 + k_eps), .fs_out = k_fs, .channels = 1, .chunk_in = 1, .chunk_out = 1}; + const double nu_in = freq_hz / k_fs; + sim.gen = [&](std::uint64_t i) { + return static_cast(k_amp * std::sin(2.0 * std::numbers::pi * nu_in * static_cast(i))); }; std::vector tail; tail.reserve(48000); @@ -46,9 +42,9 @@ namespace { }); EXPECT_EQ(asrc.status().underruns, 0u); EXPECT_EQ(asrc.status().state, srt::converter_state::locked); - const double nu_out_expected = nu_in * (1.0 + k_k_eps); + const double nu_out_expected = nu_in * (1.0 + k_eps); const auto fit = srt_test::fit_sine_tracked(tail, nu_out_expected); - EXPECT_NEAR(fit.amplitude, k_k_amp, 0.01); + EXPECT_NEAR(fit.amplitude, k_amp, 0.01); // The tracked frequency must still match the true clock ratio closely. EXPECT_NEAR(fit.freq_norm / nu_out_expected, 1.0, 2e-6); const double snr = srt_test::snr_db(fit); diff --git a/tests/test_asrc_quality_16k.cpp b/tests/test_asrc_quality_16k.cpp index 55331f5..8e49430 100644 --- a/tests/test_asrc_quality_16k.cpp +++ b/tests/test_asrc_quality_16k.cpp @@ -30,9 +30,9 @@ namespace { - constexpr double k_k_fs = 16000.0; - constexpr double k_k_eps = 200e-6; - constexpr double k_k_amp = 0.5; + constexpr double k_fs = 16000.0; + constexpr double k_eps = 200e-6; + constexpr double k_amp = 0.5; // Resamples a sine across a +200 ppm clock offset (sample-synchronous // transfer) and measures the residual after removing the fitted fundamental @@ -41,18 +41,14 @@ namespace { // from Config::forSampleRate (filter band edges, servo bandwidths and // hold times). double measure_snr_db16k(double freq_hz) { - srt::config cfg = srt::config::for_sample_rate(k_k_fs); + srt::config cfg = srt::config::for_sample_rate(k_fs); cfg.channels = 1; srt::async_sample_rate_converter asrc(cfg); - srt_test::two_clock_sim sim{.asrc = asrc, - .fs_in = k_k_fs * (1.0 + k_k_eps), - .fs_out = k_k_fs, - .channels = 1, - .chunk_in = 1, - .chunk_out = 1}; - const double nu_in = freq_hz / k_k_fs; - sim.gen = [&](std::uint64_t i) { - return static_cast(k_k_amp * std::sin(2.0 * std::numbers::pi * nu_in * static_cast(i))); + srt_test::two_clock_sim sim{ + .asrc = asrc, .fs_in = k_fs * (1.0 + k_eps), .fs_out = k_fs, .channels = 1, .chunk_in = 1, .chunk_out = 1}; + const double nu_in = freq_hz / k_fs; + sim.gen = [&](std::uint64_t i) { + return static_cast(k_amp * std::sin(2.0 * std::numbers::pi * nu_in * static_cast(i))); }; std::vector tail; tail.reserve(16000); @@ -69,9 +65,9 @@ namespace { }); EXPECT_EQ(asrc.status().underruns, 0u); EXPECT_EQ(asrc.status().state, srt::converter_state::locked); - const double nu_out_expected = nu_in * (1.0 + k_k_eps); + const double nu_out_expected = nu_in * (1.0 + k_eps); const auto fit = srt_test::fit_sine_tracked(tail, nu_out_expected); - EXPECT_NEAR(fit.amplitude, k_k_amp, 0.01); + EXPECT_NEAR(fit.amplitude, k_amp, 0.01); // The tracked frequency must still match the true clock ratio closely. EXPECT_NEAR(fit.freq_norm / nu_out_expected, 1.0, 2e-6); const double snr = srt_test::snr_db(fit); diff --git a/tests/test_fixed_point.cpp b/tests/test_fixed_point.cpp index 061ba4f..dc54971 100644 --- a/tests/test_fixed_point.cpp +++ b/tests/test_fixed_point.cpp @@ -12,8 +12,8 @@ namespace { - constexpr double k_k_fs = 48000.0; - constexpr double k_k_eps = 200e-6; + constexpr double k_fs = 48000.0; + constexpr double k_eps = 200e-6; using q15 = srt::sample_traits; using q31 = srt::sample_traits; @@ -36,28 +36,50 @@ namespace { EXPECT_EQ(q31::finalize(-(std::int64_t{1} << 60)), -2147483648LL); } + // "Note that for DC, this should get you infinite S/N ratio... for every + // phase or fractional delay, the FIR coefficients must add to 1." + // -- R. Bristow-Johnson, music-dsp. With row-sum-preserving quantization + // the property survives fixed point exactly: measured 0 LSB deviation over + // a 256-point mu sweep for both formats (tolerance 1 for safety only). TEST(FixedPoint, DcGainIsUnityQ15) { - const srt::polyphase_filter_bank bank(srt::filter_spec::balanced(), k_k_fs); + const srt::polyphase_filter_bank bank(srt::filter_spec::balanced(), k_fs); std::vector dc(bank.taps(), 32767); for (int i = 0; i < 16; ++i) { const double mu = static_cast(i) / 16.0; - // 12 LSB = a -0.003 dB branch-gain deviation: per-tap rounding of - // the compensated (rect-smoothed) rows biases a few LSB further - // than the plain Kaiser rows did — still 3x inside the +/-0.01 dB - // passband contract (37 LSB at this scale). - EXPECT_NEAR(srt::interpolate(bank, dc.data(), mu), 32767, 12) << "mu=" << mu; + EXPECT_NEAR(srt::interpolate(bank, dc.data(), mu), 32767, 1) << "mu=" << mu; } } TEST(FixedPoint, DcGainIsUnityQ31) { - const srt::polyphase_filter_bank bank(srt::filter_spec::balanced(), k_k_fs); + const srt::polyphase_filter_bank bank(srt::filter_spec::balanced(), k_fs); std::vector dc(bank.taps(), 2147483647); for (int i = 0; i < 16; ++i) { const double mu = static_cast(i) / 16.0; - EXPECT_NEAR(srt::interpolate(bank, dc.data(), mu), 2147483647.0, 256.0) << "mu=" << mu; + EXPECT_NEAR(srt::interpolate(bank, dc.data(), mu), 2147483647.0, 1.0) << "mu=" << mu; } } + // The mechanism behind the previous two tests: every quantized row sums to + // exactly k_coeff_scale (the largest-remainder correction in the bank ctor). + template + void check_row_sums_exact() { + const srt::polyphase_filter_bank bank(srt::filter_spec::balanced(), k_fs); + const auto scale = static_cast(srt::sample_traits::k_coeff_scale); + for (std::size_t p = 0; p < bank.num_phases(); ++p) { + std::int64_t sum = 0; + for (std::size_t t = 0; t < bank.taps(); ++t) { + sum += bank.phase(p)[t]; + } + ASSERT_EQ(sum, scale) << "row " << p; + } + } + TEST(FixedPoint, RowSumsAreExactQ15) { + check_row_sums_exact(); + } + TEST(FixedPoint, RowSumsAreExactQ31) { + check_row_sums_exact(); + } + // End-to-end quality across a +200 ppm clock crossing, like the float suite: // resample a sine, fit and remove the fundamental, measure the residual. template @@ -65,15 +87,11 @@ namespace { srt::config cfg; cfg.channels = 1; srt::basic_async_sample_rate_converter asrc(cfg); - srt_test::two_clock_sim_t sim{.asrc = asrc, - .fs_in = k_k_fs * (1.0 + k_k_eps), - .fs_out = k_k_fs, - .channels = 1, - .chunk_in = 1, - .chunk_out = 1}; - const double nu_in = freq_hz / k_k_fs; - const double full_scale = static_cast(std::numeric_limits::max()); - sim.gen = [&](std::uint64_t i) { + srt_test::two_clock_sim_t sim{ + .asrc = asrc, .fs_in = k_fs * (1.0 + k_eps), .fs_out = k_fs, .channels = 1, .chunk_in = 1, .chunk_out = 1}; + const double nu_in = freq_hz / k_fs; + const double full_scale = static_cast(std::numeric_limits::max()); + sim.gen = [&](std::uint64_t i) { const double v = amp * std::sin(2.0 * std::numbers::pi * nu_in * static_cast(i)); return srt::detail::round_sat(v * full_scale); }; @@ -89,7 +107,7 @@ namespace { }); EXPECT_EQ(asrc.status().underruns, 0u); EXPECT_EQ(asrc.status().state, srt::converter_state::locked); - const auto fit = srt_test::fit_sine_tracked(tail, nu_in * (1.0 + k_k_eps)); + const auto fit = srt_test::fit_sine_tracked(tail, nu_in * (1.0 + k_eps)); EXPECT_NEAR(fit.amplitude, amp, 0.01); const double snr = srt_test::snr_db(fit); std::printf("[ measured ] %5.0f Hz, %d-bit fixed: SNR %.1f dB\n", freq_hz, int(sizeof(S) * 8), snr); @@ -118,14 +136,10 @@ namespace { srt::config cfg; cfg.channels = 1; srt::async_sample_rate_converter_q15 asrc(cfg); - srt_test::two_clock_sim_t sim{.asrc = asrc, - .fs_in = k_k_fs * (1.0 + 500e-6), - .fs_out = k_k_fs, - .channels = 1, - .chunk_in = 1, - .chunk_out = 1}; - const double nu = 1000.0 / k_k_fs; - sim.gen = [&](std::uint64_t i) { + srt_test::two_clock_sim_t sim{ + .asrc = asrc, .fs_in = k_fs * (1.0 + 500e-6), .fs_out = k_fs, .channels = 1, .chunk_in = 1, .chunk_out = 1}; + const double nu = 1000.0 / k_fs; + sim.gen = [&](std::uint64_t i) { return srt::detail::round_sat( 0.99 * 32767.0 * std::sin(2.0 * std::numbers::pi * nu * static_cast(i))); }; diff --git a/tests/test_hardening.cpp b/tests/test_hardening.cpp index f37ac83..9ffc581 100644 --- a/tests/test_hardening.cpp +++ b/tests/test_hardening.cpp @@ -17,7 +17,7 @@ namespace { - constexpr double k_k_fs = 48000.0; + constexpr double k_fs = 48000.0; // Audit finding F1: with defaults, any pull block larger than the 48-frame // setpoint used to drain into a permanent underrun limit cycle (64-frame @@ -39,8 +39,8 @@ namespace { } srt::async_sample_rate_converter asrc(cfg); srt_test::two_clock_sim sim{.asrc = asrc, - .fs_in = k_k_fs * (1.0 + 200e-6), - .fs_out = k_k_fs, + .fs_in = k_fs * (1.0 + 200e-6), + .fs_out = k_fs, .channels = 1, .chunk_in = 32, .chunk_out = pull_block}; @@ -76,7 +76,7 @@ namespace { srt::config cfg; cfg.channels = 1; srt::async_sample_rate_converter asrc(cfg); - srt_test::two_clock_sim sim{.asrc = asrc, .fs_in = k_k_fs * (1.0 + 200e-6), .fs_out = k_k_fs, .channels = 1}; + srt_test::two_clock_sim sim{.asrc = asrc, .fs_in = k_fs * (1.0 + 200e-6), .fs_out = k_fs, .channels = 1}; sim.run(5.0, [](const float*, std::size_t, double) {}); // 32-frame pulls against the 48-frame default were always feasible; // the adaptation must not inflate latency for them. @@ -148,12 +148,12 @@ namespace { srt::config cfg; cfg.channels = 1; srt::async_sample_rate_converter asrc(cfg); - srt_test::two_clock_sim sim{.asrc = asrc, .fs_in = k_k_fs * (1.0 + 200e-6), .fs_out = k_k_fs, .channels = 1}; + srt_test::two_clock_sim sim{.asrc = asrc, .fs_in = k_fs * (1.0 + 200e-6), .fs_out = k_fs, .channels = 1}; sim.run(5.0, [](const float*, std::size_t, double) {}); ASSERT_EQ(asrc.status().state, srt::converter_state::locked); asrc.reset_from_consumer(); EXPECT_EQ(asrc.status().state, srt::converter_state::filling); - srt_test::two_clock_sim sim2{.asrc = asrc, .fs_in = k_k_fs * (1.0 + 200e-6), .fs_out = k_k_fs, .channels = 1}; + srt_test::two_clock_sim sim2{.asrc = asrc, .fs_in = k_fs * (1.0 + 200e-6), .fs_out = k_fs, .channels = 1}; sim2.run(5.0, [](const float*, std::size_t, double) {}); EXPECT_EQ(asrc.status().state, srt::converter_state::locked); } @@ -209,14 +209,10 @@ namespace { srt::config cfg; cfg.channels = 1; srt::async_sample_rate_converter_q15 asrc(cfg); - srt_test::two_clock_sim_t sim{.asrc = asrc, - .fs_in = k_k_fs * (1.0 + 200e-6), - .fs_out = k_k_fs, - .channels = 1, - .chunk_in = 8, - .chunk_out = 8}; - const double nu = 997.0 / k_k_fs; - sim.gen = [&](std::uint64_t i) { + srt_test::two_clock_sim_t sim{ + .asrc = asrc, .fs_in = k_fs * (1.0 + 200e-6), .fs_out = k_fs, .channels = 1, .chunk_in = 8, .chunk_out = 8}; + const double nu = 997.0 / k_fs; + sim.gen = [&](std::uint64_t i) { return srt::detail::round_sat( 0.5 * 32767.0 * std::sin(2.0 * std::numbers::pi * nu * static_cast(i))); }; @@ -244,14 +240,10 @@ namespace { srt::config cfg; cfg.channels = 1; srt::async_sample_rate_converter_q15 asrc(cfg); - srt_test::two_clock_sim_t sim{.asrc = asrc, - .fs_in = k_k_fs * (1.0 + 500e-6), - .fs_out = k_k_fs, - .channels = 1, - .chunk_in = 8, - .chunk_out = 8}; - const double nu = 1000.0 / k_k_fs; - sim.gen = [&](std::uint64_t i) { + srt_test::two_clock_sim_t sim{ + .asrc = asrc, .fs_in = k_fs * (1.0 + 500e-6), .fs_out = k_fs, .channels = 1, .chunk_in = 8, .chunk_out = 8}; + const double nu = 1000.0 / k_fs; + sim.gen = [&](std::uint64_t i) { return srt::detail::round_sat( 0.99 * 32767.0 * std::sin(2.0 * std::numbers::pi * nu * static_cast(i))); }; diff --git a/tests/test_kaiser.cpp b/tests/test_kaiser.cpp index 348f5be..ceb9d25 100644 --- a/tests/test_kaiser.cpp +++ b/tests/test_kaiser.cpp @@ -99,6 +99,30 @@ namespace { // The compensated presets must also hold their specs at scaled rates (the // 16 kHz deployment path): normalized design, same numbers. + // The k*fs transmission zeros ARE branch-DC uniformity, stated in the + // frequency domain: with exact zeros, every polyphase branch's coefficient + // sum is identical (measured spread 1.8e-15 -- machine epsilon -- vs 4.7e-6 + // for the plain fast() design, whose spread is its stopband leakage at fs). + TEST(Kaiser, CompensatedBranchSumsAreUniform) { + const auto spec = srt::filter_spec::balanced(); + const std::size_t phases = std::bit_ceil(spec.num_phases); + std::vector h(phases * spec.taps_per_phase); + design_prototype_compensated(h, phases, (spec.passband_hz + spec.stopband_hz) / 48000.0, + kaiser_beta(spec.stopband_atten_db), spec.passband_hz / 48000.0); + double lo = 1e9; + double hi = -1e9; + for (std::size_t p = 0; p < phases; ++p) { + double sum = 0.0; + for (std::size_t t = 0; t < spec.taps_per_phase; ++t) { + sum += h[t * phases + p]; + } + lo = std::min(lo, sum); + hi = std::max(hi, sum); + } + EXPECT_LT(hi - lo, 1e-12); + EXPECT_NEAR(lo, 1.0, 1e-9); + } + TEST(Kaiser, CompensatedSpecsHoldAt16k) { check_prototype_meets_spec(srt::filter_spec::balanced().scaled_to(16000.0), 16000.0); check_prototype_meets_spec(srt::filter_spec::economy().scaled_to(16000.0), 16000.0); diff --git a/tests/test_latency.cpp b/tests/test_latency.cpp index 0de35d3..d5a5460 100644 --- a/tests/test_latency.cpp +++ b/tests/test_latency.cpp @@ -8,14 +8,14 @@ namespace { - constexpr double k_k_fs = 48000.0; + constexpr double k_fs = 48000.0; TEST(Latency, ImpulseDelayMatchesDesignedLatency) { srt::config cfg; cfg.channels = 1; srt::async_sample_rate_converter asrc(cfg); srt_test::two_clock_sim sim{ - .asrc = asrc, .fs_in = k_k_fs, .fs_out = k_k_fs, .channels = 1, .chunk_in = 1, .chunk_out = 1}; + .asrc = asrc, .fs_in = k_fs, .fs_out = k_fs, .channels = 1, .chunk_in = 1, .chunk_out = 1}; const std::uint64_t impulse_index = 24000; // 0.5 s in, well past acquisition sim.gen = [&](std::uint64_t i) { return i == impulse_index ? 1.0f : 0.0f; }; std::vector out; @@ -36,7 +36,7 @@ namespace { const double refine = 0.5 * (y0 - y2) / (y0 - 2.0 * y1 + y2); const double measured_delay = static_cast(peak) + refine - static_cast(impulse_index); - const double designed_delay = asrc.designed_latency_seconds() * k_k_fs; + const double designed_delay = asrc.designed_latency_seconds() * k_fs; // The sim's event interleaving contributes up to ~1 sample of offset on // top of the designed figure. EXPECT_NEAR(measured_delay, designed_delay, 2.5); @@ -49,7 +49,7 @@ namespace { const double group_delay = asrc.filter_bank().group_delay_samples(); EXPECT_NEAR(group_delay, 24.0, 0.1); // ~T/2 for balanced (T = 48) EXPECT_NEAR(asrc.designed_latency_seconds(), - (static_cast(cfg.target_latency_frames) + group_delay) / k_k_fs, 1e-12); + (static_cast(cfg.target_latency_frames) + group_delay) / k_fs, 1e-12); // The whitepaper budget: ~1 ms core latency for the default config. EXPECT_LT(asrc.designed_latency_seconds(), 2e-3); } diff --git a/tests/test_multichannel.cpp b/tests/test_multichannel.cpp index e2bedf0..5ef5fae 100644 --- a/tests/test_multichannel.cpp +++ b/tests/test_multichannel.cpp @@ -25,9 +25,9 @@ namespace { - constexpr double k_k_fs = 48000.0; - constexpr double k_k_eps = 200e-6; - constexpr double k_k_amp = 0.4; + constexpr double k_fs = 48000.0; + constexpr double k_eps = 200e-6; + constexpr double k_amp = 0.4; /// Distinct, non-harmonically-related tone per channel, all inside the flat /// passband for up to 16 channels (max 600 + 731*15 = 11565 Hz). @@ -73,19 +73,19 @@ namespace { cfg.channels = channels; srt::basic_async_sample_rate_converter asrc(cfg); srt_test::two_clock_sim_t sim{.asrc = asrc, - .fs_in = k_k_fs * (1.0 + k_k_eps), - .fs_out = k_k_fs, + .fs_in = k_fs * (1.0 + k_eps), + .fs_out = k_fs, .channels = channels, .chunk_in = chunk, .chunk_out = chunk}; sim.gen_ch = [&](std::uint64_t i, std::size_t c) { - const double w = 2.0 * std::numbers::pi * channel_freq_hz(c) / k_k_fs; + const double w = 2.0 * std::numbers::pi * channel_freq_hz(c) / k_fs; // Per-channel phase offsets decorrelate the channel waveforms. - return make_sample(k_k_amp * std::sin(w * static_cast(i) + 0.7 * static_cast(c))); + return make_sample(k_amp * std::sin(w * static_cast(i) + 0.7 * static_cast(c))); }; std::vector tail; - tail.reserve(static_cast(window_seconds * k_k_fs + 16.0) * channels); + tail.reserve(static_cast(window_seconds * k_fs + 16.0) * channels); sim.run(total_seconds, [&](const S* x, std::size_t frames, double t) { if (t >= total_seconds - window_seconds) { tail.insert(tail.end(), x, x + frames * channels); @@ -103,7 +103,7 @@ namespace { } // Own tone: tracked fit, then exact removal of the fitted component. - const double nu_own = channel_freq_hz(c) / k_k_fs * (1.0 + k_k_eps); + const double nu_own = channel_freq_hz(c) / k_fs * (1.0 + k_eps); const auto own = srt_test::fit_sine_tracked(x, nu_own); reports[c].amplitude = own.amplitude; reports[c].snr_db = srt_test::snr_db(own); @@ -119,7 +119,7 @@ namespace { if (k == c) { continue; } - const double nu_k = channel_freq_hz(k) / k_k_fs * (1.0 + k_k_eps); + const double nu_k = channel_freq_hz(k) / k_fs * (1.0 + k_eps); const auto leak = srt_test::fit_sine(x, nu_k); const double db = 20.0 * std::log10(leak.amplitude / own.amplitude); if (db > reports[c].worst_crosstalk_db) { @@ -139,7 +139,7 @@ namespace { TEST(MultiChannel, Independence12chFloat) { const auto r = measure_independence(12, 40.0, 1.0, 1); for (const auto& ch : r) { - EXPECT_NEAR(ch.amplitude, k_k_amp, 0.01); + EXPECT_NEAR(ch.amplitude, k_amp, 0.01); EXPECT_GT(ch.snr_db, 100.0); EXPECT_LT(ch.worst_crosstalk_db, -100.0); } @@ -148,7 +148,7 @@ namespace { TEST(MultiChannel, Independence16chQ15) { const auto r = measure_independence(16, 40.0, 1.0, 1); for (const auto& ch : r) { - EXPECT_NEAR(ch.amplitude, k_k_amp, 0.01); + EXPECT_NEAR(ch.amplitude, k_amp, 0.01); EXPECT_GT(ch.snr_db, 72.0); EXPECT_LT(ch.worst_crosstalk_db, -72.0); } @@ -165,7 +165,7 @@ namespace { TEST(MultiChannelShort, Independence5chFloat) { const auto r = measure_independence(5, 4.0, 0.25, 8); for (const auto& ch : r) { - EXPECT_NEAR(ch.amplitude, k_k_amp, 0.05); + EXPECT_NEAR(ch.amplitude, k_amp, 0.05); EXPECT_GT(ch.snr_db, 35.0); EXPECT_LT(ch.worst_crosstalk_db, -50.0); } @@ -174,7 +174,7 @@ namespace { TEST(MultiChannelShort, Independence7chFloat) { const auto r = measure_independence(7, 4.0, 0.25, 8); for (const auto& ch : r) { - EXPECT_NEAR(ch.amplitude, k_k_amp, 0.05); + EXPECT_NEAR(ch.amplitude, k_amp, 0.05); EXPECT_GT(ch.snr_db, 35.0); EXPECT_LT(ch.worst_crosstalk_db, -50.0); } @@ -183,7 +183,7 @@ namespace { TEST(MultiChannelShort, Independence12chQ15) { const auto r = measure_independence(12, 4.0, 0.25, 8); for (const auto& ch : r) { - EXPECT_NEAR(ch.amplitude, k_k_amp, 0.05); + EXPECT_NEAR(ch.amplitude, k_amp, 0.05); EXPECT_GT(ch.snr_db, 35.0); EXPECT_LT(ch.worst_crosstalk_db, -45.0); } diff --git a/tests/test_polyphase.cpp b/tests/test_polyphase.cpp index 1c417af..a864230 100644 --- a/tests/test_polyphase.cpp +++ b/tests/test_polyphase.cpp @@ -9,10 +9,10 @@ namespace { - constexpr double k_k_fs = 48000.0; + constexpr double k_fs = 48000.0; TEST(Polyphase, DcGainIsUnityAcrossMu) { - const srt::polyphase_filter_bank bank(srt::filter_spec::balanced(), k_k_fs); + const srt::polyphase_filter_bank bank(srt::filter_spec::balanced(), k_fs); std::vector ones(bank.taps(), 1.0f); std::mt19937 rng(7); std::uniform_real_distribution uni(0.0, 1.0); @@ -23,7 +23,7 @@ namespace { } TEST(Polyphase, ExtraRowEqualsPhaseZeroAdvancedOneTap) { - const srt::polyphase_filter_bank bank(srt::filter_spec::balanced(), k_k_fs); + const srt::polyphase_filter_bank bank(srt::filter_spec::balanced(), k_fs); const std::size_t L = bank.num_phases(); const std::size_t T = bank.taps(); // Rows are stored tap-reversed over an oldest-first window, so "advanced @@ -39,7 +39,7 @@ namespace { // The interpolated output at fractional position mu corresponds to input time // tau = J - T/2 + mu + 1/(2L) where J is the newest sample index in the window. double max_error_db(const srt::polyphase_filter_bank& bank, double freq_hz) { - const double nu = freq_hz / k_k_fs; + const double nu = freq_hz / k_fs; const std::size_t T = bank.taps(); const double L = static_cast(bank.num_phases()); std::vector x(4 * T); @@ -61,7 +61,7 @@ namespace { } TEST(Polyphase, FractionalDelayAccuracyBalanced) { - const srt::polyphase_filter_bank bank(srt::filter_spec::balanced(), k_k_fs); + const srt::polyphase_filter_bank bank(srt::filter_spec::balanced(), k_fs); // Error budget: this absolute-error sweep sees BOTH the inter-phase // interpolation floor and the prototype's in-spec passband ripple (a // gain deviation of r dB reads here as 20*log10(10^(r/20)-1): the @@ -80,7 +80,7 @@ namespace { } TEST(Polyphase, FractionalDelayAccuracyTransparent) { - const srt::polyphase_filter_bank bank(srt::filter_spec::transparent(), k_k_fs); + const srt::polyphase_filter_bank bank(srt::filter_spec::transparent(), k_fs); EXPECT_LT(max_error_db(bank, 997.0), -104.0); EXPECT_LT(max_error_db(bank, 19000.0), -80.0); } @@ -88,7 +88,7 @@ namespace { TEST(Polyphase, MuWrapIsContinuousWithWindowShift) { // interpolate(hist, mu -> 1) must equal interpolate(hist shifted by one // newer sample, mu = 0): the whole-sample slip invariant. - const srt::polyphase_filter_bank bank(srt::filter_spec::balanced(), k_k_fs); + const srt::polyphase_filter_bank bank(srt::filter_spec::balanced(), k_fs); const std::size_t T = bank.taps(); std::vector x(2 * T); std::mt19937 rng(99); diff --git a/tests/test_servo.cpp b/tests/test_servo.cpp index 387901d..8834495 100644 --- a/tests/test_servo.cpp +++ b/tests/test_servo.cpp @@ -6,25 +6,25 @@ namespace { - constexpr double k_k_fs = 48000.0; - constexpr double k_k_target = 48.0; - constexpr std::size_t k_k_block = 32; - constexpr double k_k_dt = static_cast(k_k_block) / k_k_fs; + constexpr double k_fs = 48000.0; + constexpr double k_target = 48.0; + constexpr std::size_t k_block = 32; + constexpr double k_dt = static_cast(k_block) / k_fs; // Pure plant simulation: the FIFO integrates the rate mismatch. struct plant { - double occ = k_k_target; - void step(double eps_true, double eps_hat) { occ += (eps_true - eps_hat) * k_k_fs * k_k_dt; } + double occ = k_target; + void step(double eps_true, double eps_hat) { occ += (eps_true - eps_hat) * k_fs * k_dt; } }; TEST(Servo, LocksFromConstantOffsetAndNullsError) { - srt::pi_servo servo(srt::servo_config{}, k_k_fs, k_k_target); + srt::pi_servo servo(srt::servo_config{}, k_fs, k_target); plant plant; const double eps_true = 300e-6; bool locked_within1_5s = false; double t = 0.0; - for (; t < 30.0; t += k_k_dt) { // locked loop is 0.05 Hz: allow it to settle - const double eps = servo.update(plant.occ, 0.0, k_k_dt); + for (; t < 30.0; t += k_dt) { // locked loop is 0.05 Hz: allow it to settle + const double eps = servo.update(plant.occ, 0.0, k_dt); plant.step(eps_true, eps); if (t < 1.5 && servo.locked()) { locked_within1_5s = true; @@ -33,26 +33,26 @@ namespace { EXPECT_TRUE(locked_within1_5s); EXPECT_TRUE(servo.locked()); // Type-2 loop: constant ppm offset leaves zero standing occupancy error. - EXPECT_NEAR(plant.occ, k_k_target, 0.05); + EXPECT_NEAR(plant.occ, k_target, 0.05); EXPECT_NEAR(servo.eps_hat(), eps_true, 1e-6); // within 1 ppm } TEST(Servo, TracksSlowDriftRampWithBoundedLag) { - srt::pi_servo servo(srt::servo_config{}, k_k_fs, k_k_target); + srt::pi_servo servo(srt::servo_config{}, k_fs, k_target); plant plant; // Settle at 0 ppm first. - for (double t = 0.0; t < 5.0; t += k_k_dt) { - plant.step(0.0, servo.update(plant.occ, 0.0, k_k_dt)); + for (double t = 0.0; t < 5.0; t += k_dt) { + plant.step(0.0, servo.update(plant.occ, 0.0, k_dt)); } ASSERT_TRUE(servo.locked()); // Then ramp 1 ppm/s for 20 s (temperature-style drift). double max_err = 0.0; double eps_true = 0.0; - for (double t = 0.0; t < 20.0; t += k_k_dt) { + for (double t = 0.0; t < 20.0; t += k_dt) { eps_true = 1e-6 * t; - plant.step(eps_true, servo.update(plant.occ, 0.0, k_k_dt)); + plant.step(eps_true, servo.update(plant.occ, 0.0, k_dt)); if (t > 5.0) { - max_err = std::max(max_err, std::abs(plant.occ - k_k_target)); + max_err = std::max(max_err, std::abs(plant.occ - k_target)); } } EXPECT_TRUE(servo.locked()); @@ -63,22 +63,22 @@ namespace { } TEST(Servo, BandwidthSwitchIsTransientFree) { - srt::pi_servo servo(srt::servo_config{}, k_k_fs, k_k_target); + srt::pi_servo servo(srt::servo_config{}, k_fs, k_target); plant plant; const double eps_true = 200e-6; // Run until just locked. double t = 0.0; while (!servo.locked() && t < 5.0) { - plant.step(eps_true, servo.update(plant.occ, 0.0, k_k_dt)); - t += k_k_dt; + plant.step(eps_true, servo.update(plant.occ, 0.0, k_dt)); + t += k_dt; } ASSERT_TRUE(servo.locked()); // The narrow-bandwidth handoff keeps the integrator, so the occupancy // must not be disturbed beyond the lock threshold afterwards. double max_err = 0.0; - for (double s = 0.0; s < 10.0; s += k_k_dt) { - plant.step(eps_true, servo.update(plant.occ, 0.0, k_k_dt)); - max_err = std::max(max_err, std::abs(plant.occ - k_k_target)); + for (double s = 0.0; s < 10.0; s += k_dt) { + plant.step(eps_true, servo.update(plant.occ, 0.0, k_dt)); + max_err = std::max(max_err, std::abs(plant.occ - k_target)); } EXPECT_LT(max_err, srt::servo_config{}.lock_threshold_frames); EXPECT_TRUE(servo.locked()); @@ -87,18 +87,18 @@ namespace { TEST(Servo, ClampsToMaxDeviation) { srt::servo_config cfg; cfg.max_deviation_ppm = 100.0; - srt::pi_servo servo(cfg, k_k_fs, k_k_target); + srt::pi_servo servo(cfg, k_fs, k_target); // Huge occupancy error must saturate at 1.5x the configured range. - const double eps = servo.update(k_k_target + 10000.0, 0.0, k_k_dt); + const double eps = servo.update(k_target + 10000.0, 0.0, k_dt); EXPECT_LE(eps, 1.5 * 100e-6 + 1e-12); } TEST(Servo, DropoutResetKeepsPpmEstimate) { - srt::pi_servo servo(srt::servo_config{}, k_k_fs, k_k_target); + srt::pi_servo servo(srt::servo_config{}, k_fs, k_target); plant plant; const double eps_true = 250e-6; - for (double t = 0.0; t < 6.0; t += k_k_dt) { - plant.step(eps_true, servo.update(plant.occ, 0.0, k_k_dt)); + for (double t = 0.0; t < 6.0; t += k_dt) { + plant.step(eps_true, servo.update(plant.occ, 0.0, k_dt)); } ASSERT_NEAR(servo.eps_hat(), eps_true, 2e-6); servo.reset(true); // dropout: keep the integrator