diff --git a/CMakeLists.txt b/CMakeLists.txt index 2e8cf8c4..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 @@ -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) 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/clutter/WienerHopf.cpp b/src/process/clutter/WienerHopf.cpp index 9fa9ddc5..5f028d57 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,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(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[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..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.cpp b/src/process/meta/FftLength.cpp new file mode 100644 index 00000000..d76dd170 --- /dev/null +++ b/src/process/meta/FftLength.cpp @@ -0,0 +1,166 @@ +#include "FftLength.h" + +#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(); + 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; +} + +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))); + 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 + // 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; +} +} diff --git a/src/process/meta/FftLength.h b/src/process/meta/FftLength.h new file mode 100644 index 00000000..8d0eacb1 --- /dev/null +++ b/src/process/meta/FftLength.h @@ -0,0 +1,61 @@ +#pragma once + +#include +#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) { + 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); +} + +// 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); + +// 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 new file mode 100644 index 00000000..52c13898 --- /dev/null +++ b/test/unit/process/clutter/TestClutterFft.cpp @@ -0,0 +1,213 @@ +#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; + +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); + // 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) { + 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"); + + // 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"); + } + + { + // 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; + 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; + } +}