From 5f24889bcf609b44e94689df004cac659444fa4c Mon Sep 17 00:00:00 2001 From: Daniel Gustainis Date: Fri, 15 Dec 2023 16:55:00 +1030 Subject: [PATCH] Added option to round range FFT lengths to the next Hamming number --- src/process/ambiguity/Ambiguity.cpp | 72 ++++++++++++++++++++++-- src/process/ambiguity/Ambiguity.h | 22 +++++++- src/process/ambiguity/test_ambiguity.cpp | 50 +++++++++++++--- 3 files changed, 130 insertions(+), 14 deletions(-) diff --git a/src/process/ambiguity/Ambiguity.cpp b/src/process/ambiguity/Ambiguity.cpp index b2a42eb..05ddf3e 100644 --- a/src/process/ambiguity/Ambiguity.cpp +++ b/src/process/ambiguity/Ambiguity.cpp @@ -5,9 +5,10 @@ #include #include #include +#include // constructor -Ambiguity::Ambiguity(int32_t delayMin, int32_t delayMax, int32_t dopplerMin, int32_t dopplerMax, uint32_t fs, uint32_t n) +Ambiguity::Ambiguity(int32_t delayMin, int32_t delayMax, int32_t dopplerMin, int32_t dopplerMax, uint32_t fs, uint32_t n, bool roundHamming) : delayMin_{delayMin} , delayMax_{delayMax} , dopplerMin_{dopplerMin} @@ -20,7 +21,6 @@ Ambiguity::Ambiguity(int32_t delayMin, int32_t delayMax, int32_t dopplerMin, int // doppler calculations std::deque doppler; double resolutionDoppler = 1.0 / (static_cast(n) / static_cast(fs)); - doppler_res_ = resolutionDoppler; doppler.push_back(dopplerMiddle_); int i = 1; while (dopplerMiddle_ + (i * resolutionDoppler) <= dopplerMax) @@ -55,7 +55,10 @@ Ambiguity::Ambiguity(int32_t delayMin, int32_t delayMax, int32_t dopplerMin, int } // other setup - nfft_ = (2 * nCorr_) - 1; + nfft_ = 2 * nCorr_ - 1; + if (roundHamming) { + nfft_ = next_hamming(nfft_); + } dataCorr_.resize(2 * nDelayBins_ + 1); // compute FFTW plans in constructor @@ -84,6 +87,10 @@ Ambiguity::~Ambiguity() Map> *Ambiguity::process(IqData *x, IqData *y) { + using Timer = std::chrono::steady_clock; + auto t0{Timer::now()}; + Timer::duration range_fft_dur{}; + // shift reference if not 0 centered if (dopplerMiddle_ != 0) { @@ -109,9 +116,10 @@ Map> *Ambiguity::process(IqData *x, IqData *y) dataYi_[j] = {0, 0}; } + auto t1{Timer::now()}; fftw_execute(fftXi_); fftw_execute(fftYi_); - + range_fft_dur += Timer::now() - t1; // compute correlation for (int j = 0; j < nfft_; j++) @@ -119,7 +127,9 @@ Map> *Ambiguity::process(IqData *x, IqData *y) dataZi_[j] = (dataYi_[j] * std::conj(dataXi_[j])) / (double)nfft_; } + t1 = Timer::now(); fftw_execute(fftZi_); + range_fft_dur += Timer::now() - t1; // extract center of corr for (int j = 0; j < nDelayBins_; j++) @@ -142,6 +152,7 @@ Map> *Ambiguity::process(IqData *x, IqData *y) } // doppler processing + auto t1{Timer::now()}; for (int i = 0; i < nDelayBins_; i++) { delayProfile_ = map_->get_col(i); @@ -161,5 +172,58 @@ Map> *Ambiguity::process(IqData *x, IqData *y) map_->set_col(i, corr_); } + auto to_ms = [] (const Timer::duration& dur) { + return std::chrono::duration_cast>(dur).count(); + }; + + latest_performance_.process_time_ms = to_ms(Timer::now() - t0); + latest_performance_.doppler_fft_time_ms = to_ms(Timer::now() - t1); + latest_performance_.range_fft_time_ms = to_ms(range_fft_dur); + return map_.get(); } + +/** + * @brief Hamming number generator + * + * @author Nigel Galloway + * @cite https://rosettacode.org/wiki/Hamming_numbers + * @todo Can this be done with constexpr??? + */ +class HammingGenerator { +private: + std::vector _H, _hp, _hv, _x; +public: + bool operator!=(const HammingGenerator &other) const { return true; } + HammingGenerator begin() const { return *this; } + HammingGenerator end() const { return *this; } + unsigned int operator*() const { return _x.back(); } + HammingGenerator(const std::vector &pfs) : _H(pfs), _hp(pfs.size(), 0), _hv({pfs}), _x({1}) {} + const HammingGenerator &operator++() + { + for (int i = 0; i < _H.size(); i++) + for (; _hv[i] <= _x.back(); _hv[i] = _x[++_hp[i]] * _H[i]) + ; + _x.push_back(_hv[0]); + for (int i = 1; i < _H.size(); i++) + if (_hv[i] < _x.back()) + _x.back() = _hv[i]; + return *this; + } +}; + + +uint32_t next_hamming(uint32_t value) { + for (auto i : HammingGenerator({2,3,5})) { + if (i > value) { + return i; + } + } + return 0; +} + +std::ostream& operator<<(std::ostream& str, const Ambiguity::PerformanceStats& stats) { + return str << "Total time: " << stats.process_time_ms << "ms\n" << + "Range FFT time: " << stats.range_fft_time_ms << "ms\n" << + "Doppler FFT time: " << stats.doppler_fft_time_ms << "ms"; +} \ No newline at end of file diff --git a/src/process/ambiguity/Ambiguity.h b/src/process/ambiguity/Ambiguity.h index 9137f59..5c0dc06 100644 --- a/src/process/ambiguity/Ambiguity.h +++ b/src/process/ambiguity/Ambiguity.h @@ -21,6 +21,13 @@ public: using Complex = std::complex; + struct PerformanceStats { + double process_time_ms{0}; + double range_fft_time_ms{0}; + double doppler_fft_time_ms{0}; + }; + + /// @brief Constructor. /// @param delayMin Minimum delay (bins). /// @param delayMax Maximum delay (bins). @@ -28,8 +35,9 @@ public: /// @param dopplerMax Maximum Doppler (Hz). /// @param fs Sampling frequency (Hz). /// @param n Number of samples. + /// @param roundHamming Round the correlation FFT length to a Hamming number for performance /// @return The object. - Ambiguity(int32_t delayMin, int32_t delayMax, int32_t dopplerMin, int32_t dopplerMax, uint32_t fs, uint32_t n); + Ambiguity(int32_t delayMin, int32_t delayMax, int32_t dopplerMin, int32_t dopplerMax, uint32_t fs, uint32_t n, bool roundHamming = false); /// @brief Destructor. /// @return Void. @@ -53,7 +61,7 @@ public: uint32_t fft_bin_count() const { return nfft_; } - double doppler_res_; + PerformanceStats get_latest_performance() const { return latest_performance_; } private: /// @brief Minimum delay (bins). int32_t delayMin_; @@ -115,4 +123,12 @@ private: /// @brief Map to store result. std::unique_ptr> map_; -}; \ No newline at end of file + PerformanceStats latest_performance_; +}; + +/// @brief Calculate the next 5-smooth Hamming Number larger than value +/// @param value Value to round +/// @return value rounded to Hamming number +uint32_t next_hamming(uint32_t value); + +std::ostream& operator<<(std::ostream& str, const Ambiguity::PerformanceStats& stats); \ No newline at end of file diff --git a/src/process/ambiguity/test_ambiguity.cpp b/src/process/ambiguity/test_ambiguity.cpp index 60cbd71..9aa237b 100644 --- a/src/process/ambiguity/test_ambiguity.cpp +++ b/src/process/ambiguity/test_ambiguity.cpp @@ -3,11 +3,10 @@ #include "Ambiguity.h" #include +#include std::random_device g_rd; -using namespace Catch::literals; - // Have to use out ref parameter because there's no copy/move ctors void random_iq(IqData& iq_data) { std::mt19937 gen(g_rd()); @@ -72,8 +71,8 @@ TEST_CASE("Constructor", "[constructor]") CHECK(ambiguity.fft_bin_count() == 6643); } -// Make sure process produces an output -TEST_CASE("Process_Simple", "[process]") +// Make sure the constructor is calculating the parameters correctly with rounded FFT length +TEST_CASE("Constructor_Round", "[constructor]") { int32_t delay_min{-10}; int32_t delay_max{300}; @@ -84,7 +83,30 @@ TEST_CASE("Process_Simple", "[process]") float cpi_s{0.5}; uint32_t n_samples = cpi_s * fs; // narrow on purpose - Ambiguity ambiguity(delay_min,delay_max,doppler_min,doppler_max,fs,n_samples); + Ambiguity ambiguity(delay_min,delay_max,doppler_min,doppler_max,fs,n_samples,true); + + CHECK_THAT(ambiguity.cpi_length_seconds(), Catch::Matchers::WithinAbs(cpi_s, 0.02)); + CHECK(ambiguity.doppler_middle() == 0); + CHECK(ambiguity.corr_samples_per_pulse() == 3322); + CHECK(ambiguity.delay_bin_count() == delay_max + std::abs(delay_min) + 1); + CHECK(ambiguity.doppler_bin_count() == 301); + CHECK(ambiguity.fft_bin_count() == 6750); +} + +TEST_CASE("Process_Simple", "[process]") +{ + auto round_hamming = GENERATE(true, false); + + int32_t delay_min{-10}; + int32_t delay_max{300}; + int32_t doppler_min{-300}; + int32_t doppler_max{300}; + + uint32_t fs{2'000'000}; + float cpi_s{0.5}; + uint32_t n_samples = cpi_s * fs; // narrow on purpose + + Ambiguity ambiguity(delay_min,delay_max,doppler_min,doppler_max,fs,n_samples, round_hamming); IqData x{n_samples}; IqData y{n_samples}; @@ -95,11 +117,15 @@ TEST_CASE("Process_Simple", "[process]") map->set_metrics(); CHECK(map->maxPower > 0.0); CHECK(map->noisePower > 0.0); + + std::cout << "Process_Simple with" << (round_hamming ? " hamming\n" : "out hamming\n") + << ambiguity.get_latest_performance() << "\n-----------" << std::endl; } -// Sanity check that we're getting numbers close to the baseline ambiguity processing function. TEST_CASE("Process_File", "[process]") { + auto round_hamming = GENERATE(true, false); + int32_t delay_min{-10}; int32_t delay_max{300}; int32_t doppler_min{-300}; @@ -109,7 +135,7 @@ TEST_CASE("Process_File", "[process]") float cpi_s{0.5}; uint32_t n_samples = cpi_s * fs; // narrow on purpose - Ambiguity ambiguity(delay_min,delay_max,doppler_min,doppler_max,fs,n_samples); + Ambiguity ambiguity(delay_min,delay_max,doppler_min,doppler_max,fs,n_samples, round_hamming); IqData x{n_samples}; IqData y{n_samples}; @@ -120,4 +146,14 @@ TEST_CASE("Process_File", "[process]") map->set_metrics(); CHECK_THAT(map->maxPower, Catch::Matchers::WithinAbs(30.2816, 0.001)); CHECK_THAT(map->noisePower, Catch::Matchers::WithinAbs(76.918, 0.001)); + + std::cout << "Process_File with" << (round_hamming ? " hamming\n" : "out hamming\n") + << ambiguity.get_latest_performance() << "\n-----------" << std::endl; +} + +TEST_CASE("Next_Hamming", "[hamming]") +{ + CHECK(next_hamming(104) == 108); + CHECK(next_hamming(3322) == 3375); + CHECK(next_hamming(19043) == 19200); } \ No newline at end of file