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; +}