Choose the clutter convolution FFT length by measurement

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
@@ -112,6 +113,21 @@ 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/process/meta/FftLength.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)

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);

src/process/clutter/WienerHopf.cpp

@@ -1,4 +1,5 @@
 #include "WienerHopf.h"
+#include "process/meta/FftLength.h"
 #include <complex>
 #include <iostream>
 #include <vector>
@@ -9,8 +10,16 @@ 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<uint32_t>(taps);
   nSamples = _nSamples;
+  // 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);
@@ -25,9 +34,9 @@ WienerHopf::WienerHopf(int32_t _delayMin, int32_t _delayMax, uint32_t _nSamples)
   dataOutY = new std::complex<double>[nSamples];
   dataA = new std::complex<double>[nSamples];
   dataB = new std::complex<double>[nSamples];
-  filtX = new std::complex<double>[nBins + nSamples + 1];
-  filtW = new std::complex<double>[nBins + nSamples + 1];
-  filt = new std::complex<double>[nBins + nSamples + 1];
+  filtX = new std::complex<double>[nFilter];
+  filtW = new std::complex<double>[nFilter];
+  filt = new std::complex<double>[nFilter];
   fftX = fftw_plan_dft_1d(nSamples, reinterpret_cast<fftw_complex *>(dataX),
                           reinterpret_cast<fftw_complex *>(dataOutX), FFTW_FORWARD, FFTW_ESTIMATE);
   fftY = fftw_plan_dft_1d(nSamples, reinterpret_cast<fftw_complex *>(dataY),
@@ -36,11 +45,11 @@ WienerHopf::WienerHopf(int32_t _delayMin, int32_t _delayMax, uint32_t _nSamples)
                           reinterpret_cast<fftw_complex *>(dataA), FFTW_BACKWARD, FFTW_ESTIMATE);
   fftB = fftw_plan_dft_1d(nSamples, reinterpret_cast<fftw_complex *>(dataB),
                           reinterpret_cast<fftw_complex *>(dataB), FFTW_BACKWARD, FFTW_ESTIMATE);
-  fftFiltX = fftw_plan_dft_1d(nBins + nSamples + 1, reinterpret_cast<fftw_complex *>(filtX),
+  fftFiltX = fftw_plan_dft_1d(nFilter, reinterpret_cast<fftw_complex *>(filtX),
                               reinterpret_cast<fftw_complex *>(filtX), FFTW_FORWARD, FFTW_ESTIMATE);
-  fftFiltW = fftw_plan_dft_1d(nBins + nSamples + 1, reinterpret_cast<fftw_complex *>(filtW),
+  fftFiltW = fftw_plan_dft_1d(nFilter, reinterpret_cast<fftw_complex *>(filtW),
                               reinterpret_cast<fftw_complex *>(filtW), FFTW_FORWARD, FFTW_ESTIMATE);
-  fftFilt = fftw_plan_dft_1d(nBins + nSamples + 1, reinterpret_cast<fftw_complex *>(filt),
+  fftFilt = fftw_plan_dft_1d(nFilter, reinterpret_cast<fftw_complex *>(filt),
                              reinterpret_cast<fftw_complex *>(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;

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();

src/process/meta/FftLength.cpp

@@ -0,0 +1,166 @@
+#include "FftLength.h"
+
+#include <algorithm>
+#include <chrono>
+#include <cstdlib>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <sstream>
+#include <utility>
+
+#include <fftw3.h>
+
+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();
+}
+
+// "<length> <key>", with the key last so an FFTW version string containing
+// spaces still round-trips.
+std::vector<std::pair<uint32_t, std::string>> readCache(const std::string& path) {
+  std::vector<std::pair<uint32_t, std::string>> 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<int>::max()))
+      continue;  // a corrupt line only ever costs a re-measurement
+    entries.emplace_back(static_cast<uint32_t>(length), line.substr(split + 1));
+  }
+  return entries;
+}
+
+void writeCache(const std::string& path,
+                std::vector<std::pair<uint32_t, std::string>> entries,
+                uint32_t length, const std::string& key) {
+  entries.erase(std::remove_if(entries.begin(), entries.end(),
+                               [&key](const std::pair<uint32_t, std::string>& 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<uint32_t, std::string>& entry : entries)
+    file << entry.first << ' ' << entry.second << '\n';
+}
+
+}  // namespace
+
+std::vector<uint32_t> fftLengthCandidates(uint64_t minimum, double slack) {
+  constexpr uint64_t limit = std::numeric_limits<int>::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<uint64_t>(double(minimum) * slack));
+
+  // The unpadded length always stays in the running; padding is only ever an
+  // optimisation, never a correctness requirement.
+  std::vector<uint32_t> lengths{static_cast<uint32_t>(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<uint32_t>(d));
+
+  std::sort(lengths.begin(), lengths.end());
+  lengths.erase(std::unique(lengths.begin(), lengths.end()), lengths.end());
+  return lengths;
+}
+
+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<uint32_t> lengths = fftLengthCandidates(minimum, slack);
+  const std::string path = fftLengthCachePath();
+  const std::string key = cacheKey(minimum, slack, threads);
+  const std::vector<std::pair<uint32_t, std::string>> cached = readCache(path);
+
+  for (const std::pair<uint32_t, std::string>& 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_complex*>(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<double>(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
+  // without poisoning the cache with a length we never measured.
+  if (bestSeconds < 0.0) return nextFastFftLength(minimum);
+
+  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;
+}
+}

src/process/meta/FftLength.h

@@ -0,0 +1,61 @@
+#pragma once
+
+#include <cstdint>
+#include <initializer_list>
+#include <limits>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+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) {
+  constexpr uint64_t limit = std::numeric_limits<int>::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<uint32_t>(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<uint32_t> fftLengthCandidates(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);
+}