Added ambiguity processing unit test and refactored for code quality and memory management

CMakeLists.txt

@@ -36,6 +36,8 @@ add_library(rapidjson ${PROJECT_LIB_DIR}/rapidjson-1.1.0/)
 add_library(sdrplay /usr/local/include/sdrplay_api.h)
 add_library(asio ${PROJECT_LIB_DIR}/asio-1.26.0/asio.hpp)
 add_library(httplib ${PROJECT_LIB_DIR}/cpp-httplib-0.12.2/httplib.h)
+add_library(catch2 ${PROJECT_LIB_DIR}/catch2/catch_amalgamated.cpp)
+target_include_directories(catch2 PUBLIC ${PROJECT_LIB_DIR}/catch2)
 
 include_directories("${PROJECT_LIB_DIR}/rapidjson-1.1.0/")
 set_target_properties(rapidjson PROPERTIES LINKER_LANGUAGE CXX)
@@ -71,3 +73,10 @@ include_directories("${PROJECT_SOURCE_DIR}/process/detection/")
 include_directories("${PROJECT_SOURCE_DIR}/process/spectrum/")
 include_directories("${PROJECT_SOURCE_DIR}/data/")
 include_directories("${PROJECT_SOURCE_DIR}/data/meta/")
+
+add_executable(test_ambiguity
+  ${PROJECT_SOURCE_DIR}/data/IqData.cpp
+  ${PROJECT_SOURCE_DIR}/data/Map.cpp
+  ${PROJECT_SOURCE_DIR}/process/ambiguity/Ambiguity.cpp
+  ${PROJECT_SOURCE_DIR}/process/ambiguity/test_ambiguity.cpp)
+target_link_libraries(test_ambiguity catch2 fftw3 fftw3_threads)

src/process/ambiguity/Ambiguity.cpp

@@ -3,174 +3,163 @@
 #include <iostream>
 #include <deque>
 #include <vector>
+#include <numeric>
 #include <math.h>
 
 // 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)
+  : delayMin_{delayMin}
+  , delayMax_{delayMax}
+  , dopplerMin_{dopplerMin}
+  , dopplerMax_{dopplerMax}
+  , fs_{fs}
+  , nSamples_{n}
+  , nDelayBins_{static_cast<uint16_t>(delayMax - delayMin + 1)} // If delayMin > delayMax = trouble, what's the exception policy?
+  , dopplerMiddle_{(dopplerMin_ + dopplerMax_) / 2.0}
 {
-  // input
-  delayMin = _delayMin;
-  delayMax = _delayMax;
-  dopplerMin = _dopplerMin;
-  dopplerMax = _dopplerMax;
-  fs = _fs;
-  n = _n;
-
-  // delay calculations
-  std::deque<int> delay;
-  nDelayBins = delayMax - delayMin + 1;
-  for (int i = 0; i < nDelayBins; i++)
-  {
-    delay.push_back(delayMin + i);
-  }
-
   // doppler calculations
   std::deque<double> doppler;
-  double resolutionDoppler = (double)1 / ((double)n / (double)fs);
-  dopplerMiddle = (dopplerMin + dopplerMax) / 2;
-  doppler.push_back(dopplerMiddle);
+  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)
+  while (dopplerMiddle_ + (i * resolutionDoppler) <= dopplerMax)
   {
-    doppler.push_back(dopplerMiddle + (i * resolutionDoppler));
-    doppler.push_front(dopplerMiddle - (i * resolutionDoppler));
+    doppler.push_back(dopplerMiddle_ + (i * resolutionDoppler));
+    doppler.push_front(dopplerMiddle_ - (i * resolutionDoppler));
     i++;
   }
-  nDopplerBins = doppler.size();
+  nDopplerBins_ = doppler.size();
 
   // batches constants
-  nCorr = (int)(n / nDopplerBins);
-  cpi = ((double)nCorr * (double)nDopplerBins) / fs;
+  nCorr_ = n / nDopplerBins_;
+  cpi_ = (static_cast<double>(nCorr_) * nDopplerBins_) / fs;
 
   // update doppler bins to true cpi time
-  resolutionDoppler = 1 / cpi;
-  doppler.clear();
-  doppler.push_back(dopplerMiddle);
+  resolutionDoppler = 1.0 / cpi_;
+
+  // create ambiguity map
+  map_ = std::make_unique<Map<Complex>>(nDopplerBins_, nDelayBins_);
+
+  // delay calculations
+  map_->delay.resize(nDelayBins_);
+  std::iota(map_->delay.begin(), map_->delay.end(), delayMin_);
+
+  map_->doppler.push_front(dopplerMiddle_);
   i = 1;
-  while (doppler.size() < nDopplerBins)
+  while (map_->doppler.size() < nDopplerBins_)
   {
-    doppler.push_back(dopplerMiddle + (i * resolutionDoppler));
-    doppler.push_front(dopplerMiddle - (i * resolutionDoppler));
+    map_->doppler.push_back(dopplerMiddle_ + (i * resolutionDoppler));
+    map_->doppler.push_front(dopplerMiddle_ - (i * resolutionDoppler));
     i++;
   }
 
   // other setup
-  nfft = (2 * nCorr) - 1;
-  dataCorr = new std::complex<double>[2 * nDelayBins + 1];
+  nfft_ = (2 * nCorr_) - 1;
+  dataCorr_.resize(2 * nDelayBins_ + 1);
 
   // compute FFTW plans in constructor
-  dataXi = new std::complex<double>[nfft];
-  dataYi = new std::complex<double>[nfft];
-  dataZi = new std::complex<double>[nfft];
-  dataDoppler = new std::complex<double>[nfft];
-  fftXi = fftw_plan_dft_1d(nfft, reinterpret_cast<fftw_complex *>(dataXi),
-                           reinterpret_cast<fftw_complex *>(dataXi), FFTW_FORWARD, FFTW_ESTIMATE);
-  fftYi = fftw_plan_dft_1d(nfft, reinterpret_cast<fftw_complex *>(dataYi),
-                           reinterpret_cast<fftw_complex *>(dataYi), FFTW_FORWARD, FFTW_ESTIMATE);
-  fftZi = fftw_plan_dft_1d(nfft, reinterpret_cast<fftw_complex *>(dataZi),
-                           reinterpret_cast<fftw_complex *>(dataZi), FFTW_BACKWARD, FFTW_ESTIMATE);
-  fftDoppler = fftw_plan_dft_1d(nDopplerBins, reinterpret_cast<fftw_complex *>(dataDoppler),
-                                reinterpret_cast<fftw_complex *>(dataDoppler), FFTW_FORWARD, FFTW_ESTIMATE);
+  dataXi_.resize(nfft_);
+  dataYi_.resize(nfft_);
+  dataZi_.resize(nfft_);
+  dataDoppler_.resize(nfft_);
+  fftXi_ = fftw_plan_dft_1d(nfft_, reinterpret_cast<fftw_complex *>(dataXi_.data()),
+                           reinterpret_cast<fftw_complex *>(dataXi_.data()), FFTW_FORWARD, FFTW_ESTIMATE);
+  fftYi_ = fftw_plan_dft_1d(nfft_, reinterpret_cast<fftw_complex *>(dataYi_.data()),
+                           reinterpret_cast<fftw_complex *>(dataYi_.data()), FFTW_FORWARD, FFTW_ESTIMATE);
+  fftZi_ = fftw_plan_dft_1d(nfft_, reinterpret_cast<fftw_complex *>(dataZi_.data()),
+                           reinterpret_cast<fftw_complex *>(dataZi_.data()), FFTW_BACKWARD, FFTW_ESTIMATE);
+  fftDoppler_ = fftw_plan_dft_1d(nDopplerBins_, reinterpret_cast<fftw_complex *>(dataDoppler_.data()),
+                                reinterpret_cast<fftw_complex *>(dataDoppler_.data()), FFTW_FORWARD, FFTW_ESTIMATE);
 
-  // create ambiguity map
-  map = new Map<std::complex<double>>(nDopplerBins, nDelayBins);
-
-  // store map parameters
-  for (int i = 0; i < nDelayBins; i++)
-  {
-    map->delay.push_back(delay[i]);
-  }
-  for (int i = 0; i < nDopplerBins; i++)
-  {
-    map->doppler.push_back(doppler[i]);
-  }
 }
 
 Ambiguity::~Ambiguity()
 {
-  fftw_destroy_plan(fftXi);
-  fftw_destroy_plan(fftYi);
-  fftw_destroy_plan(fftZi);
-  fftw_destroy_plan(fftDoppler);
+  fftw_destroy_plan(fftXi_);
+  fftw_destroy_plan(fftYi_);
+  fftw_destroy_plan(fftZi_);
+  fftw_destroy_plan(fftDoppler_);
 }
 
 Map<std::complex<double>> *Ambiguity::process(IqData *x, IqData *y)
 {
   // shift reference if not 0 centered
-  if (dopplerMiddle != 0)
+  if (dopplerMiddle_ != 0)
   {
     std::complex<double> j = {0, 1};
     for (int i = 0; i < x->get_length(); i++)
     {
-      x->push_back(x->pop_front() * std::exp(1.0 * j * 2.0 * M_PI * dopplerMiddle * ((double)i / fs)));
+      x->push_back(x->pop_front() * std::exp(1.0 * j * 2.0 * M_PI * dopplerMiddle_ * ((double)i / fs_)));
     }
   }
 
   // range processing
-  for (int i = 0; i < nDopplerBins; i++)
+  for (int i = 0; i < nDopplerBins_; i++)
   {
-    for (int j = 0; j < nCorr; j++)
+    for (int j = 0; j < nCorr_; j++)
     {
-      dataXi[j] = x->pop_front();
-      dataYi[j] = y->pop_front();
+      dataXi_[j] = x->pop_front();
+      dataYi_[j] = y->pop_front();
     }
 
-    for (int j = nCorr; j < nfft; j++)
+    for (int j = nCorr_; j < nfft_; j++)
     {
-      dataXi[j] = {0, 0};
-      dataYi[j] = {0, 0};
+      dataXi_[j] = {0, 0};
+      dataYi_[j] = {0, 0};
     }
 
-    fftw_execute(fftXi);
-    fftw_execute(fftYi);
+    fftw_execute(fftXi_);
+    fftw_execute(fftYi_);
+
 
     // compute correlation
-    for (int j = 0; j < nfft; j++)
+    for (int j = 0; j < nfft_; j++)
     {
-      dataZi[j] = (dataYi[j] * std::conj(dataXi[j])) / (double)nfft;
+      dataZi_[j] = (dataYi_[j] * std::conj(dataXi_[j])) / (double)nfft_;
     }
 
-    fftw_execute(fftZi);
+    fftw_execute(fftZi_);
 
     // extract center of corr
-    for (int j = 0; j < nDelayBins; j++)
+    for (int j = 0; j < nDelayBins_; j++)
     {
-      dataCorr[j] = dataZi[nfft - nDelayBins + j];
+      dataCorr_[j] = dataZi_[nfft_ - nDelayBins_ + j];
     }
-    for (int j = 0; j < nDelayBins + 1; j++)
+    for (int j = 0; j < nDelayBins_ + 1; j++)
     {
-      dataCorr[j + nDelayBins] = dataZi[j];
+      dataCorr_[j + nDelayBins_] = dataZi_[j];
     }
 
     // cast from std::complex to std::vector
-    corr.clear();
-    for (int j = 0; j < nDelayBins; j++)
+    corr_.clear();
+    for (int j = 0; j < nDelayBins_; j++)
     {
-      corr.push_back(dataCorr[nDelayBins + delayMin + j - 1 + 1]);
+      corr_.push_back(dataCorr_[nDelayBins_ + delayMin_ + j - 1 + 1]);
     }
 
-    map->set_row(i, corr);
+    map_->set_row(i, corr_);
   }
 
   // doppler processing
-  for (int i = 0; i < nDelayBins; i++)
+  for (int i = 0; i < nDelayBins_; i++)
   {
-    delayProfile = map->get_col(i);
-    for (int j = 0; j < nDopplerBins; j++)
+    delayProfile_ = map_->get_col(i);
+    for (int j = 0; j < nDopplerBins_; j++)
     {
-      dataDoppler[j] = {delayProfile[j].real(), delayProfile[j].imag()};
+      dataDoppler_[j] = {delayProfile_[j].real(), delayProfile_[j].imag()};
     }
 
-    fftw_execute(fftDoppler);
+    fftw_execute(fftDoppler_);
 
-    corr.clear();
-    for (int j = 0; j < nDopplerBins; j++)
+    corr_.clear();
+    for (int j = 0; j < nDopplerBins_; j++)
     {
-      corr.push_back(dataDoppler[(j + int(nDopplerBins / 2) + 1) % nDopplerBins]);
+      corr_.push_back(dataDoppler_[(j + int(nDopplerBins_ / 2) + 1) % nDopplerBins_]);
     }
 
-    map->set_col(i, corr);
+    map_->set_col(i, corr_);
   }
 
-  return map;
+  return map_.get();
 }

src/process/ambiguity/Ambiguity.h

@@ -6,72 +6,21 @@
 /// @author 30hours
 /// @todo Ambiguity maps are still offset by 1 bin.
 
-#ifndef AMBIGUITY_H
-#define AMBIGUITY_H
+#pragma once
 
 #include <IqData.h>
 #include <Map.h>
 #include <stdint.h>
 #include <fftw3.h>
+#include <memory>
 
 class Ambiguity
 {
-private:
-  /// @brief Minimum delay (bins).
-  int32_t delayMin;
-
-  /// @brief Maximum delay (bins).
-  int32_t delayMax;
-
-  /// @brief Minimum Doppler (Hz).
-  int32_t dopplerMin;
-
-  /// @brief Maximum Doppler (Hz).
-  int32_t dopplerMax;
-
-  /// @brief Sampling frequency (Hz).
-  uint32_t fs;
-
-  /// @brief Number of samples.
-  uint32_t n;
-
-  /// @brief Center of Doppler bins (Hz).
-  double dopplerMiddle;
-
-  /// @brief Number of delay bins.
-  uint16_t nDelayBins;
-
-  /// @brief Number of Doppler bins.
-  uint16_t nDopplerBins;
-
-  /// @brief Number of correlation samples per pulse.
-  uint16_t nCorr;
-
-  /// @brief True CPI time (s).
-  double cpi;
-
-  /// @brief FFTW plans for ambiguity processing.
-  /// @{
-  fftw_plan fftXi, fftYi, fftZi, fftDoppler;
-  /// @}
-
-  /// @brief FFTW storage for ambiguity processing.
-  /// @{
-  std::complex<double> *dataXi, *dataYi, *dataZi, *dataCorr, *dataDoppler;
-  /// @}
-
-  /// @brief Number of samples to perform FFT per pulse.
-  uint32_t nfft;
-
-  /// @brief Vector storage for ambiguity processing
-  /// @{
-  std::vector<std::complex<double>> corr, delayProfile;
-  /// @}
-
-  /// @brief Pointer to map to store result.
-  Map<std::complex<double>> *map;
 
 public:
+
+  using Complex = std::complex<double>;
+
   /// @brief Constructor.
   /// @param delayMin Minimum delay (bins).
   /// @param delayMax Maximum delay (bins).
@@ -90,7 +39,80 @@ public:
   /// @param x Reference samples.
   /// @param y Surveillance samples.
   /// @return Ambiguity map data of IQ samples.
-  Map<std::complex<double>> *process(IqData *x, IqData *y);
-};
+  Map<Complex> *process(IqData *x, IqData *y);
 
-#endif
+  double doppler_middle() const { return dopplerMiddle_; }
+
+  uint16_t delay_bin_count() const { return nDelayBins_; }
+
+  uint16_t doppler_bin_count() const { return nDopplerBins_; }
+
+  uint16_t corr_samples_per_pulse() const { return nCorr_; }
+
+  double cpi_length_seconds() const { return cpi_; }
+
+  uint32_t fft_bin_count() const { return nfft_; }
+
+  double doppler_res_;
+private:
+  /// @brief Minimum delay (bins).
+  int32_t delayMin_;
+
+  /// @brief Maximum delay (bins).
+  int32_t delayMax_;
+
+  /// @brief Minimum Doppler (Hz).
+  int32_t dopplerMin_;
+
+  /// @brief Maximum Doppler (Hz).
+  int32_t dopplerMax_;
+
+  /// @brief Sampling frequency (Hz).
+  uint32_t fs_;
+
+  /// @brief Number of samples.
+  uint32_t nSamples_;
+
+  /// @brief Center of Doppler bins (Hz).
+  double dopplerMiddle_;
+
+  /// @brief Number of delay bins.
+  uint16_t nDelayBins_;
+
+  /// @brief Number of Doppler bins.
+  uint16_t nDopplerBins_;
+
+  /// @brief Number of correlation samples per pulse.
+  uint16_t nCorr_;
+
+  /// @brief True CPI time (s).
+  double cpi_;
+
+  /// @brief FFTW plans for ambiguity processing.
+  fftw_plan fftXi_;
+  fftw_plan fftYi_;
+  fftw_plan fftZi_;
+  fftw_plan fftDoppler_;
+
+  /// @brief FFTW storage for ambiguity processing.
+  /// @{
+  std::vector<Complex> dataXi_;
+  std::vector<Complex> dataYi_;
+  std::vector<Complex> dataZi_;
+  std::vector<Complex> dataCorr_;
+  std::vector<Complex> dataDoppler_;
+  /// @}
+
+  /// @brief Number of samples to perform FFT per pulse.
+  uint32_t nfft_;
+
+  /// @brief Vector storage for ambiguity processing
+  /// @{
+  std::vector<Complex> corr_;
+  std::vector<Complex> delayProfile_;
+  /// @}
+
+  /// @brief Map to store result.
+  std::unique_ptr<Map<Complex>> map_;
+
+};

src/process/ambiguity/test_ambiguity.cpp

@@ -0,0 +1,123 @@
+#define CATCH_CONFIG_MAIN
+#include "catch_amalgamated.hpp"
+
+#include "Ambiguity.h"
+#include <random>
+
+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());
+    std::uniform_real_distribution<> dist(-100.0, 100.0);
+
+    for (uint32_t i = 0; i < iq_data.get_n(); ++i) {
+        iq_data.push_back({dist(gen), dist(gen)});
+    }
+}
+
+void read_file(IqData& buffer1, IqData& buffer2, const std::string& file)
+{
+  short i1, q1, i2, q2;
+  auto file_replay = fopen(file.c_str(), "rb");
+  if (!file_replay) {
+    return;
+  }
+
+  auto read_short = [](short& v, FILE* fid) {
+    auto rv{fread(&v, 1, sizeof(short), fid)};
+    return rv == sizeof(short);
+  };
+
+  while (!feof(file_replay))
+  {
+    if (!read_short(i1, file_replay)) break;
+    if (!read_short(q1, file_replay)) break;
+    if (!read_short(i2, file_replay)) break;
+    if (!read_short(q2, file_replay)) break;
+
+    buffer1.push_back({(double)i1, (double)q1});
+    buffer2.push_back({(double)i2, (double)q2});
+
+    // Only read for the buffer length - this class is very poorly designed.
+    if (buffer1.get_length() == buffer1.get_n()) {
+        break;
+    }
+  }
+
+  fclose(file_replay);
+}
+
+// Make sure the constructor is calculating the parameters correctly.
+TEST_CASE("Constructor", "[constructor]")
+{
+    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);
+
+    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() == 6643);
+}
+
+// Make sure process produces an output
+TEST_CASE("Process_Simple", "[process]")
+{
+    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);
+
+    IqData x{n_samples};
+    IqData y{n_samples};
+
+    random_iq(x);
+    random_iq(y);
+    auto map{ambiguity.process(&x, &y)};
+    map->set_metrics();
+    CHECK(map->maxPower > 0.0);
+    CHECK(map->noisePower > 0.0);
+}
+
+// Sanity check that we're getting numbers close to the baseline ambiguity processing function.
+TEST_CASE("Process_File", "[process]")
+{
+    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);
+    IqData x{n_samples};
+    IqData y{n_samples};
+
+    read_file(x, y, "20231214-230611.rspduo");
+    REQUIRE(x.get_length() == x.get_n());
+
+    auto map{ambiguity.process(&x ,&y)};
+    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));
+}