From dee78bdba6f9f069ddac360f73e815d6ae8f94f6 Mon Sep 17 00:00:00 2001 From: Josh Poole Date: Mon, 14 Sep 2026 10:12:30 +0100 Subject: [PATCH 1/3] 20260914 - Pad the clutter convolution to an FFT-friendly length WienerHopf planned fftFiltX/fftFiltW/fftFilt at nBins + nSamples + 1. At the shipped config (fs 2e6, cpi 0.5, delay -10..400) that is 1000411 = 269 x 3719, both prime, so FFTW falls back to Rader instead of a fast codelet. The length only has to be at least the linear-convolution length (nSamples + nBins - 1) for the result to be alias-free, and only the first nSamples outputs are read, so it is free to round up. nextFastFftLength() picks the smallest length built from 2/3/5/7 with at most one factor of 11 or 13, giving 1002375 for the shipped config. This is structural rather than bad luck: the length depends on CPI and delay span, so filter cost jumps unpredictably whenever either is touched. Ambiguity already solves the same problem with next_hamming(); the clutter filter never got the treatment. Also fixes a latent index bug the ported test caught: `i - delayMin` promotes to unsigned, so a positive delayMin wraps at 2^32 and reads the wrong sample. Negative and zero delayMin cancel exactly, so the shipped config (-10) was unaffected. Measured on x86, production geometry, pinned to 4 cores, median of 6: clutter filter 281.12 ms -> 200.63 ms (1.40x). Output differs only by floating-point reassociation: RMS relative difference 6.05e-14, max relative error 2.95e-15, far below the 1e-6 regression tolerance. Ported from mickeyslaven/blah2-VectorWarp de45d9b (MIT). Co-Authored-By: Claude Opus 5 (1M context) --- CMakeLists.txt | 14 +++ src/process/clutter/WienerHopf.cpp | 33 ++++-- src/process/clutter/WienerHopf.h | 5 + src/process/meta/FftLength.h | 26 +++++ test/unit/process/clutter/TestClutterFft.cpp | 117 +++++++++++++++++++ 5 files changed, 183 insertions(+), 12 deletions(-) create mode 100644 src/process/meta/FftLength.h create mode 100644 test/unit/process/clutter/TestClutterFft.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 2e8cf8c4..c401c75c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -112,6 +112,20 @@ target_link_libraries(testHammingNumber PRIVATE set_target_properties(testHammingNumber PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${PROJECT_BINARY_TEST_UNIT_DIR}") +add_executable(testClutterFft + test/unit/process/clutter/TestClutterFft.cpp + src/process/clutter/WienerHopf.cpp + src/data/IqData.cpp +) +target_link_libraries(testClutterFft PRIVATE + armadillo + fftw3 + fftw3_threads +) +set_target_properties(testClutterFft PROPERTIES + RUNTIME_OUTPUT_DIRECTORY "${PROJECT_BINARY_TEST_UNIT_DIR}") + # TODO: Unsure if will be using CTest. add_test(NAME testAmbiguity COMMAND testAmbiguity) add_test(NAME testTracker COMMAND testTracker) +add_test(NAME testClutterFft COMMAND testClutterFft) diff --git a/src/process/clutter/WienerHopf.cpp b/src/process/clutter/WienerHopf.cpp index 9fa9ddc5..647a309f 100644 --- a/src/process/clutter/WienerHopf.cpp +++ b/src/process/clutter/WienerHopf.cpp @@ -1,4 +1,5 @@ #include "WienerHopf.h" +#include "process/meta/FftLength.h" #include #include #include @@ -9,8 +10,13 @@ WienerHopf::WienerHopf(int32_t _delayMin, int32_t _delayMax, uint32_t _nSamples) // input delayMin = _delayMin; delayMax = _delayMax; - nBins = delayMax - delayMin; + const int64_t taps = int64_t(delayMax) - delayMin; + if (!_nSamples || taps <= 0 || uint64_t(taps) > _nSamples) + throw std::invalid_argument("Clutter filter needs a non-empty half-open delay range no longer than the CPI"); + nBins = static_cast(taps); nSamples = _nSamples; + // Pad only the linear convolution; keep taps and circular correlations unchanged. + nFilter = blah2::nextFastFftLength(uint64_t(nSamples) + nBins + 1); // initialise data A = arma::cx_mat(nBins, nBins); @@ -25,9 +31,9 @@ WienerHopf::WienerHopf(int32_t _delayMin, int32_t _delayMax, uint32_t _nSamples) dataOutY = new std::complex[nSamples]; dataA = new std::complex[nSamples]; dataB = new std::complex[nSamples]; - filtX = new std::complex[nBins + nSamples + 1]; - filtW = new std::complex[nBins + nSamples + 1]; - filt = new std::complex[nBins + nSamples + 1]; + filtX = new std::complex[nFilter]; + filtW = new std::complex[nFilter]; + filt = new std::complex[nFilter]; fftX = fftw_plan_dft_1d(nSamples, reinterpret_cast(dataX), reinterpret_cast(dataOutX), FFTW_FORWARD, FFTW_ESTIMATE); fftY = fftw_plan_dft_1d(nSamples, reinterpret_cast(dataY), @@ -36,11 +42,11 @@ WienerHopf::WienerHopf(int32_t _delayMin, int32_t _delayMax, uint32_t _nSamples) reinterpret_cast(dataA), FFTW_BACKWARD, FFTW_ESTIMATE); fftB = fftw_plan_dft_1d(nSamples, reinterpret_cast(dataB), reinterpret_cast(dataB), FFTW_BACKWARD, FFTW_ESTIMATE); - fftFiltX = fftw_plan_dft_1d(nBins + nSamples + 1, reinterpret_cast(filtX), + fftFiltX = fftw_plan_dft_1d(nFilter, reinterpret_cast(filtX), reinterpret_cast(filtX), FFTW_FORWARD, FFTW_ESTIMATE); - fftFiltW = fftw_plan_dft_1d(nBins + nSamples + 1, reinterpret_cast(filtW), + fftFiltW = fftw_plan_dft_1d(nFilter, reinterpret_cast(filtW), reinterpret_cast(filtW), FFTW_FORWARD, FFTW_ESTIMATE); - fftFilt = fftw_plan_dft_1d(nBins + nSamples + 1, reinterpret_cast(filt), + fftFilt = fftw_plan_dft_1d(nFilter, reinterpret_cast(filt), reinterpret_cast(filt), FFTW_BACKWARD, FFTW_ESTIMATE); } @@ -64,7 +70,10 @@ bool WienerHopf::process(IqData *x, IqData *y) // change deque to std::complex for (i = 0; i < nSamples; i++) { - dataX[i] = xData[(((i - delayMin) % nSamples) + nSamples) % nSamples]; + // Signed arithmetic: `i - delayMin` promotes to unsigned, so a positive + // delayMin wraps at 2^32 and lands on the wrong sample. + const int64_t shifted = (int64_t(i) - delayMin) % int64_t(nSamples); + dataX[i] = xData[shifted < 0 ? shifted + nSamples : shifted]; dataY[i] = yData[i]; } @@ -126,7 +135,7 @@ bool WienerHopf::process(IqData *x, IqData *y) { filtX[i] = dataX[i]; } - for (i = nSamples; i < nBins + nSamples + 1; i++) + for (i = nSamples; i < nFilter; i++) { filtX[i] = {0, 0}; } @@ -136,7 +145,7 @@ bool WienerHopf::process(IqData *x, IqData *y) { filtW[i] = w[i]; } - for (i = nBins; i < nBins + nSamples + 1; i++) + for (i = nBins; i < nFilter; i++) { filtW[i] = {0, 0}; } @@ -146,7 +155,7 @@ bool WienerHopf::process(IqData *x, IqData *y) fftw_execute(fftFiltW); // compute convolution/filter - for (i = 0; i < nBins + nSamples + 1; i++) + for (i = 0; i < nFilter; i++) { filt[i] = (filtW[i] * filtX[i]); } @@ -156,7 +165,7 @@ bool WienerHopf::process(IqData *x, IqData *y) y->clear(); for (i = 0; i < nSamples; i++) { - y->push_back(dataY[i] - (filt[i] / (double)(nBins + nSamples + 1))); + y->push_back(dataY[i] - (filt[i] / (double)nFilter)); } return true; diff --git a/src/process/clutter/WienerHopf.h b/src/process/clutter/WienerHopf.h index 5f9ddc3b..e215df84 100644 --- a/src/process/clutter/WienerHopf.h +++ b/src/process/clutter/WienerHopf.h @@ -29,6 +29,9 @@ class WienerHopf /// @brief Number of samples per CPI. uint32_t nSamples; + /// @brief Zero-padded linear-convolution FFT length; does not change taps. + uint32_t nFilter; + /// @brief True if clutter filter processing is successful. bool success; @@ -67,6 +70,8 @@ class WienerHopf /// @return The object. WienerHopf(int32_t delayMin, int32_t delayMax, uint32_t nSamples); + uint32_t filter_fft_length() const { return nFilter; } + /// @brief Destructor. /// @return Void. ~WienerHopf(); diff --git a/src/process/meta/FftLength.h b/src/process/meta/FftLength.h new file mode 100644 index 00000000..918ed760 --- /dev/null +++ b/src/process/meta/FftLength.h @@ -0,0 +1,26 @@ +#pragma once + +#include +#include +#include +#include + +namespace blah2 { +// FFTW supports these small factors efficiently. Permit at most one factor of +// 11 or 13. Enumerate bounded candidates instead of an unbounded integer scan. +inline uint32_t nextFastFftLength(uint64_t minimum) { + constexpr uint64_t limit = std::numeric_limits::max(); + if (!minimum || minimum > limit) + throw std::invalid_argument("FFT length must fit a positive FFTW int"); + uint64_t best = limit + 1; + for (uint64_t extra : {1u, 11u, 13u}) + for (uint64_t a = extra; a <= limit && a < best; a *= 2) + for (uint64_t b = a; b <= limit && b < best; b *= 3) + for (uint64_t c = b; c <= limit && c < best; c *= 5) + for (uint64_t d = c; d <= limit && d < best; d *= 7) + if (d >= minimum) best = d; + if (best > limit) + throw std::invalid_argument("No supported padded FFT length fits FFTW"); + return static_cast(best); +} +} diff --git a/test/unit/process/clutter/TestClutterFft.cpp b/test/unit/process/clutter/TestClutterFft.cpp new file mode 100644 index 00000000..171d53fc --- /dev/null +++ b/test/unit/process/clutter/TestClutterFft.cpp @@ -0,0 +1,117 @@ +#include "process/clutter/WienerHopf.h" +#include "process/meta/FftLength.h" +#include +#include +#include + +using Complex = std::complex; + +static void require(bool condition, const char* message) { + if (!condition) throw std::runtime_error(message); +} + +static bool fast(uint32_t value) { + for (uint32_t factor : {2u, 3u, 5u, 7u}) + while (value % factor == 0) value /= factor; + return value == 1 || value == 11 || value == 13; +} + +static void run(unsigned samples, int first, int last) { + const unsigned taps = last - first; + IqData reference(samples), surveillance(samples); + WienerHopf filter(first, last, samples); + require(filter.filter_fft_length() == blah2::nextFastFftLength(samples + taps + 1), + "Filter did not use the selected padded length"); + std::mt19937 rng(9211); + std::normal_distribution random; + for (int repeat = 0; repeat < 3; ++repeat) { + reference.clear(); + surveillance.clear(); + arma::cx_vec x(samples), y(samples); + for (unsigned i = 0; i < samples; ++i) + reference.push_back({random(rng), random(rng)}); + const auto originalReference = reference.get_data(); + for (unsigned i = 0; i < samples; ++i) { + const int64_t shifted = (int64_t(i) - first) % samples; + x[i] = originalReference[shifted < 0 ? shifted + samples : shifted]; + y[i] = x[i] * Complex(.7, .2) + Complex(.2 * random(rng), .2 * random(rng)); + surveillance.push_back(y[i]); + } + + // Independent circular correlations and direct linear convolution, not + // another FFT implementation. Reuse the filter on three different CPIs. + arma::cx_vec a(taps, arma::fill::zeros), b(taps, arma::fill::zeros); + for (unsigned lag = 0; lag < taps; ++lag) + for (unsigned i = 0; i < samples; ++i) { + a[lag] += std::conj(x[(i + lag) % samples]) * x[i]; + b[lag] += y[(i + lag) % samples] * std::conj(x[i]); + } + arma::cx_mat matrix = arma::toeplitz(a); + for (unsigned row = 0; row < taps; ++row) + for (unsigned col = 0; col < row; ++col) + matrix(row, col) = std::conj(matrix(row, col)); + const arma::cx_vec weights = arma::solve(matrix, b); + require(filter.process(&reference, &surveillance), "Full-rank fixture rejected"); + require(surveillance.get_length() == samples, "Output sample count changed"); + const auto filtered = surveillance.get_data(); + for (unsigned i = 0; i < samples; ++i) { + Complex expected = y[i]; + for (unsigned tap = 0; tap < taps && tap <= i; ++tap) + expected -= weights[tap] * x[i - tap]; + require(std::abs(expected - filtered[i]) < 1e-9, + "Padded clutter differs from direct convolution"); + } + require(reference.get_data() == originalReference, "Reference mutated"); + } + std::cout << "PASS samples=" << samples << " taps=" << taps + << " first=" << first << " fft=" << filter.filter_fft_length() << '\n'; +} + +int main() { + try { + // Independently scan a bounded range to check both admissibility and + // minimality, including lengths containing repeated factors of 11/13. + for (uint32_t minimum = 1; minimum <= 4096; ++minimum) { + uint32_t expected = minimum; + while (!fast(expected)) ++expected; + require(blah2::nextFastFftLength(minimum) == expected, + "Selected FFT length is not the smallest allowed size"); + } + require(blah2::nextFastFftLength(480211) == 481140, + "200 ms convolution geometry changed"); + require(blah2::nextFastFftLength(1200211) == 1200500, + "500 ms convolution geometry changed"); + for (uint64_t invalid : {uint64_t(0), uint64_t(INT32_MAX), + uint64_t(INT32_MAX) + 1, UINT64_MAX}) { + bool rejected = false; + try { (void)blah2::nextFastFftLength(invalid); } + catch (const std::invalid_argument&) { rejected = true; } + require(rejected, "Unrepresentable FFT geometry accepted"); + } + for (unsigned samples : {64u, 127u, 257u}) + for (int first : {-3, 0, 2}) run(samples, first, first + 8); + run(128, -2, -1); + run(31, -2, 2); // Already-fast convolution length (36). + run(32, 0, 32); // Tap count equal to CPI sample count. + for (const auto& bounds : {std::pair{0, 0}, {5, 2}, {0, 65}, + {INT32_MIN, INT32_MAX}}) { + bool rejected = false; + try { WienerHopf invalid(bounds.first, bounds.second, 64); } + catch (const std::invalid_argument&) { rejected = true; } + require(rejected, "Invalid clutter tap range accepted"); + } + bool rejected = false; + try { WienerHopf invalid(0, 1, 0); } + catch (const std::invalid_argument&) { rejected = true; } + require(rejected, "Empty CPI accepted"); + IqData x(64), y(64); + for (unsigned i = 0; i < 64; ++i) { x.push_back({0, 0}); y.push_back({1, 0}); } + const auto original = y.get_data(); + WienerHopf filter(0, 8, 64); + require(!filter.process(&x, &y) && y.get_data() == original, + "Singular clutter input was accepted or changed"); + } catch (const std::exception& error) { + std::cerr << error.what() << '\n'; + return 1; + } +} From c4b3db91530ad22c0d856da6a4902e6df2ec3ba8 Mon Sep 17 00:00:00 2001 From: Josh Poole Date: Mon, 14 Sep 2026 11:29:01 +0100 Subject: [PATCH 2/3] 20260914 - Time the clutter FFT length instead of deriving it The previous commit padded the clutter convolution to the smallest length built from 2/3/5/7 (plus at most one 11 or 13). On x86 that measured 1.40x faster. On the Pi 5 it is a regression, and an A-B-A live swap on owl-ded9 against v0.4.3 measured it: A1 v0.4.3 n=76 clutter_filter 763.2 ms cpi 1211.5 ms B padded n=47 clutter_filter 828.8 ms cpi 1276.5 ms A2 v0.4.3 n=46 clutter_filter 758.7 ms cpi 1206.1 ms That is +67.9 ms (0.918x), Welch t = 12.7, against a -4.5 ms A1/A2 drift on identical code. v0.4.3 to this branch's parent changes no C++ at all, so the comparison isolates the padding. The cause is that blah2.cpp calls fftw_plan_with_nthreads(4), so every plan is a 4-thread plan. FFTW parallelises across a Cooley-Tukey factor, so how evenly the factors divide the thread count matters more than how small they are. Benchmarked on a Pi 5 against the image's own NEON FFTW, 2 forward plus 1 backward, best of 5: length 1 thread 4 threads 1016064 = 2^8*3^4*7^2 94.24 ms 1036800 = 2^9*3^4*5^2 140.23 ms 1048576 = 2^20 248.07 ms 151.18 ms 1000411 = 269*3719 466.03 ms 255.86 ms (unpadded) 1003520 = 2^12*5*7^2 172.30 ms 266.23 ms 1002375 = 3^6*5^3*11 273.55 ms 393.63 ms (what we padded to) Single-threaded the padding is a 1.7x win; threaded it is a 1.54x loss. The ordering is not predictable from the factorisation either way: 2^8*3^4*7^2 beats 2^20 by 1.6x, and 2^13*5^3 is worse than two large primes. Any static rule is guesswork, so stop writing one. fastestFftLength() times the admissible candidates within 2% of the minimum, plus the unpadded length itself, and keeps the quickest. It runs once in a constructor, after the thread count is set, in a process that ran 139 days without a restart on jn1. nextFastFftLength() stays: it is the fallback when nothing can be planned or allocated, and the test still pins its arithmetic. The clutter test now asserts the properties the convolution depends on (at least the alias-free length, and one of the offered candidates) rather than a fixed length that varies by machine, and covers the candidate enumeration directly. Co-Authored-By: Claude Opus 5 (1M context) --- CMakeLists.txt | 2 + src/process/clutter/WienerHopf.cpp | 7 +- src/process/meta/FftLength.cpp | 84 ++++++++++++++++++++ src/process/meta/FftLength.h | 18 +++++ test/unit/process/clutter/TestClutterFft.cpp | 49 +++++++++++- 5 files changed, 156 insertions(+), 4 deletions(-) create mode 100644 src/process/meta/FftLength.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index c401c75c..51f92509 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -53,6 +53,7 @@ add_executable(blah2 src/process/tracker/Tracker.cpp src/process/spectrum/SpectrumAnalyser.cpp src/process/meta/HammingNumber.cpp + src/process/meta/FftLength.cpp src/process/utility/Socket.cpp src/data/IqData.cpp src/data/Map.cpp @@ -115,6 +116,7 @@ set_target_properties(testHammingNumber PROPERTIES add_executable(testClutterFft test/unit/process/clutter/TestClutterFft.cpp src/process/clutter/WienerHopf.cpp + src/process/meta/FftLength.cpp src/data/IqData.cpp ) target_link_libraries(testClutterFft PRIVATE diff --git a/src/process/clutter/WienerHopf.cpp b/src/process/clutter/WienerHopf.cpp index 647a309f..5f028d57 100644 --- a/src/process/clutter/WienerHopf.cpp +++ b/src/process/clutter/WienerHopf.cpp @@ -15,8 +15,11 @@ WienerHopf::WienerHopf(int32_t _delayMin, int32_t _delayMax, uint32_t _nSamples) throw std::invalid_argument("Clutter filter needs a non-empty half-open delay range no longer than the CPI"); nBins = static_cast(taps); nSamples = _nSamples; - // Pad only the linear convolution; keep taps and circular correlations unchanged. - nFilter = blah2::nextFastFftLength(uint64_t(nSamples) + nBins + 1); + // Pad only the linear convolution; keep taps and circular correlations + // unchanged. The length is timed rather than derived: under + // fftw_plan_with_nthreads() the quickest size does not follow from the + // factorisation, and the ranking differs between machines and thread counts. + nFilter = blah2::fastestFftLength(uint64_t(nSamples) + nBins + 1); // initialise data A = arma::cx_mat(nBins, nBins); diff --git a/src/process/meta/FftLength.cpp b/src/process/meta/FftLength.cpp new file mode 100644 index 00000000..d0cd373f --- /dev/null +++ b/src/process/meta/FftLength.cpp @@ -0,0 +1,84 @@ +#include "FftLength.h" + +#include +#include +#include + +#include + +namespace blah2 { + +std::vector fftLengthCandidates(uint64_t minimum, double slack) { + constexpr uint64_t limit = std::numeric_limits::max(); + if (!minimum || minimum > limit) + throw std::invalid_argument("FFT length must fit a positive FFTW int"); + if (!(slack >= 0.0)) + throw std::invalid_argument("FFT length slack must not be negative"); + + const uint64_t ceiling = + std::min(limit, minimum + static_cast(double(minimum) * slack)); + + // The unpadded length always stays in the running; padding is only ever an + // optimisation, never a correctness requirement. + std::vector lengths{static_cast(minimum)}; + for (uint64_t extra : {1u, 11u, 13u}) + for (uint64_t a = extra; a <= ceiling; a *= 2) + for (uint64_t b = a; b <= ceiling; b *= 3) + for (uint64_t c = b; c <= ceiling; c *= 5) + for (uint64_t d = c; d <= ceiling; d *= 7) + if (d >= minimum) lengths.push_back(static_cast(d)); + + std::sort(lengths.begin(), lengths.end()); + lengths.erase(std::unique(lengths.begin(), lengths.end()), lengths.end()); + return lengths; +} + +uint32_t fastestFftLength(uint64_t minimum, double slack) { + const std::vector lengths = fftLengthCandidates(minimum, slack); + + uint32_t best = lengths.front(); + double bestSeconds = -1.0; + for (uint32_t n : lengths) { + auto *buffer = + static_cast(fftw_malloc(sizeof(fftw_complex) * size_t(n))); + if (buffer == nullptr) continue; + + fftw_plan forward = + fftw_plan_dft_1d(int(n), buffer, buffer, FFTW_FORWARD, FFTW_ESTIMATE); + fftw_plan backward = + fftw_plan_dft_1d(int(n), buffer, buffer, FFTW_BACKWARD, FFTW_ESTIMATE); + if (forward != nullptr && backward != nullptr) { + for (uint32_t i = 0; i < n; i++) { + buffer[i][0] = double(i) / double(n); + buffer[i][1] = 0.0; + } + fftw_execute(forward); // discard the first pass, which warms the caches + + const auto start = std::chrono::steady_clock::now(); + fftw_execute(forward); + fftw_execute(forward); + fftw_execute(backward); // the two-forward, one-backward mix of a CPI + const double seconds = + std::chrono::duration(std::chrono::steady_clock::now() - start) + .count(); + + if (bestSeconds < 0.0 || seconds < bestSeconds) { + bestSeconds = seconds; + best = n; + } + } + + if (forward != nullptr) fftw_destroy_plan(forward); + if (backward != nullptr) fftw_destroy_plan(backward); + fftw_free(buffer); + } + + // Nothing could be planned or allocated; fall back to the static choice. + if (bestSeconds < 0.0) return nextFastFftLength(minimum); + + std::cout << "Clutter filter FFT length " << best << " (from " << lengths.size() + << " candidates >= " << minimum << ", " << bestSeconds * 1000.0 + << " ms per CPI of transforms)" << std::endl; + return best; +} +} diff --git a/src/process/meta/FftLength.h b/src/process/meta/FftLength.h index 918ed760..6b7696d2 100644 --- a/src/process/meta/FftLength.h +++ b/src/process/meta/FftLength.h @@ -4,6 +4,7 @@ #include #include #include +#include namespace blah2 { // FFTW supports these small factors efficiently. Permit at most one factor of @@ -23,4 +24,21 @@ inline uint32_t nextFastFftLength(uint64_t minimum) { throw std::invalid_argument("No supported padded FFT length fits FFTW"); return static_cast(best); } + +// Padded lengths worth considering: the admissible sizes within `slack` of the +// minimum, plus the minimum itself, which is often quickest despite factoring +// badly. Sorted ascending. +std::vector fftLengthCandidates(uint64_t minimum, double slack = 0.02); + +// Times every candidate with the planner's current thread count and returns the +// quickest. Once fftw_plan_with_nthreads() is in play, transform cost stops +// being predictable from the factorisation: FFTW parallelises across a +// Cooley-Tukey factor, so how evenly the factors divide the thread count +// matters more than how small they are. Measured on a Pi 5 at the shipped +// clutter geometry, 2 forward plus 1 backward, the smallest admissible length +// (1002375 = 3^6*5^3*11) took 393.63 ms against 255.86 ms for the unpadded +// 1000411 = 269*3719, while 1016064 = 2^8*3^4*7^2 took 94.24 ms. Reversing the +// thread count reverses the ranking, so the length has to be measured on the +// machine and thread count that will run it. +uint32_t fastestFftLength(uint64_t minimum, double slack = 0.02); } diff --git a/test/unit/process/clutter/TestClutterFft.cpp b/test/unit/process/clutter/TestClutterFft.cpp index 171d53fc..c1dd09e3 100644 --- a/test/unit/process/clutter/TestClutterFft.cpp +++ b/test/unit/process/clutter/TestClutterFft.cpp @@ -1,8 +1,10 @@ #include "process/clutter/WienerHopf.h" #include "process/meta/FftLength.h" +#include #include #include #include +#include using Complex = std::complex; @@ -20,8 +22,16 @@ static void run(unsigned samples, int first, int last) { const unsigned taps = last - first; IqData reference(samples), surveillance(samples); WienerHopf filter(first, last, samples); - require(filter.filter_fft_length() == blah2::nextFastFftLength(samples + taps + 1), - "Filter did not use the selected padded length"); + // The length is timed at construction, so it is whichever candidate was + // quickest here rather than a fixed value. Both properties the convolution + // relies on must still hold. + const uint64_t minimum = uint64_t(samples) + taps + 1; + const std::vector allowed = blah2::fftLengthCandidates(minimum); + require(filter.filter_fft_length() >= minimum, + "Filter FFT length is below the alias-free convolution length"); + require(std::find(allowed.begin(), allowed.end(), filter.filter_fft_length()) + != allowed.end(), + "Filter did not use one of the candidate padded lengths"); std::mt19937 rng(9211); std::normal_distribution random; for (int repeat = 0; repeat < 3; ++repeat) { @@ -81,6 +91,41 @@ int main() { "200 ms convolution geometry changed"); require(blah2::nextFastFftLength(1200211) == 1200500, "500 ms convolution geometry changed"); + + // The timed selector may return any candidate, so pin the properties the + // convolution depends on rather than a length that varies by machine. + for (uint64_t minimum : {uint64_t(97), uint64_t(4096), uint64_t(480211), + uint64_t(1000411)}) { + const std::vector candidates = blah2::fftLengthCandidates(minimum); + require(!candidates.empty(), "No candidate FFT lengths offered"); + require(std::is_sorted(candidates.begin(), candidates.end()), + "Candidate FFT lengths are not sorted"); + require(std::adjacent_find(candidates.begin(), candidates.end()) + == candidates.end(), + "Candidate FFT lengths contain duplicates"); + require(candidates.front() == minimum, + "Unpadded length is not among the candidates"); + for (uint32_t candidate : candidates) { + require(candidate >= minimum, "Candidate is shorter than the minimum"); + require(candidate == minimum || fast(candidate), + "Padded candidate is not an admissible FFT length"); + } + require(std::find(candidates.begin(), candidates.end(), + blah2::nextFastFftLength(minimum)) != candidates.end(), + "Smallest admissible length is not among the candidates"); + const uint32_t timed = blah2::fastestFftLength(minimum); + require(std::find(candidates.begin(), candidates.end(), timed) + != candidates.end(), + "Timed FFT length is not one of the candidates"); + } + require(blah2::fftLengthCandidates(1000411, 0.0).size() >= 1, + "Zero slack must still offer the unpadded length"); + { + bool rejected = false; + try { (void)blah2::fftLengthCandidates(1000411, -0.5); } + catch (const std::invalid_argument&) { rejected = true; } + require(rejected, "Negative FFT length slack accepted"); + } for (uint64_t invalid : {uint64_t(0), uint64_t(INT32_MAX), uint64_t(INT32_MAX) + 1, UINT64_MAX}) { bool rejected = false; From c453c57ac73ec05d2ea1df4a8482478caa523447 Mon Sep 17 00:00:00 2001 From: Josh Poole Date: Mon, 14 Sep 2026 11:51:33 +0100 Subject: [PATCH 3/3] 20260914 - Remember the measured clutter FFT length between runs Timing the candidate lengths costs about 7 s on a Pi 5 (15 candidates, each planned and run four times at ~1e6 points), paid on every start on top of the existing sleep(5). That roughly doubles time to first CPI, which matters most on the nodes that restart most. The answer only changes when the geometry, the thread count or the FFTW build changes, so cache it against exactly those three. The file lives in the save directory, which is a host mount, so it survives a container restart; BLAH2_FFT_CACHE overrides the path. Only the first start after a change pays the sweep. A few entries are kept rather than one, so moving a node between configurations and back does not force a re-measurement each way. The cache is an optimisation and never a contract. A corrupt line, an unreadable file, a read-only or missing save mount, or a remembered length shorter than the geometry now needs, all fall back to measuring. Nothing about correctness rides on it: every candidate is at least the alias-free convolution length, so a stale entry costs speed, never accuracy. fftw_plan_with_nthreads(4) becomes kPlannerThreads, shared with the cache key, since the ranking inverts between 1 and 4 threads and the key would otherwise silently lie. Tests cover the round trip (and assert the cached path is quicker than measuring), isolation between geometries, no eviction when a second geometry is added, recovery from a corrupt file, and an unwritable path. Built and run on a Pi 5 against the deployed FFTW: all geometries pass, including the positive delayMin case that exercises the signed-index fix. Co-Authored-By: Claude Opus 5 (1M context) --- src/blah2.cpp | 3 +- src/process/meta/FftLength.cpp | 96 ++++++++++++++++++-- src/process/meta/FftLength.h | 39 +++++--- test/unit/process/clutter/TestClutterFft.cpp | 51 +++++++++++ 4 files changed, 170 insertions(+), 19 deletions(-) diff --git a/src/blah2.cpp b/src/blah2.cpp index 8a130562..121b0398 100644 --- a/src/blah2.cpp +++ b/src/blah2.cpp @@ -13,6 +13,7 @@ #include "process/detection/CfarDetector1D.h" #include "process/detection/Centroid.h" #include "process/detection/Interpolate.h" +#include "process/meta/FftLength.h" #include "process/spectrum/SpectrumAnalyser.h" #include "process/tracker/Tracker.h" #include "process/utility/Socket.h" @@ -124,7 +125,7 @@ int main(int argc, char **argv) std::cout << "Error in FFTW multithreading." << std::endl; return -1; } - fftw_plan_with_nthreads(4); + fftw_plan_with_nthreads(blah2::kPlannerThreads); // setup socket sleep(5); diff --git a/src/process/meta/FftLength.cpp b/src/process/meta/FftLength.cpp index d0cd373f..d76dd170 100644 --- a/src/process/meta/FftLength.cpp +++ b/src/process/meta/FftLength.cpp @@ -2,11 +2,71 @@ #include #include +#include +#include +#include #include +#include +#include #include namespace blah2 { +namespace { + +// Keep a handful of geometries so that moving a node between configurations, +// or back again, does not force a re-measurement each time. +constexpr size_t kMaxCacheEntries = 8; + +// Everything that can change which length wins. The FFTW build is in here +// because the codelets it ships decide the ranking; the thread count is in here +// because the ranking inverts between 1 and 4 threads. +std::string cacheKey(uint64_t minimum, double slack, int threads) { + std::ostringstream key; + key << minimum << ' ' << std::fixed << std::setprecision(6) << slack << ' ' + << threads << ' ' << fftw_version; + return key.str(); +} + +// " ", with the key last so an FFTW version string containing +// spaces still round-trips. +std::vector> readCache(const std::string& path) { + std::vector> entries; + std::ifstream file(path); + std::string line; + while (entries.size() < kMaxCacheEntries && std::getline(file, line)) { + const size_t split = line.find(' '); + if (split == std::string::npos) continue; + uint64_t length = 0; + std::istringstream parse(line.substr(0, split)); + if (!(parse >> length) || length == 0 || + length > uint64_t(std::numeric_limits::max())) + continue; // a corrupt line only ever costs a re-measurement + entries.emplace_back(static_cast(length), line.substr(split + 1)); + } + return entries; +} + +void writeCache(const std::string& path, + std::vector> entries, + uint32_t length, const std::string& key) { + entries.erase(std::remove_if(entries.begin(), entries.end(), + [&key](const std::pair& entry) { + return entry.second == key; + }), + entries.end()); + entries.emplace_back(length, key); + if (entries.size() > kMaxCacheEntries) + entries.erase(entries.begin(), entries.end() - kMaxCacheEntries); + + // A missing or read-only save directory just means measuring every start. + std::ofstream file(path, std::ios::trunc); + if (!file) return; + for (const std::pair& entry : entries) + file << entry.first << ' ' << entry.second << '\n'; +} + +} // namespace std::vector fftLengthCandidates(uint64_t minimum, double slack) { constexpr uint64_t limit = std::numeric_limits::max(); @@ -33,14 +93,33 @@ std::vector fftLengthCandidates(uint64_t minimum, double slack) { return lengths; } -uint32_t fastestFftLength(uint64_t minimum, double slack) { +std::string fftLengthCachePath() { + if (const char* override = std::getenv("BLAH2_FFT_CACHE")) + if (*override != '\0') return override; + return "/opt/blah2/save/fft-length.cache"; +} + +uint32_t fastestFftLength(uint64_t minimum, double slack, int threads) { const std::vector lengths = fftLengthCandidates(minimum, slack); + const std::string path = fftLengthCachePath(); + const std::string key = cacheKey(minimum, slack, threads); + const std::vector> cached = readCache(path); + + for (const std::pair& entry : cached) { + // Only trust a remembered length that is still long enough to filter + // correctly, whatever else may have changed. + if (entry.second == key && entry.first >= minimum) { + std::cout << "Clutter filter FFT length " << entry.first << " (cached in " + << path << ")" << std::endl; + return entry.first; + } + } uint32_t best = lengths.front(); double bestSeconds = -1.0; for (uint32_t n : lengths) { - auto *buffer = - static_cast(fftw_malloc(sizeof(fftw_complex) * size_t(n))); + auto* buffer = + static_cast(fftw_malloc(sizeof(fftw_complex) * size_t(n))); if (buffer == nullptr) continue; fftw_plan forward = @@ -73,12 +152,15 @@ uint32_t fastestFftLength(uint64_t minimum, double slack) { fftw_free(buffer); } - // Nothing could be planned or allocated; fall back to the static choice. + // Nothing could be planned or allocated; fall back to the static choice + // without poisoning the cache with a length we never measured. if (bestSeconds < 0.0) return nextFastFftLength(minimum); - std::cout << "Clutter filter FFT length " << best << " (from " << lengths.size() - << " candidates >= " << minimum << ", " << bestSeconds * 1000.0 - << " ms per CPI of transforms)" << std::endl; + writeCache(path, cached, best, key); + std::cout << "Clutter filter FFT length " << best << " (measured from " + << lengths.size() << " candidates >= " << minimum << ", " + << bestSeconds * 1000.0 << " ms per CPI of transforms, cached in " + << path << ")" << std::endl; return best; } } diff --git a/src/process/meta/FftLength.h b/src/process/meta/FftLength.h index 6b7696d2..8d0eacb1 100644 --- a/src/process/meta/FftLength.h +++ b/src/process/meta/FftLength.h @@ -4,9 +4,14 @@ #include #include #include +#include #include namespace blah2 { +// Threads handed to the FFTW planner. Which transform length is quickest +// depends on this, so it belongs anywhere the choice is made or remembered. +inline constexpr int kPlannerThreads = 4; + // FFTW supports these small factors efficiently. Permit at most one factor of // 11 or 13. Enumerate bounded candidates instead of an unbounded integer scan. inline uint32_t nextFastFftLength(uint64_t minimum) { @@ -30,15 +35,27 @@ inline uint32_t nextFastFftLength(uint64_t minimum) { // badly. Sorted ascending. std::vector fftLengthCandidates(uint64_t minimum, double slack = 0.02); -// Times every candidate with the planner's current thread count and returns the -// quickest. Once fftw_plan_with_nthreads() is in play, transform cost stops -// being predictable from the factorisation: FFTW parallelises across a -// Cooley-Tukey factor, so how evenly the factors divide the thread count -// matters more than how small they are. Measured on a Pi 5 at the shipped -// clutter geometry, 2 forward plus 1 backward, the smallest admissible length -// (1002375 = 3^6*5^3*11) took 393.63 ms against 255.86 ms for the unpadded -// 1000411 = 269*3719, while 1016064 = 2^8*3^4*7^2 took 94.24 ms. Reversing the -// thread count reverses the ranking, so the length has to be measured on the -// machine and thread count that will run it. -uint32_t fastestFftLength(uint64_t minimum, double slack = 0.02); +// Where measured lengths are remembered between runs. BLAH2_FFT_CACHE overrides +// the default, which sits in the save directory so it survives a container +// restart. +std::string fftLengthCachePath(); + +// Returns the quickest candidate, measuring only when the answer is not already +// cached for this geometry, thread count and FFTW build. +// +// Once fftw_plan_with_nthreads() is in play, transform cost stops being +// predictable from the factorisation: FFTW parallelises across a Cooley-Tukey +// factor, so how evenly the factors divide the thread count matters more than +// how small they are. Measured on a Pi 5 at the shipped clutter geometry, 2 +// forward plus 1 backward, the smallest admissible length (1002375 = +// 3^6*5^3*11) took 393.63 ms against 255.86 ms for the unpadded 1000411 = +// 269*3719, while 1016064 = 2^8*3^4*7^2 took 94.24 ms. Reversing the thread +// count reverses the ranking, so the length has to be measured on the machine +// and thread count that will run it. +// +// A stale or unreadable cache costs speed, never correctness: every candidate +// is at least the alias-free convolution length, so any of them filters +// correctly. +uint32_t fastestFftLength(uint64_t minimum, double slack = 0.02, + int threads = kPlannerThreads); } diff --git a/test/unit/process/clutter/TestClutterFft.cpp b/test/unit/process/clutter/TestClutterFft.cpp index c1dd09e3..52c13898 100644 --- a/test/unit/process/clutter/TestClutterFft.cpp +++ b/test/unit/process/clutter/TestClutterFft.cpp @@ -1,9 +1,15 @@ #include "process/clutter/WienerHopf.h" #include "process/meta/FftLength.h" #include +#include +#include +#include +#include #include #include #include +#include +#include #include using Complex = std::complex; @@ -126,6 +132,51 @@ int main() { catch (const std::invalid_argument&) { rejected = true; } require(rejected, "Negative FFT length slack accepted"); } + + { + // Cache round-trip. Measuring is what makes startup slow, so a second + // call with the same geometry must reuse the stored answer rather than + // repeat the sweep. + const std::string cache = "/tmp/blah2-fft-length-test.cache"; + std::remove(cache.c_str()); + setenv("BLAH2_FFT_CACHE", cache.c_str(), 1); + require(blah2::fftLengthCachePath() == cache, + "BLAH2_FFT_CACHE did not override the cache path"); + + const auto timed = [](uint64_t minimum) { + const auto start = std::chrono::steady_clock::now(); + const uint32_t length = blah2::fastestFftLength(minimum); + return std::pair{length, + std::chrono::duration(std::chrono::steady_clock::now() - start) + .count()}; + }; + + const std::pair cold = timed(40009); + const std::pair warm = timed(40009); + require(warm.first == cold.first, "Cached FFT length differs from measured"); + require(warm.second < cold.second, + "Cached lookup was no quicker than measuring"); + require(std::ifstream(cache).good(), "Cache file was not written"); + + // A different geometry must not collide with the stored entry. + const uint32_t other = blah2::fastestFftLength(50021); + require(other >= 50021, "Cache returned a length below a new minimum"); + require(blah2::fastestFftLength(40009) == cold.first, + "Adding a geometry evicted the earlier one"); + + // Garbage must be survivable: a cache is an optimisation, not a contract. + { std::ofstream(cache, std::ios::trunc) << "not-a-number\n\nx y z\n"; } + require(blah2::fastestFftLength(40009) >= 40009, + "Corrupt cache was not recovered from"); + + // An unwritable location must not stop the radar starting. + setenv("BLAH2_FFT_CACHE", "/nonexistent-directory/fft.cache", 1); + require(blah2::fastestFftLength(40009) >= 40009, + "Unwritable cache path was not tolerated"); + + unsetenv("BLAH2_FFT_CACHE"); + std::remove(cache.c_str()); + } for (uint64_t invalid : {uint64_t(0), uint64_t(INT32_MAX), uint64_t(INT32_MAX) + 1, UINT64_MAX}) { bool rejected = false;