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..791fe571 100644 --- a/src/process/clutter/WienerHopf.cpp +++ b/src/process/clutter/WienerHopf.cpp @@ -1,5 +1,7 @@ #include "WienerHopf.h" +#include "process/meta/FftLength.h" #include +#include #include #include @@ -9,8 +11,15 @@ 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; taps and circular correlations are + // unchanged. Any length at or above nSamples + nBins - 1 is alias-free and + // only the first nSamples outputs are read, so rounding up is free. + nFilter = blah2::clutterFftLength(uint64_t(nSamples) + nBins + 1); // initialise data A = arma::cx_mat(nBins, nBins); @@ -25,9 +34,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 +45,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 +73,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 +138,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 +148,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 +158,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 +168,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..8d92f362 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 Padded length of the convolution FFTs. + uint32_t nFilter; + /// @brief True if clutter filter processing is successful. bool success; @@ -67,6 +70,10 @@ class WienerHopf /// @return The object. WienerHopf(int32_t delayMin, int32_t delayMax, uint32_t nSamples); + /// @brief Padded length the convolution FFTs are planned at. + /// @return Number of points. + 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..23ba1a61 --- /dev/null +++ b/src/process/meta/FftLength.h @@ -0,0 +1,56 @@ +#pragma once + +#include +#include +#include + +namespace blah2 { + +// The length the clutter convolution is padded to. +// +// 1016064 = 2^8 * 3^4 * 7^2. Measured, not derived. The shipped geometry +// (fs 2e6, cpi 0.5, clutter delay -10..400) needs at least 1000411 points, and +// that factors as 269 * 3719, both prime, so FFTW falls back to Rader and the +// transform costs roughly three times what it should. +// +// This value was picked by timing every admissible length within 2% of the +// minimum on the target hardware. It came first out of fifteen candidates on +// four different Pi 5 boards, across both memory variants, at every thread +// count from one to four, idle and under load, by a margin of 2.25x to 3.8x +// over not padding. Best and median times agreed throughout. +// +// It is deliberately a constant rather than something derived. A formula was +// tried first and was worse than doing nothing: picking the *smallest* +// admissible length gives 1002375 = 3^6 * 5^3 * 11, which measured 1.54x +// SLOWER than unpadded on ARM at four threads while looking like a 1.40x win +// on x86. Transform cost is not predictable from the factorisation, so the +// only honest way to choose is to measure, and the only way to keep a measured +// answer stable is to write it down. +// +// Re-measure if the geometry, the FFTW build or the target SoC changes. +inline constexpr uint32_t kClutterFftLength = 1016064; + +// The padded length to use for a convolution needing at least `minimum` points. +// +// Padding is only ever an optimisation: any length at or above the alias-free +// minimum produces the identical linear convolution, and only the first +// nSamples outputs are read. So the constant is used only where it genuinely +// covers the geometry without being wastefully larger, and anything else falls +// back to the unpadded length, which is always correct and is exactly what this +// code did before the constant existed. +inline uint32_t clutterFftLength(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"); + + // Same 2% band the candidates were drawn from. At the shipped geometry the + // constant sits 1.56% above the minimum, so it applies; change cpi or the + // delay span far enough and it quietly stops applying rather than aliasing. + if (kClutterFftLength >= minimum && kClutterFftLength <= minimum + minimum / 50) + return kClutterFftLength; + + return static_cast(minimum); +} + +} diff --git a/test/unit/process/clutter/TestClutterFft.cpp b/test/unit/process/clutter/TestClutterFft.cpp new file mode 100644 index 00000000..3efa7d62 --- /dev/null +++ b/test/unit/process/clutter/TestClutterFft.cpp @@ -0,0 +1,145 @@ +// The padded clutter convolution. +// +// Two things need holding down. That the padded length is always at least the +// alias-free convolution length, since anything shorter wraps and corrupts the +// start of the output. And that the filter still computes the right answer, +// checked against an independent direct convolution rather than against another +// FFT, so a fault in the transform path cannot hide behind itself. +#include "process/clutter/WienerHopf.h" +#include "process/meta/FftLength.h" + +#include +#include +#include +#include +#include +#include +#include + +using Complex = std::complex; + +static void require(bool condition, const char* message) { + if (!condition) throw std::runtime_error(message); +} + +static void run(unsigned samples, int first, int last) { + const unsigned taps = last - first; + IqData reference(samples), surveillance(samples); + WienerHopf filter(first, last, samples); + + // The only property the convolution depends on: long enough to be alias-free. + const uint64_t minimum = uint64_t(samples) + taps + 1; + require(filter.filter_fft_length() >= minimum, + "Filter FFT length is below the alias-free convolution 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 a direct linear convolution, not + // another FFT implementation. The filter is reused across three 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 { + // The shipped geometry is the one the constant was measured for: 2 MS/s, + // 0.5 s CPI, clutter delay -10..400, so 1000000 + 410 + 1 points. + const uint64_t shipped = 1000411; + require(blah2::clutterFftLength(shipped) == blah2::kClutterFftLength, + "Shipped geometry no longer uses the measured length"); + require(blah2::kClutterFftLength >= shipped, + "Measured length is shorter than the shipped geometry needs"); + std::cout << "PASS shipped geometry uses " << blah2::kClutterFftLength << '\n'; + + // The guard. The constant is only right for the geometry it was measured + // at, so anything it does not cover falls back to the unpadded length, + // which is always alias-free and is what the code did before it existed. + require(blah2::clutterFftLength(blah2::kClutterFftLength) == blah2::kClutterFftLength, + "Exact fit rejected"); + require(blah2::clutterFftLength(blah2::kClutterFftLength + 1) + == blah2::kClutterFftLength + 1, + "Geometry above the constant must not use it, that would alias"); + require(blah2::clutterFftLength(200411) == 200411, + "A much smaller geometry must not be padded fivefold"); + require(blah2::clutterFftLength(2000411) == 2000411, + "A larger geometry must fall back to unpadded"); + for (uint64_t minimum : {uint64_t(1), uint64_t(97), uint64_t(200411), + uint64_t(999999), shipped, uint64_t(2000411)}) + require(blah2::clutterFftLength(minimum) >= minimum, + "Padded length below the alias-free minimum"); + std::cout << "PASS geometry guard\n"; + + // A length outside what FFTW can index is a programming error, not a + // silently truncated transform. + for (uint64_t invalid : {uint64_t(0), uint64_t(INT32_MAX) + 1, UINT64_MAX}) { + bool rejected = false; + try { (void)blah2::clutterFftLength(invalid); } + catch (const std::invalid_argument&) { rejected = true; } + require(rejected, "Unrepresentable FFT geometry accepted"); + } + std::cout << "PASS invalid geometry rejected\n"; + + // Correctness across a spread of geometries, including positive delayMin, + // which used to read the wrong sample: `i - delayMin` promoted to unsigned + // and wrapped at 2^32. Negative and zero cancelled exactly, so the shipped + // -10 was unaffected and the fault stayed hidden. + 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); + run(32, 0, 32); // tap count equal to the CPI sample count + + for (const auto& bounds : {std::pair{0, 0}, {5, 2}, {0, 65}, + {INT32_MIN, INT32_MAX}}) { + bool rejected = false; + try { WienerHopf reject(bounds.first, bounds.second, 64); } + catch (const std::invalid_argument&) { rejected = true; } + require(rejected, "Degenerate delay range accepted"); + } + std::cout << "PASS degenerate delay ranges rejected\n"; + + std::cout << "All clutter FFT checks passed\n"; + } catch (const std::exception& error) { + std::cerr << "FAIL: " << error.what() << '\n'; + return 1; + } + return 0; +}