From aaea9552f14005983dd7644f19168fee350ed04f Mon Sep 17 00:00:00 2001 From: Josh Poole Date: Mon, 14 Sep 2026 14:08:56 +0100 Subject: [PATCH] 20260914 - Process each CPI in two overlapped stages The radar loop ran every stage of a CPI back to back, so throughput was the sum of the stages: 828 ms on owl-ded9 against a 500 ms CPI, and 688 ms on fairforest-b. Splitting the loop at the clutter filter makes throughput the slowest stage instead. The clutter filter holds no state across CPIs, so the split itself needs no numerical argument. Slots carry one CPI's channels and timing from the front stage to the back, two of them so the stages overlap fully, and everything that does hold state across CPIs (the tracker) stays inside the back stage, which processes CPIs one at a time in capture order. Three details the split has to get right: - The spectrum analyser moves from before the clutter filter to after it, so the slow half is extract plus filter rather than extract plus spectrum plus filter. It reads the reference channel, which the filter only reads, so its output does not change. TestStageOrder proves that on a signal the filter cancels by 58.8 dB, rather than by reading the source. - A live retune is latched against the CPI it arrives on instead of acted on where it is seen, because the tracker it resets now runs in the other stage. A CPI dropped by a failed filter hands the retune back rather than eating it. - FFTW bakes the thread count into the plan, so the stages plan with 2 threads each rather than both claiming all 4. The clutter transforms prefer 2 to 4 anyway (118.2 ms against 146.8 per CPI at 1e6 points), so the front stage gains; ambiguity gives up about 21 ms. Concurrent threaded executes on distinct plans were checked for deadlock and for identical results first. That last one is the only part that is not bit-identical: FFTW parallelises across a Cooley-Tukey factor, so a different thread count sums in a different order. Measured on aarch64, the disagreement is 2.4e-16 of peak at 1e6 points and 6.4e-21 at 200k, against the 6e-14 the FFT length work already accepted. On x86 the two thread counts agree exactly. TestStageOrder holds the bound. "cpi" still reports the work done on a CPI across both stages, excluding the queue wait, so it stays the quantity the serial loop reported. Throughput is now a separate "cpi_interval", and "pipeline_wait" says which stage is the bottleneck. Tests compile and pass on aarch64 as well as x86. Running the equivalence suite on ARM is what caught that IqData serialises its never-assigned min, max and mean, so whole-JSON comparison of two IqData objects is allocator luck; the test compares the spectrum and frequency arrays instead. That bug is pre-existing and still unfixed. Not yet measured on hardware. Expected throughput is the front stage: about 497 ms on owl-ded9 and 412 ms on fairforest-b, against 828 and 688 serial. That is not a reliable 500 ms on owl-ded9, whose clutter filter swings 13.6% between identical runs, so it needs the clutter filter work to follow. Co-Authored-By: Claude Opus 5 (1M context) --- CMakeLists.txt | 29 ++ src/blah2.cpp | 306 +++++++++++------- src/process/clutter/WienerHopf.cpp | 5 +- src/process/clutter/WienerHopf.h | 7 +- src/process/meta/FftLength.h | 9 + src/process/utility/CpiPipeline.h | 107 ++++++ test/unit/process/clutter/TestStageOrder.cpp | 206 ++++++++++++ test/unit/process/utility/TestCpiPipeline.cpp | 164 ++++++++++ 8 files changed, 718 insertions(+), 115 deletions(-) create mode 100644 src/process/utility/CpiPipeline.h create mode 100644 test/unit/process/clutter/TestStageOrder.cpp create mode 100644 test/unit/process/utility/TestCpiPipeline.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 2a7326f1..2e6d18a4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -166,6 +166,33 @@ target_link_libraries(testCfarHoist PRIVATE set_target_properties(testCfarHoist PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${PROJECT_BINARY_TEST_UNIT_DIR}") +add_executable(testCpiPipeline + test/unit/process/utility/TestCpiPipeline.cpp + src/data/IqData.cpp +) +# armadillo is linked only to pick up the vcpkg include path, which is +# where rapidjson lives; RAPIDJSON_INCLUDE_DIRS is NOTFOUND in this build. +target_link_libraries(testCpiPipeline PRIVATE + armadillo +) +set_target_properties(testCpiPipeline PROPERTIES + RUNTIME_OUTPUT_DIRECTORY "${PROJECT_BINARY_TEST_UNIT_DIR}") + +add_executable(testStageOrder + test/unit/process/clutter/TestStageOrder.cpp + src/process/clutter/WienerHopf.cpp + src/process/spectrum/SpectrumAnalyser.cpp + src/process/meta/FftLength.cpp + src/data/IqData.cpp +) +target_link_libraries(testStageOrder PRIVATE + armadillo + fftw3 + fftw3_threads +) +set_target_properties(testStageOrder 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) @@ -173,3 +200,5 @@ add_test(NAME testClutterFft COMMAND testClutterFft) add_test(NAME testJsonKm COMMAND testJsonKm) add_test(NAME testAmbiguityIndexing COMMAND testAmbiguityIndexing) add_test(NAME testCfarHoist COMMAND testCfarHoist) +add_test(NAME testCpiPipeline COMMAND testCpiPipeline) +add_test(NAME testStageOrder COMMAND testStageOrder) diff --git a/src/blah2.cpp b/src/blah2.cpp index 9e52879f..ef72f64e 100644 --- a/src/blah2.cpp +++ b/src/blah2.cpp @@ -16,6 +16,7 @@ #include "process/meta/FftLength.h" #include "process/spectrum/SpectrumAnalyser.h" #include "process/tracker/Tracker.h" +#include "process/utility/CpiPipeline.h" #include "process/utility/Socket.h" #include "process/utility/TuneState.h" #include "data/meta/Constants.h" @@ -111,13 +112,19 @@ int main(int argc, char **argv) // setup process CPI uint32_t nSamples = fs * tCpi; - IqData *x = new IqData(nSamples); - IqData *y = new IqData(nSamples); - Map> *map; - std::unique_ptr detection; - std::unique_ptr detection1; - std::unique_ptr detection2; - std::unique_ptr track; + + // Slots in flight between the two processing stages. Two is enough for full + // overlap: the front stage fills one while the back stage drains the other, + // and a stage that runs ahead blocks on the free list. + constexpr size_t kPipelineDepth = 2; + std::vector> slots; + blah2::BlockingQueue freeSlots; + blah2::BlockingQueue filteredSlots; + for (size_t i = 0; i < kPipelineDepth; i++) + { + slots.push_back(std::make_unique(nSamples)); + freeSlots.push(slots.back().get()); + } // setup fftw multithread if (fftw_init_threads() == 0) @@ -125,7 +132,6 @@ int main(int argc, char **argv) std::cout << "Error in FFTW multithreading." << std::endl; return -1; } - fftw_plan_with_nthreads(blah2::kPlannerThreads); // setup socket sleep(5); @@ -154,14 +160,20 @@ int main(int argc, char **argv) tree["process"]["ambiguity"]["delayMax"] >> delayMax; tree["process"]["ambiguity"]["dopplerMin"] >> dopplerMin; tree["process"]["ambiguity"]["dopplerMax"] >> dopplerMax; - Ambiguity *ambiguity = new Ambiguity(delayMin, delayMax, + // FFTW bakes the thread count into the plan, so the two stages get their own + // share of the cores by planning with different counts. Ambiguity and the + // spectrum analyser run in the back stage, the clutter filter in the front. + fftw_plan_with_nthreads(blah2::kBackStageThreads); + Ambiguity *ambiguity = new Ambiguity(delayMin, delayMax, dopplerMin, dopplerMax, fs, nSamples, roundHamming); // setup process clutter int32_t delayMinClutter, delayMaxClutter; tree["process"]["clutter"]["delayMin"] >> delayMinClutter; tree["process"]["clutter"]["delayMax"] >> delayMaxClutter; - WienerHopf *filter = new WienerHopf(delayMinClutter, delayMaxClutter, nSamples); + fftw_plan_with_nthreads(blah2::kFrontStageThreads); + WienerHopf *filter = new WienerHopf(delayMinClutter, delayMaxClutter, nSamples, + blah2::kFrontStageThreads); // setup process detection double pfa, minDoppler; @@ -194,6 +206,7 @@ int main(int argc, char **argv) // setup process spectrum analyser double spectrumBandwidth = 2000; + fftw_plan_with_nthreads(blah2::kBackStageThreads); SpectrumAnalyser *spectrumAnalyser = new SpectrumAnalyser(nSamples, spectrumBandwidth); // process options @@ -226,139 +239,208 @@ int main(int argc, char **argv) // setup output timing uint64_t tStart = current_time_ms(); Timing *timing = new Timing(tStart); - std::vector timing_name; - std::vector timing_time; - std::string jsonTiming; - std::vector time; - - // setup output json - std::string mapJson, detectionJson, jsonTracker, jsonIqData; - // run process + // Front stage: pull a CPI out of the capture buffers and clutter filter it. + // This is the slow half, so it sets the throughput of the whole radar. std::thread t2([&]{ while (true) { - buffer1->lock(); - buffer2->lock(); - if ((buffer1->get_length() > nSamples) && (buffer2->get_length() > nSamples)) + blah2::CpiSlot *slot = freeSlots.pop(); + + // wait for a full CPI to land in the capture buffers + while (true) { - time.push_back(current_time_us()); - // extract data from buffer - for (uint32_t i = 0; i < nSamples; i++) + buffer1->lock(); + buffer2->lock(); + if ((buffer1->get_length() > nSamples) && (buffer2->get_length() > nSamples)) { - x->push_back(buffer1->pop_front()); - y->push_back(buffer2->pop_front()); + break; } buffer1->unlock(); buffer2->unlock(); - timing_helper(timing_name, timing_time, time, "extract_buffer"); + // short delay to prevent tight looping + std::this_thread::sleep_for(std::chrono::milliseconds(1)); + } - // live retune: refresh wavelength and drop stale tracks on fc change - if (g_tuneState.fcChanged.exchange(false)) - { - fc = g_tuneState.currentFc.load(); - lambda = (double)Constants::c/fc; - tracker->set_lambda(lambda); - tracker->reset(); - } + slot->reset(); + slot->time.push_back(current_time_us()); + // extract data from buffer + for (uint32_t i = 0; i < nSamples; i++) + { + slot->x->push_back(buffer1->pop_front()); + slot->y->push_back(buffer2->pop_front()); + } + buffer1->unlock(); + buffer2->unlock(); + timing_helper(slot->timingName, slot->timingTime, slot->time, "extract_buffer"); + + // Latch a live retune against the CPI it applies to rather than acting + // on it here: the tracker it resets runs in the back stage. + slot->fcChanged = g_tuneState.fcChanged.exchange(false); + if (slot->fcChanged) + { + slot->fc = g_tuneState.currentFc.load(); + } - // spectrum - spectrumAnalyser->process(x); - timing_helper(timing_name, timing_time, time, "spectrum"); - - // clutter filter - if (isClutter) + // clutter filter + if (isClutter) + { + if (!filter->process(slot->x.get(), slot->y.get())) { - if (!filter->process(x, y)) + // Drop the CPI as the serial loop did, but hand the retune back so + // a failed filter cannot swallow it. + if (slot->fcChanged) { - continue; + g_tuneState.fcChanged.store(true); } - timing_helper(timing_name, timing_time, time, "clutter_filter"); - } - - // ambiguity process - map = ambiguity->process(x, y); - map->set_metrics(); - timing_helper(timing_name, timing_time, time, "ambiguity_processing"); - - // detection process - if (isDetection) - { - detection1 = cfarDetector1D->process(map); - detection2 = centroid->process(detection1.get()); - detection = interpolate->process(detection2.get(), map); - timing_helper(timing_name, timing_time, time, "detector"); + slot->x->clear(); + slot->y->clear(); + freeSlots.push(slot); + continue; } + timing_helper(slot->timingName, slot->timingTime, slot->time, "clutter_filter"); + } - // tracker process - if (isTracker) - { - track = tracker->process(detection.get(), time[0]/1000); - timing_helper(timing_name, timing_time, time, "tracker"); - } + filteredSlots.push(slot); + } + }); - // output IqData meta data - jsonIqData = x->to_json(time[0]/1000); - socket_iqdata.sendData(jsonIqData); + // Back stage: everything downstream of the clutter filter, in capture order. + // The tracker's state stays here, so it still sees every CPI exactly once and + // in sequence. + std::thread t3([&]{ + Map> *map; + std::unique_ptr detection; + std::unique_ptr detection1; + std::unique_ptr detection2; + std::unique_ptr track; + std::string mapJson, detectionJson, jsonTracker, jsonIqData, jsonTiming; + uint64_t previousCpiEnd = 0; - // output map data - mapJson = map->to_json_km(time[0]/1000, fs); - if (saveMap) - { - map->save(mapJson, saveMapPath); - } - socket_map.sendData(mapJson); + while (true) + { + blah2::CpiSlot *slot = filteredSlots.pop(); + std::vector &timing_name = slot->timingName; + std::vector &timing_time = slot->timingTime; + std::vector &time = slot->time; - // output detection data - if (isDetection) - { - detectionJson = detection->to_json_km(time[0]/1000, fs); - socket_detection.sendData(detectionJson); - } + // Idle time waiting on the front stage. Near zero means this stage is + // the bottleneck; large means the clutter filter is. + timing_helper(timing_name, timing_time, time, "pipeline_wait"); - // output tracker data - if (isTracker) - { - jsonTracker = track->to_json(time[0]/1000); - socket_track.sendData(jsonTracker); - } + // live retune: refresh wavelength and drop stale tracks on fc change + if (slot->fcChanged) + { + fc = slot->fc; + lambda = (double)Constants::c/fc; + tracker->set_lambda(lambda); + tracker->reset(); + } - // output radar data timer - timing_helper(timing_name, timing_time, time, "output_radar_data"); + // spectrum. Reads the reference channel, which the clutter filter only + // reads too, so running it after the filter instead of before leaves + // its output unchanged. It has to stay ahead of ambiguity, which drains + // the channel it reads. + spectrumAnalyser->process(slot->x.get()); + timing_helper(timing_name, timing_time, time, "spectrum"); - // cpi timer - time.push_back(current_time_us()); - double delta_ms = (double)(time.back()-time[0]) / 1000; - timing_name.push_back("cpi"); - timing_time.push_back(delta_ms); - if (verbose) - { - std::cout << "CPI time (ms): " << delta_ms << std::endl; - } + // ambiguity process + map = ambiguity->process(slot->x.get(), slot->y.get()); + map->set_metrics(); + timing_helper(timing_name, timing_time, time, "ambiguity_processing"); + + // detection process + if (isDetection) + { + detection1 = cfarDetector1D->process(map); + detection2 = centroid->process(detection1.get()); + detection = interpolate->process(detection2.get(), map); + timing_helper(timing_name, timing_time, time, "detector"); + } + + // tracker process + if (isTracker) + { + track = tracker->process(detection.get(), time[0]/1000); + timing_helper(timing_name, timing_time, time, "tracker"); + } - // output timing data - timing->update(time[0]/1000, timing_time, timing_name); - jsonTiming = timing->to_json(); - socket_timing.sendData(jsonTiming); - timing_time.clear(); - timing_name.clear(); + // output IqData meta data + jsonIqData = slot->x->to_json(time[0]/1000); + socket_iqdata.sendData(jsonIqData); - // output CPI timestamp for updating data - std::string t0_string = std::to_string(time[0]/1000); - socket_timestamp.sendData(t0_string); - time.clear(); + // output map data + mapJson = map->to_json_km(time[0]/1000, fs); + if (saveMap) + { + map->save(mapJson, saveMapPath); + } + socket_map.sendData(mapJson); + // output detection data + if (isDetection) + { + detectionJson = detection->to_json_km(time[0]/1000, fs); + socket_detection.sendData(detectionJson); } - else + + // output tracker data + if (isTracker) { - buffer1->unlock(); - buffer2->unlock(); - // short delay to prevent tight looping - std::this_thread::sleep_for(std::chrono::milliseconds(1)); + jsonTracker = track->to_json(time[0]/1000); + socket_track.sendData(jsonTracker); } + + // output radar data timer + timing_helper(timing_name, timing_time, time, "output_radar_data"); + + // cpi timer. Work done on this CPI across both stages, excluding the + // queue wait, so it stays the same quantity the serial loop reported. + double delta_ms = 0; + for (size_t i = 0; i < timing_time.size(); i++) + { + if (timing_name[i] != "pipeline_wait") + { + delta_ms += timing_time[i]; + } + } + timing_name.push_back("cpi"); + timing_time.push_back(delta_ms); + + // Wall clock between CPIs leaving the pipeline. With the stages + // overlapped this, not "cpi", is the throughput the radar achieves. + uint64_t cpiEnd = current_time_us(); + if (previousCpiEnd != 0) + { + timing_name.push_back("cpi_interval"); + timing_time.push_back((double)(cpiEnd - previousCpiEnd) / 1000); + } + previousCpiEnd = cpiEnd; + + if (verbose) + { + std::cout << "CPI time (ms): " << delta_ms << std::endl; + } + + // output timing data + timing->update(time[0]/1000, timing_time, timing_name); + jsonTiming = timing->to_json(); + socket_timing.sendData(jsonTiming); + + // output CPI timestamp for updating data + std::string t0_string = std::to_string(time[0]/1000); + socket_timestamp.sendData(t0_string); + + // Ambiguity leaves a partial batch behind; clearing gives the front + // stage the same empty queue the serial loop's eviction produced. + slot->x->clear(); + slot->y->clear(); + freeSlots.push(slot); } }); + t2.join(); + t3.join(); t1.join(); return 0; diff --git a/src/process/clutter/WienerHopf.cpp b/src/process/clutter/WienerHopf.cpp index cf06a513..233193f1 100644 --- a/src/process/clutter/WienerHopf.cpp +++ b/src/process/clutter/WienerHopf.cpp @@ -5,7 +5,8 @@ #include // constructor -WienerHopf::WienerHopf(int32_t _delayMin, int32_t _delayMax, uint32_t _nSamples) +WienerHopf::WienerHopf(int32_t _delayMin, int32_t _delayMax, uint32_t _nSamples, + int _plannerThreads) { // input delayMin = _delayMin; @@ -19,7 +20,7 @@ WienerHopf::WienerHopf(int32_t _delayMin, int32_t _delayMax, uint32_t _nSamples) // 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); + nFilter = blah2::fastestFftLength(uint64_t(nSamples) + nBins + 1, 0.02, _plannerThreads); // initialise data A = arma::cx_mat(nBins, nBins); diff --git a/src/process/clutter/WienerHopf.h b/src/process/clutter/WienerHopf.h index d7e6ccdb..35e906e1 100644 --- a/src/process/clutter/WienerHopf.h +++ b/src/process/clutter/WienerHopf.h @@ -10,6 +10,7 @@ #define WIENERHOPF_H #include "data/IqData.h" +#include "process/meta/FftLength.h" #include #include #include @@ -66,8 +67,12 @@ class WienerHopf /// @param delayMin Minimum clutter filter delay (bins). /// @param delayMax Maximum clutter filter delay (bins). /// @param nSamples Number of samples per CPI. + /// @param plannerThreads Threads the filter's plans will run with. The + /// quickest transform length depends on the thread count, so the length + /// search has to be told what the plans will actually use. /// @return The object. - WienerHopf(int32_t delayMin, int32_t delayMax, uint32_t nSamples); + WienerHopf(int32_t delayMin, int32_t delayMax, uint32_t nSamples, + int plannerThreads = blah2::kPlannerThreads); uint32_t filter_fft_length() const { return nFilter; } diff --git a/src/process/meta/FftLength.h b/src/process/meta/FftLength.h index 8d0eacb1..ace49ca3 100644 --- a/src/process/meta/FftLength.h +++ b/src/process/meta/FftLength.h @@ -12,6 +12,15 @@ namespace blah2 { // depends on this, so it belongs anywhere the choice is made or remembered. inline constexpr int kPlannerThreads = 4; +// The two processing stages run concurrently, so their plans have to share the +// four cores rather than each claiming all of them. Measured per-plan on a Pi 5 +// at the shipped geometry, the clutter transforms prefer 2 threads to 4 anyway +// (1e6 points: 118.2 ms against 146.8 ms per CPI, and 89.5 against 96.8 for the +// filter length), so the front stage gains from halving. The ambiguity batch +// transforms are the ones that pay, giving up about 21 ms. +inline constexpr int kFrontStageThreads = 2; +inline constexpr int kBackStageThreads = 2; + // 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) { diff --git a/src/process/utility/CpiPipeline.h b/src/process/utility/CpiPipeline.h new file mode 100644 index 00000000..e3e742e2 --- /dev/null +++ b/src/process/utility/CpiPipeline.h @@ -0,0 +1,107 @@ +/// @file CpiPipeline.h +/// @brief Slot and handoff machinery for processing CPIs in a pipeline. +/// @details The radar loop used to run every stage of a CPI back to back, so +/// throughput was the sum of the stages. Splitting it at the clutter filter +/// makes throughput the slowest stage instead. Output is unchanged: each stage +/// is deterministic, slots carry one CPI's state end to end, and everything +/// holding state across CPIs (the tracker) stays inside one serial stage, so +/// CPIs are still finished in the order they were captured. +/// @author Josh Poole + +#ifndef CPIPIPELINE_H +#define CPIPIPELINE_H + +#include "data/IqData.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace blah2 +{ + +/// @brief One CPI in flight between pipeline stages. +struct CpiSlot +{ + explicit CpiSlot(uint32_t nSamples) + : x(std::make_unique(nSamples)), y(std::make_unique(nSamples)) + { + } + + /// @brief Reference channel for this CPI. + std::unique_ptr x; + + /// @brief Surveillance channel for this CPI, overwritten by the clutter filter. + std::unique_ptr y; + + /// @brief Stage boundary timestamps (us). time[0] is the extract start, and + /// is the timestamp every output for this CPI is stamped with. + std::vector time; + + /// @brief Names of the stage timings accumulated so far. + std::vector timingName; + + /// @brief Stage timings (ms), in step with timingName. + std::vector timingTime; + + /// @brief True if a retune landed on this CPI. + /// @details Latched by the front stage at extract rather than acted on there, + /// because the tracker it resets lives in the back stage. Consuming the flag + /// in one stage and acting on it in the other would apply the reset to + /// whichever CPI happened to be in flight. + bool fcChanged = false; + + /// @brief Centre frequency to adopt, valid when fcChanged. + uint32_t fc = 0; + + /// @brief Ready the slot for a fresh CPI. + void reset() + { + time.clear(); + timingName.clear(); + timingTime.clear(); + fcChanged = false; + fc = 0; + } +}; + +/// @brief Blocking handoff between pipeline stages. +/// @details No capacity limit is needed: the number of slots in circulation is +/// fixed at startup, so a stage that runs ahead blocks on the free list. The +/// radar loop runs until the process is killed, so there is no close path. +template +class BlockingQueue +{ +public: + void push(T value) + { + { + std::lock_guard lock(mutex); + queue.push_back(std::move(value)); + } + condition.notify_one(); + } + + T pop() + { + std::unique_lock lock(mutex); + condition.wait(lock, [this] { return !queue.empty(); }); + T value = std::move(queue.front()); + queue.pop_front(); + return value; + } + +private: + std::mutex mutex; + std::condition_variable condition; + std::deque queue; +}; + +} + +#endif diff --git a/test/unit/process/clutter/TestStageOrder.cpp b/test/unit/process/clutter/TestStageOrder.cpp new file mode 100644 index 00000000..acd99f5f --- /dev/null +++ b/test/unit/process/clutter/TestStageOrder.cpp @@ -0,0 +1,206 @@ +// Splitting the radar loop at the clutter filter moves the spectrum stage from +// before the filter to after it, because the front stage is the slow half and +// the spectrum reads only the reference channel. That is output-neutral if and +// only if the clutter filter leaves the reference channel untouched, so this +// proves that on data rather than by reading the source. +// +// Also checks the slot recycle. Ambiguity drains all but a partial batch, and +// the serial loop pushed the next CPI on top and let IqData::push_back evict +// the remainder; the back stage now clears the slot instead. Those have to +// leave the same samples behind. + +#include "data/IqData.h" +#include "process/clutter/WienerHopf.h" +#include "process/spectrum/SpectrumAnalyser.h" + +#include +#include +#include +#include +#include +#include +#include + +using Complex = std::complex; + +static int failures = 0; + +static void require(bool condition, const char* message) +{ + std::printf("%-64s %s\n", message, condition ? "ok" : "FAIL"); + if (!condition) failures++; +} + +// A CPI with strong correlated clutter at several delays plus a weak echo, so +// the filter has real work to do and the solve is not degenerate. +static void fill(IqData* x, IqData* y, uint32_t n) +{ + std::vector ref(n); + for (uint32_t i = 0; i < n; i++) + { + ref[i] = {0.7 * std::sin(i * 0.013) + 0.3 * std::sin(i * 0.0007 + 1.1), + 0.7 * std::cos(i * 0.011) + 0.3 * std::cos(i * 0.0009 + 0.3)}; + } + for (uint32_t i = 0; i < n; i++) + { + Complex s = 0.9 * ref[i]; + if (i >= 5) s += 0.45 * ref[i - 5]; + if (i >= 40) s += 0.20 * ref[i - 40]; + if (i >= 137) s += 0.01 * ref[i - 137]; + s += Complex(1e-4 * std::sin(i * 0.37), 1e-4 * std::cos(i * 0.41)); + x->push_back(ref[i]); + y->push_back(s); + } +} + +static bool same(const std::deque& a, const std::deque& b) +{ + if (a.size() != b.size()) return false; + for (size_t i = 0; i < a.size(); i++) + if (a[i].real() != b[i].real() || a[i].imag() != b[i].imag()) return false; + return true; +} + +// The spectrum and frequency arrays, which is what the spectrum stage +// produces. IqData::to_json also emits min, max and mean, which are declared +// but never assigned anywhere in the codebase, so every IqData serialises +// whatever happened to be in those bytes. Two objects therefore differ there +// for reasons that have nothing to do with this test: comparing whole JSON +// passes on x86 and fails on aarch64 purely on allocator luck. +static std::string spectrum_of(const std::string& json) +{ + const size_t start = json.find("\"frequency\""); + return start == std::string::npos ? json : json.substr(start); +} + +int main() +{ + const uint32_t n = 200000; + const int32_t delayMin = -10; + const int32_t delayMax = 400; + const double spectrumBandwidth = 2000; + + // The clutter filter must not touch the reference channel. + { + IqData x(n), y(n); + fill(&x, &y, n); + const std::deque before = x.view_data(); + + WienerHopf filter(delayMin, delayMax, n, blah2::kFrontStageThreads); + require(filter.process(&x, &y), "clutter filter converged on the test signal"); + require(same(before, x.view_data()), + "reference channel bit-identical after the clutter filter"); + + // Confirm the filter actually cancelled, so the check above is not vacuous. + IqData xRaw(n), yRaw(n); + fill(&xRaw, &yRaw, n); + double residual = 0, signal = 0; + const auto& filtered = y.view_data(); + const auto& raw = yRaw.view_data(); + for (uint32_t i = n / 2; i < n; i++) + { + residual += std::norm(filtered[i]); + signal += std::norm(raw[i]); + } + const double cancelDb = 10.0 * std::log10(signal / residual); + std::printf(" clutter cancellation on the test signal: %.1f dB\n", cancelDb); + require(cancelDb > 20.0, "filter cancelled enough that the test is not vacuous"); + } + + // Spectrum output is the same whichever side of the filter it runs on. + { + IqData xa(n), ya(n); + fill(&xa, &ya, n); + SpectrumAnalyser serialOrder(n, spectrumBandwidth); + serialOrder.process(&xa); + const std::string jsonBefore = xa.to_json(1234567890ULL); + WienerHopf filterA(delayMin, delayMax, n, blah2::kFrontStageThreads); + filterA.process(&xa, &ya); + + IqData xb(n), yb(n); + fill(&xb, &yb, n); + WienerHopf filterB(delayMin, delayMax, n, blah2::kFrontStageThreads); + filterB.process(&xb, &yb); + SpectrumAnalyser pipelineOrder(n, spectrumBandwidth); + pipelineOrder.process(&xb); + const std::string jsonAfter = xb.to_json(1234567890ULL); + + require(!spectrum_of(jsonBefore).empty(), "spectrum was actually produced"); + require(spectrum_of(jsonBefore) == spectrum_of(jsonAfter), + "spectrum byte-identical whichever side of the filter it runs"); + require(same(ya.view_data(), yb.view_data()), + "filtered surveillance channel identical under the reordering"); + } + + // Recycling a slot by clearing it matches the serial refill-with-eviction. + { + const uint32_t leftover = 78; // what ambiguity's partial batch leaves + IqData serial(n), pipelined(n); + for (uint32_t i = 0; i < n; i++) + { + const Complex v = {(double)i, -(double)i}; + serial.push_back(v); + pipelined.push_back(v); + } + for (uint32_t i = 0; i < n - leftover; i++) + { + serial.pop_front(); + pipelined.pop_front(); + } + require(serial.get_length() == leftover, "leftover partial batch reproduced"); + + pipelined.clear(); + for (uint32_t i = 0; i < n; i++) + { + const Complex v = {(double)i * 3.5, (double)i * -0.25}; + serial.push_back(v); + pipelined.push_back(v); + } + require(same(serial.view_data(), pipelined.view_data()), + "clear-then-refill matches the serial refill-with-eviction"); + } + + // The two stages plan at 2 threads each rather than both claiming 4, so the + // transforms are not the ones the serial build planned. FFTW parallelises + // across a Cooley-Tukey factor, so a different thread count sums in a + // different order and the numbers move. Bound how far: the FFT length work + // already accepted 6e-14 relative, about 264 dB below signal, so this has to + // be no worse. + { + fftw_init_threads(); + const uint32_t m = 200000; + std::vector seed(m); + for (uint32_t i = 0; i < m; i++) + seed[i] = {std::sin(i * 0.001), std::cos(i * 0.0007)}; + + auto transform = [&](int threads) { + std::vector buffer(seed); + fftw_plan_with_nthreads(threads); + fftw_plan plan = fftw_plan_dft_1d( + m, reinterpret_cast(buffer.data()), + reinterpret_cast(buffer.data()), FFTW_FORWARD, FFTW_ESTIMATE); + fftw_execute(plan); + fftw_destroy_plan(plan); + return buffer; + }; + + const std::vector serial = transform(blah2::kPlannerThreads); + const std::vector pipelined = transform(blah2::kBackStageThreads); + + double peak = 0, worst = 0; + for (uint32_t i = 0; i < m; i++) + { + peak = std::max(peak, std::abs(serial[i])); + worst = std::max(worst, std::abs(serial[i] - pipelined[i])); + } + const double relative = worst / peak; + std::printf(" planner %d threads vs %d: %.3g of peak (%.0f dB below)\n", + blah2::kPlannerThreads, blah2::kBackStageThreads, relative, + 20.0 * std::log10(peak / worst)); + require(relative < 6e-14, + "per-stage thread count moves an FFT less than the accepted 6e-14"); + } + + std::printf("\n%s\n", failures ? "FAILURES" : "all checks passed"); + return failures ? 1 : 0; +} diff --git a/test/unit/process/utility/TestCpiPipeline.cpp b/test/unit/process/utility/TestCpiPipeline.cpp new file mode 100644 index 00000000..3cc6d23d --- /dev/null +++ b/test/unit/process/utility/TestCpiPipeline.cpp @@ -0,0 +1,164 @@ +// The radar loop processes a CPI in two overlapped stages rather than one +// serial run, so throughput is the slowest stage instead of their sum. That is +// only output-neutral if the handoff preserves what the serial loop guaranteed: +// every CPI reaches the back stage exactly once, in capture order, with its own +// timing and retune state, and dropped CPIs return their slot without losing a +// pending retune. This exercises those properties with the stages replaced by +// sleeps of the measured owl-ded9 durations. + +#include "process/utility/CpiPipeline.h" + +#include +#include +#include +#include +#include +#include +#include + +using blah2::BlockingQueue; +using blah2::CpiSlot; + +static int failures = 0; + +static void require(bool condition, const char* message) +{ + std::printf("%-62s %s\n", message, condition ? "ok" : "FAIL"); + if (!condition) failures++; +} + +static uint64_t now_us() +{ + return std::chrono::duration_cast( + std::chrono::steady_clock::now().time_since_epoch()) + .count(); +} + +static void work(double ms) +{ + std::this_thread::sleep_for(std::chrono::microseconds((long)(ms * 1000))); +} + +// Measured stage times on owl-ded9, scaled so the test runs in seconds. Only +// the ratio between the two stages matters to what is being checked. +static constexpr double kScale = 0.02; +static constexpr double kExtract = 36.0 * kScale; +static constexpr double kClutter = 461.3 * kScale; +static constexpr double kBackStage = (68.5 + 240.4 + 3.3 + 22.0) * kScale; + +static constexpr int kCpis = 40; +static constexpr uint64_t kFailingCpi = 7; + +int main() +{ + constexpr size_t kPipelineDepth = 2; + std::vector> slots; + BlockingQueue freeSlots; + BlockingQueue filteredSlots; + for (size_t i = 0; i < kPipelineDepth; i++) + { + slots.push_back(std::make_unique(1)); + freeSlots.push(slots.back().get()); + } + + std::atomic fcChanged{false}; + std::atomic stop{false}; + std::vector backOrder; + std::vector retunesApplied; + uint64_t nextCpiId = 0; + uint64_t t0 = now_us(); + + // Front stage: extract, then clutter filter. CPI kFailingCpi fails the + // filter, with a retune raised just before it. + std::thread front([&] { + while (!stop.load()) + { + CpiSlot* slot = freeSlots.pop(); + if (stop.load()) break; + slot->reset(); + slot->time.push_back(now_us()); + uint64_t id = nextCpiId++; + work(kExtract); + slot->timingName.push_back("extract_buffer"); + slot->timingTime.push_back(kExtract); + + if (id == kFailingCpi) fcChanged.store(true); + slot->fcChanged = fcChanged.exchange(false); + if (slot->fcChanged) slot->fc = 500000000u + (uint32_t)id; + + work(kClutter); + if (id == kFailingCpi) + { + if (slot->fcChanged) fcChanged.store(true); + slot->x->clear(); + slot->y->clear(); + freeSlots.push(slot); + continue; + } + slot->timingName.push_back("clutter_filter"); + slot->timingTime.push_back(kClutter); + slot->time.push_back(id); // stands in for the CPI's payload + filteredSlots.push(slot); + if ((int)id >= kCpis) stop.store(true); + } + }); + + // Back stage: everything downstream, in capture order, holding the state the + // tracker would. + std::thread back([&] { + uint64_t lastSeen = 0; + bool first = true; + for (int seen = 0; seen < kCpis - 1; seen++) + { + CpiSlot* slot = filteredSlots.pop(); + uint64_t id = slot->time[1]; + backOrder.push_back(id); + if (slot->fcChanged) retunesApplied.push_back(id); + if (!first && id <= lastSeen) failures++; + lastSeen = id; + first = false; + work(kBackStage); + slot->x->clear(); + slot->y->clear(); + freeSlots.push(slot); + } + }); + + back.join(); + uint64_t elapsed = now_us() - t0; + stop.store(true); + freeSlots.push(slots[0].get()); // let the front stage observe stop and exit + front.join(); + + const double perCpiMs = (double)elapsed / 1000.0 / (kCpis - 1); + const double serialMs = kExtract + kClutter + kBackStage; + const double frontMs = kExtract + kClutter; + + std::printf("throughput %.2f ms/CPI, serial would be %.2f, front stage %.2f\n\n", + perCpiMs, serialMs, frontMs); + + require(backOrder.size() == (size_t)(kCpis - 1), "back stage saw every CPI exactly once"); + + bool ordered = true; + for (size_t i = 1; i < backOrder.size(); i++) + if (backOrder[i] <= backOrder[i - 1]) ordered = false; + require(ordered, "CPIs reached the back stage in capture order"); + + bool contiguous = true; + for (size_t i = 0; i < backOrder.size(); i++) + { + const uint64_t expect = i < kFailingCpi ? i : i + 1; + if (backOrder[i] != expect) contiguous = false; + } + require(contiguous, "only the dropped CPI is missing, no others"); + + require(retunesApplied.size() == 1 && retunesApplied[0] == kFailingCpi + 1, + "retune survived the dropped CPI and landed on the next"); + + require(perCpiMs < serialMs * 0.85, "throughput beats the serial sum of stages"); + require(perCpiMs > frontMs * 0.85 && perCpiMs < frontMs * 1.35, + "throughput tracks the slow stage, not the sum"); + + std::printf("\n%s\n", failures ? "FAILURES" : "all checks passed"); + return failures ? 1 : 0; +}