Added option to round range FFT lengths to the next Hamming number

src/process/ambiguity/Ambiguity.cpp

@@ -5,9 +5,10 @@
 #include <vector>
 #include <numeric>
 #include <math.h>
+#include <chrono>
 
 // 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<double> doppler;
   double resolutionDoppler = 1.0 / (static_cast<double>(n) / static_cast<double>(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<std::complex<double>> *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<std::complex<double>> *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<std::complex<double>> *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<std::complex<double>> *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<std::complex<double>> *Ambiguity::process(IqData *x, IqData *y)
     map_->set_col(i, corr_);
   }
 
+  auto to_ms = [] (const Timer::duration& dur) {
+    return std::chrono::duration_cast<std::chrono::duration<double, std::milli>>(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<unsigned int> _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<unsigned int> &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";
+}

src/process/ambiguity/Ambiguity.h

@@ -21,6 +21,13 @@ public:
 
   using Complex = std::complex<double>;
 
+  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<Complex>> map_;
 
-};
+  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);

src/process/ambiguity/test_ambiguity.cpp

@@ -3,11 +3,10 @@
 
 #include "Ambiguity.h"
 #include <random>
+#include <iostream>
 
 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);
 }