calculate channel based on mvue

include/iq_stream_decoder.hpp

@@ -13,6 +13,7 @@
 #include "fixed_queue.hpp"
 #include "l2/lower_mac.hpp"
 #include "streaming_ordered_output_thread_pool_executor.hpp"
+#include <armadillo>
 #include <complex>
 #include <memory>
 
@@ -23,6 +24,8 @@
  */
 class IQStreamDecoder {
   public:
+    using QueueT = FixedQueue<std::complex<float>, 300>;
+
     IQStreamDecoder(
         const std::shared_ptr<StreamingOrderedOutputThreadPoolExecutor<LowerMac::return_type>>& lower_mac_worker_queue,
         const std::shared_ptr<LowerMac>& lower_mac, const std::shared_ptr<BitStreamDecoder>& bit_stream_decoder,
@@ -31,19 +34,18 @@ class IQStreamDecoder {
 
     void process_complex(std::complex<float> symbol) noexcept;
 
-  private:
-    using QueueT = FixedQueue<std::complex<float>, 300>;
+    template <std::size_t ChannelSize>
+    static auto solve_channel(const std::vector<std::complex<float>>& pilots, const QueueT& signal_queue,
+                              std::size_t signal_offset) -> arma::cx_fvec;
 
+  private:
     static std::complex<float> hard_decision(std::complex<float> const& symbol);
 
     template <class iterator_type> static void symbols_to_bitstream(iterator_type it, uint8_t* bits, std::size_t len);
 
     static void abs_convolve_same_length(const QueueT& queueA, std::size_t offsetA, const std::complex<float>* itb,
                                          std::size_t len, float* res);
 
-    static auto solve_channel(const std::vector<std::complex<float>>& pilots,
-                          const QueueT& signal_queue, const std::size_t signal_offset);
-
     std::vector<std::complex<float>> channel_estimation(std::vector<std::complex<float>> const& stream,
                                                         std::vector<std::complex<float>> const& pilots);
 

src/experiments/channel_estimation.cpp

@@ -6,10 +6,39 @@
  *   Marenz Schmidl
  */
 
+#include "iq_stream_decoder.hpp"
 #include <iostream>
 
-auto main(int argc, char** argv) -> int {
-    auto decoder = IQStreamDecoder(nullptr, nullptr, nullptr, /*is_uplink=*/true);
+auto main(int /*argc*/, char** /*argv*/) -> int {
+    const std::vector<std::complex<float>> training_seq_x = {
+        {1.0, -1.0}, {-1.0, 1.0}, {-1.0, -1.0}, {-1.0, 1.0}, {1.0, 1.0}, {1.0, 1.0}, {-1.0, -1.0}, {1.0, -1.0},
+        {1.0, -1.0}, {-1.0, 1.0}, {-1.0, -1.0}, {-1.0, 1.0}, {1.0, 1.0}, {1.0, 1.0}, {-1.0, -1.0}};
+
+    // Fill the queue with the training sequence
+    IQStreamDecoder::QueueT queue;
+    for (auto i = 0; i < training_seq_x.size(); i++) {
+        std::complex<float> v = {0, 0};
+        // if (i > 0) {
+        //     std::complex<float> k = {-0.5, -0.5};
+        //     v += training_seq_x[i - 1] * k;
+        // }
+        // if (i < training_seq_x.size() - 1) {
+        //     std::complex<float> k = {-0.5, -0.5};
+        //     v += training_seq_x[i + 1] * k;
+        // }
+        std::complex<float> k = {0.5, 0.2};
+        v += training_seq_x[i] * k;
+
+        queue.push(v);
+    }
+    // for (const auto& elem : training_seq_x) {
+    //     queue.push(elem);
+    // }
+
+    auto result =
+        IQStreamDecoder::solve_channel<3>(training_seq_x, queue, /*signal_offset=*/300 - training_seq_x.size());
+
+    std::cout << result << '\n';
 
     return EXIT_SUCCESS;
 }

src/iq_stream_decoder.cpp

@@ -95,17 +95,50 @@ std::vector<std::complex<float>> IQStreamDecoder::channel_estimation(std::vector
     return stream;
 }
 
- auto IQStreamDecoder::solve_channel(const std::vector<std::complex<float>>& pilots,
-                          const QueueT& signal_queue, const std::size_t signal_offset)
-    -> arma::cx_fvec {
-    auto arma_pilots = arma::cx_fvec(pilots);
-    auto arma_signal = arma::cx_fvec(pilots.size());
-    for (auto i = 0; i < arma_signal.size(); i++) {
-        arma_signal[i] = signal_queue[signal_offset + i];
+template <std::size_t ChannelSize>
+auto IQStreamDecoder::solve_channel(const std::vector<std::complex<float>>& pilots, const QueueT& signal_queue,
+                                    const std::size_t signal_offset) -> arma::cx_fvec {
+    // Calculate the minimum variance unbiased estimator
+    auto h = arma::cx_fmat(/*n_rows=*/pilots.size(), /*n_cols=*/ChannelSize, arma::fill::zeros);
+
+    for (auto row = 0; row < h.n_rows; row++) {
+        auto index = row;
+        for (auto col = 0; col < h.n_cols; col++) {
+            if (index >= 0) {
+                h.row(row).col(col) = signal_queue[signal_offset + index];
+            }
+            index--;
+        }
     }
-    auto arma_conj_pilots = arma::conj(arma_pilots);
-    arma::cx_fvec h_vec = arma::solve(arma_conj_pilots.t() * arma_pilots, arma_conj_pilots.t() * arma::conj(arma_signal).t());
-    return h_vec;
+
+    // std::cout << h << '\n';
+
+    auto signal = arma::cx_fvec(pilots.size());
+    for (auto i = 0; i < signal.size(); i++) {
+        signal[i] = signal_queue[signal_offset + i];
+    }
+
+    // std::cout << signal << '\n';
+
+    auto h_hermitian = arma::cx_fmat(/*n_rows=*/ChannelSize, /*n_cols=*/pilots.size(), arma::fill::zeros);
+
+    h_hermitian = h.t();
+    h_hermitian = arma::conj(h_hermitian);
+    // std::cout << h_hermitian << '\n';
+
+    auto h_hermitian_h = h_hermitian * h;
+    // std::cout << h_hermitian_h << '\n';
+
+    // auto h_hermitian_h_inv = arma::inv(h_hermitian_h);
+    // std::cout << h_hermitian_h_inv << '\n';
+
+    auto h_hermitian_signal = h_hermitian * signal;
+    // std::cout << h_hermitian_signal << '\n';
+
+    auto c_mvue = arma::solve(h_hermitian_h, h_hermitian_signal);
+    // std::cout << c_mvue << '\n';
+
+    return c_mvue;
 }
 
 void IQStreamDecoder::process_complex(std::complex<float> symbol) noexcept {
@@ -137,7 +170,7 @@ void IQStreamDecoder::process_complex(std::complex<float> symbol) noexcept {
         if (detectedX >= SEQUENCE_DETECTION_THRESHOLD) {
             // std::cout << "Potential CUB found" << std::endl;
 
-            auto channel = solve_channel(training_seq_x_, symbol_buffer_hard_decision_, 44);
+            auto channel = solve_channel<3>(training_seq_x_, symbol_buffer_hard_decision_, 44);
             std::cout << channel << std::endl;
 
             auto len = 103;