Added complex noise function

This was an attempt at getting the divide by zeros to stop, but if this didn't do it then it might be that the issue isn't a divide by zero.  Will need to create a MATLAB test to find out what's actually happening

gnuradio/gr-droneid/CMakeLists.txt

@@ -25,6 +25,8 @@ cmake_minimum_required(VERSION 3.8)
 project(gr-droneid CXX C)
 enable_testing()
 
+set(MATLAB_PATH "" CACHE STRING "Path to MATLAB installation")
+
 # Install to PyBOMBS target prefix if defined
 if(DEFINED ENV{PYBOMBS_PREFIX})
     set(CMAKE_INSTALL_PREFIX $ENV{PYBOMBS_PREFIX})

gnuradio/gr-droneid/grc/CMakeLists.txt

@@ -25,5 +25,7 @@ install(FILES
     droneid_misc_utils.block.yml
     droneid_lte_decode.block.yml
     droneid_decode.block.yml
-    droneid_normalized_xcorr_estimate.block.yml DESTINATION share/gnuradio/grc/blocks
+    droneid_normalized_xcorr_estimate.block.yml
+    droneid_variance.block.yml
+    droneid_dot_prod.block.yml DESTINATION share/gnuradio/grc/blocks
 )

gnuradio/gr-droneid/grc/droneid_dot_prod.block.yml

@@ -0,0 +1,42 @@
+id: droneid_dot_prod
+label: dot_prod
+category: '[droneid]'
+
+templates:
+  imports: import droneid
+  make: droneid.dot_prod(${taps})
+
+#  Make one 'parameters' list entry for every parameter you want settable from the GUI.
+#     Keys include:
+#     * id (makes the value accessible as \$keyname, e.g. in the make entry)
+#     * label (label shown in the GUI)
+#     * dtype (e.g. int, float, complex, byte, short, xxx_vector, ...)
+parameters:
+- id: taps
+  label: Filter Taps
+  dtype: complex_vector
+
+#  Make one 'inputs' list entry per input and one 'outputs' list entry per output.
+#  Keys include:
+#      * label (an identifier for the GUI)
+#      * domain (optional - stream or message. Default is stream)
+#      * dtype (e.g. int, float, complex, byte, short, xxx_vector, ...)
+#      * vlen (optional - data stream vector length. Default is 1)
+#      * optional (optional - set to 1 for optional inputs. Default is 0)
+inputs:
+- label: in
+  domain: stream
+  dtype: complex
+  vlen: 1
+  optional: 0
+
+outputs:
+- label: out
+  domain: stream
+  dtype: complex
+  vlen: 1
+  optional: 0
+
+#  'file_format' specifies the version of the GRC yml format used in the file
+#  and should usually not be changed.
+file_format: 1

gnuradio/gr-droneid/grc/droneid_variance.block.yml

@@ -0,0 +1,42 @@
+id: droneid_variance
+label: variance
+category: '[droneid]'
+
+templates:
+  imports: import droneid
+  make: droneid.variance(${window_size})
+
+#  Make one 'parameters' list entry for every parameter you want settable from the GUI.
+#     Keys include:
+#     * id (makes the value accessible as \$keyname, e.g. in the make entry)
+#     * label (label shown in the GUI)
+#     * dtype (e.g. int, float, complex, byte, short, xxx_vector, ...)
+parameters:
+- id: window_size
+  label: Window Size
+  dtype: int
+
+#  Make one 'inputs' list entry per input and one 'outputs' list entry per output.
+#  Keys include:
+#      * label (an identifier for the GUI)
+#      * domain (optional - stream or message. Default is stream)
+#      * dtype (e.g. int, float, complex, byte, short, xxx_vector, ...)
+#      * vlen (optional - data stream vector length. Default is 1)
+#      * optional (optional - set to 1 for optional inputs. Default is 0)
+inputs:
+- label: in
+  domain: stream
+  dtype: complex
+  vlen: 1
+  optional: 0
+
+outputs:
+- label: out
+  domain: stream
+  dtype: float
+  vlen: 1
+  optional: 0
+
+#  'file_format' specifies the version of the GRC yml format used in the file
+#  and should usually not be changed.
+file_format: 1

gnuradio/gr-droneid/include/droneid/CMakeLists.txt

@@ -31,5 +31,7 @@ install(FILES
     lte_decode.h
     decode.h
     normalized_xcorr.h
-    normalized_xcorr_estimate.h DESTINATION include/droneid
+    normalized_xcorr_estimate.h
+    variance.h
+    dot_prod.h DESTINATION include/droneid
 )

gnuradio/gr-droneid/include/droneid/dot_prod.h

@@ -0,0 +1,54 @@
+/* -*- c++ -*- */
+/*
+ * Copyright 2022 gr-droneid author.
+ *
+ * This is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3, or (at your option)
+ * any later version.
+ *
+ * This software is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this software; see the file COPYING.  If not, write to
+ * the Free Software Foundation, Inc., 51 Franklin Street,
+ * Boston, MA 02110-1301, USA.
+ */
+
+#ifndef INCLUDED_DRONEID_DOT_PROD_H
+#define INCLUDED_DRONEID_DOT_PROD_H
+
+#include <droneid/api.h>
+#include <gnuradio/block.h>
+
+namespace gr {
+    namespace droneid {
+
+        /*!
+         * \brief <+description of block+>
+         * \ingroup droneid
+         *
+         */
+        class DRONEID_API dot_prod : virtual public gr::block {
+        public:
+            typedef boost::shared_ptr<dot_prod> sptr;
+
+            /*!
+             * \brief Return a shared_ptr to a new instance of droneid::dot_prod.
+             *
+             * To avoid accidental use of raw pointers, droneid::dot_prod's
+             * constructor is in a private implementation
+             * class. droneid::dot_prod::make is the public interface for
+             * creating new instances.
+             */
+            static sptr make(const std::vector<gr_complex> & /*taps*/);
+        };
+
+    } // namespace droneid
+} // namespace gr
+
+#endif /* INCLUDED_DRONEID_DOT_PROD_H */
+

gnuradio/gr-droneid/include/droneid/misc_utils.h

@@ -43,6 +43,8 @@ namespace gr {
             static uint32_t get_short_cp_len(double sample_rate);
             static uint32_t get_fft_size(double sample_rate);
 
+            static std::vector<std::complex<float>> create_gaussian_noise(uint32_t sample_count);
+
             static std::vector<std::complex<float>> create_zc_sequence(double sample_rate, uint32_t root);
             static std::vector<std::complex<float>> conj(const std::vector<std::complex<float>> & input);
 

gnuradio/gr-droneid/include/droneid/variance.h

@@ -0,0 +1,55 @@
+/* -*- c++ -*- */
+/*
+ * Copyright 2022 gr-droneid author.
+ *
+ * This is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3, or (at your option)
+ * any later version.
+ *
+ * This software is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this software; see the file COPYING.  If not, write to
+ * the Free Software Foundation, Inc., 51 Franklin Street,
+ * Boston, MA 02110-1301, USA.
+ */
+
+#ifndef INCLUDED_DRONEID_VARIANCE_H
+#define INCLUDED_DRONEID_VARIANCE_H
+
+#include <droneid/api.h>
+#include <gnuradio/sync_block.h>
+
+namespace gr {
+  namespace droneid {
+
+    /*!
+     * \brief <+description of block+>
+     * \ingroup droneid
+     *
+     */
+    class DRONEID_API variance : virtual public gr::block
+    {
+     public:
+      typedef boost::shared_ptr<variance> sptr;
+
+      /*!
+       * \brief Return a shared_ptr to a new instance of droneid::variance.
+       *
+       * To avoid accidental use of raw pointers, droneid::variance's
+       * constructor is in a private implementation
+       * class. droneid::variance::make is the public interface for
+       * creating new instances.
+       */
+      static sptr make(uint32_t window_size);
+    };
+
+  } // namespace droneid
+} // namespace gr
+
+#endif /* INCLUDED_DRONEID_VARIANCE_H */
+

gnuradio/gr-droneid/lib/CMakeLists.txt

@@ -33,6 +33,8 @@ list(APPEND droneid_sources
     decode_impl.cc
     normalized_xcorr.cc
     normalized_xcorr_estimate_impl.cc
+    variance_impl.cc
+    dot_prod_impl.cc
 )
 
 set(droneid_sources "${droneid_sources}" PARENT_SCOPE)
@@ -80,9 +82,32 @@ include(GrTest)
 #include_directories()
 # List all files that contain Boost.UTF unit tests here
 list(APPEND test_droneid_sources
-        qa_normalized_xcorr.cc
-        qa_normalized_xcorr_estimate.cc
+    qa_normalized_xcorr_estimate.cc
 )
+
+if (MATLAB_PATH)
+    message(STATUS "Using MATLAB path '${MATLAB_PATH}'")
+
+    if (NOT IS_DIRECTORY ${MATLAB_PATH})
+        message(FATAL_ERROR "MATLAB path '${MATLAB_PATH}' does not exist or is not a directory.  Please fix or remove")
+    endif()
+
+    # Add paths to MATLAB libraries and headers
+    target_link_directories(gnuradio-droneid PUBLIC ${MATLAB_PATH}/extern/bin/glnxa64)
+    target_include_directories(gnuradio-droneid PUBLIC PUBLIC ${MATLAB_PATH}/extern/include)
+
+    # Add MATLAB library dependencies
+    target_link_libraries(gnuradio-droneid MatlabEngine MatlabDataArray)
+
+    list(APPEND test_droneid_sources
+        qa_variance.cc
+        qa_normalized_xcorr.cc
+        qa_dot_prod.cc
+    )
+else()
+    message(WARNING "No MATLAB path specified, so some tests will be skipped")
+endif()
+
 # Anything we need to link to for the unit tests go here
 list(APPEND GR_TEST_TARGET_DEPS gnuradio-droneid)
 

gnuradio/gr-droneid/lib/dot_prod_impl.cc

@@ -0,0 +1,179 @@
+/* -*- c++ -*- */
+/*
+ * Copyright 2022 gr-droneid author.
+ *
+ * This is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3, or (at your option)
+ * any later version.
+ *
+ * This software is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this software; see the file COPYING.  If not, write to
+ * the Free Software Foundation, Inc., 51 Franklin Street,
+ * Boston, MA 02110-1301, USA.
+ */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <gnuradio/io_signature.h>
+#include "dot_prod_impl.h"
+#include <volk/volk.h>
+#include <numeric>
+#include <droneid/misc_utils.h>
+
+namespace gr {
+    namespace droneid {
+
+        dot_prod::sptr
+        dot_prod::make(const std::vector<gr_complex> &taps) {
+            return gnuradio::get_initial_sptr
+                    (new dot_prod_impl(taps));
+        }
+
+
+        /*
+         * The private constructor
+         */
+        dot_prod_impl::dot_prod_impl(const std::vector<gr_complex> &taps)
+                : gr::block("dot_prod",
+                            gr::io_signature::make(1, 1, sizeof(gr_complex)),
+                            gr::io_signature::make(1, 1, sizeof(gr_complex))),
+                  taps_(taps), window_size_(taps.size()) {
+
+            // Remove the mean from the taps, conjugate the taps, and calculate the variance ahead of time
+            const auto mean =
+                    std::accumulate(taps_.begin(), taps_.end(), gr_complex{0, 0}) / static_cast<float>(taps_.size());
+
+            for (auto & tap : taps_) {
+                tap = std::conj(tap) - mean;
+            }
+
+            taps_var_ = misc_utils::var_no_mean(taps_);
+
+            // Create some constants to enable the use of multiplies instead of divides later
+            window_size_recip_ = 1.0f / static_cast<float>(window_size_);
+            window_size_recip_complex_ = gr_complex{window_size_recip_, 0};
+        }
+
+        /*
+         * Our virtual destructor.
+         */
+        dot_prod_impl::~dot_prod_impl() {
+        }
+
+        int
+        dot_prod_impl::general_work(int noutput_items,
+                                    gr_vector_int &ninput_items,
+                                    gr_vector_const_void_star &input_items,
+                                    gr_vector_void_star &output_items) {
+            const auto *in = (const gr_complex *) input_items[0];
+            auto *out = (gr_complex *) output_items[0];
+
+            consume_each(noutput_items);
+
+            // This is how the remaining samples are buffered between calls.  It's important to realize that this algo
+            // needs <window_size> samples to be able to produce one output value.  This means that there will always
+            // be unused samples at the end of each function call that need to be held onto until the next call.  The
+            // hope was that set_history() took care of this, but it does not.  So, the remaining samples from the last
+            // call are stored in <buffer_>.  The <in> buffer can't hold more samples (it's not known how many samples
+            // wide the buffer is) so in order to use the old samples without jumping through very slow hoops, the new
+            // samples are appended to the old samples.
+            buffer_.insert(buffer_.end(), in, in + noutput_items);
+
+            // Exit early if there aren't enough samples to process.
+            if (buffer_.size() < window_size_) {
+                return 0;
+            }
+
+            // Figure out how many windows worth of data can be processed.  It's possible that this specific call
+            // doesn't have enough storage in its output buffer to hold all the samples that could be processed.  For
+            // this reason the min of the available output buffer space and number of windows that could be processed
+            // must be used.
+            const auto num_steps = std::min(static_cast<uint64_t>(noutput_items), buffer_.size() - window_size_);
+
+            // Resize the buffers as needed
+            if (sums_.size() < num_steps) {
+                sums_.resize(num_steps);
+                abs_squared_.resize(num_steps + window_size_);
+                vars_.resize(num_steps);
+            }
+
+            // TODO(24June2022): There are <window_size-1> extra operations happening on each call.  This comes from the
+            //                   fact that some of these computations are being done on samples that are going to be
+            //                   used again on the next function call.  Would be a good idea to buffer the abs squared
+            //                   and maybe the running variance average.
+
+            // What is happening below is roughly the following:
+            //
+            // for idx = 1:length(buffer_) - window_size_
+            //   window = buffer_(idx:idx + window_size_ - 1);
+            //   variance = sum(abs(window).^2) / window_size_;
+            //   dot_prod = sum(window .* taps_) / window_size_;
+            //   out(idx) = dot_prod / sqrt(variance * taps_var_);
+            // end
+            //
+            // But the variance is calculated as a running sum.  The first variance has to be calculated the hard way,
+            // and then every iteration of the loop will subtract off the left-most element of the window that just
+            // dropped off, and adds on the new right-most element in the window.
+            //
+            // Doing this calculation of the first element outside the loop prevents needing a conditional in the
+            // critical section
+
+            // Calculate the first variance the hard way
+            volk_32fc_magnitude_squared_32f(&abs_squared_[0], &buffer_[0], num_steps + window_size_);
+            auto running_var = std::accumulate(abs_squared_.begin(), abs_squared_.begin() + window_size_, 0.f);
+            vars_[0] = running_var;
+
+            // Calculate the first dot product
+            volk_32fc_x2_dot_prod_32fc(&out[0], &buffer_[0], &taps_[0], window_size_);
+
+            // Calculate the running abs value sum and dot product for the remaining samples
+            for (uint32_t idx = 1; idx < num_steps; idx++) {
+                // sum(abs(window).^2)
+                running_var = running_var - abs_squared_[idx - 1] + abs_squared_[idx + window_size_];
+                vars_[idx] = running_var;
+
+                // Compute tue dot product of the current window and the filter taps
+                // sum(window .* taps_)
+                volk_32fc_x2_dot_prod_32fc(&out[idx], &buffer_[idx], &taps_[0], window_size_);
+            }
+
+            // Scale the dot product down
+            volk_32fc_s32fc_multiply_32fc(&out[0], &out[0], window_size_recip_complex_, num_steps);
+
+            // Scale the variance sums down
+            volk_32f_s32f_multiply_32f(&vars_[0], &vars_[0], window_size_recip_, num_steps);
+
+            // Multiply each variance by the tap variances then take the reciprocal
+            volk_32f_s32f_multiply_32f(&vars_[0], &vars_[0], taps_var_, num_steps);
+
+            // Take the square root of the product of the two variances
+            volk_32f_sqrt_32f(&vars_[0], &vars_[0], num_steps);
+
+            // There's no VOLK function for the reciprocal operation.  This is being done so that a multiply can be
+            // used next to divide the dot product results by the sqrt calculated above
+            for (auto & var : vars_) {
+                var = 1.0f / var;
+            }
+
+            // Divide by the square root above
+            volk_32fc_32f_multiply_32fc(&out[0], &out[0], &vars_[0], num_steps);
+
+            // Remove all the samples that have been processed from the buffer.  Leaving just the last <window_size_-1>
+            // samples for the next call
+            buffer_.erase(buffer_.begin(), buffer_.begin() + num_steps);
+
+            // Tell runtime system how many output items we produced.
+            return num_steps;
+        }
+
+    } /* namespace droneid */
+} /* namespace gr */
+

gnuradio/gr-droneid/lib/dot_prod_impl.h

@@ -0,0 +1,61 @@
+/* -*- c++ -*- */
+/*
+ * Copyright 2022 gr-droneid author.
+ *
+ * This is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3, or (at your option)
+ * any later version.
+ *
+ * This software is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this software; see the file COPYING.  If not, write to
+ * the Free Software Foundation, Inc., 51 Franklin Street,
+ * Boston, MA 02110-1301, USA.
+ */
+
+#ifndef INCLUDED_DRONEID_DOT_PROD_IMPL_H
+#define INCLUDED_DRONEID_DOT_PROD_IMPL_H
+
+#include <droneid/dot_prod.h>
+#include <queue>
+
+namespace gr {
+    namespace droneid {
+
+        class dot_prod_impl : public dot_prod {
+        private:
+            const uint32_t window_size_;
+            float taps_var_;
+            float window_size_recip_;
+            gr_complex window_size_recip_complex_;
+            std::vector<gr_complex> taps_;
+            std::vector<gr_complex> sums_;
+            std::vector<float> vars_;
+            std::vector<float> abs_squared_;
+
+
+            std::vector<gr_complex> buffer_;
+            // Nothing to declare in this block.
+
+        public:
+            dot_prod_impl(const std::vector<gr_complex> & taps);
+
+            ~dot_prod_impl();
+
+            int general_work(int noutput_items,
+                             gr_vector_int &ninput_items,
+                             gr_vector_const_void_star &input_items,
+                             gr_vector_void_star &output_items);
+
+        };
+
+    } // namespace droneid
+} // namespace gr
+
+#endif /* INCLUDED_DRONEID_DOT_PROD_IMPL_H */
+

gnuradio/gr-droneid/lib/misc_utils.cc

@@ -28,6 +28,7 @@
 #include <droneid/misc_utils.h>
 #include <gnuradio/fft/fft.h>
 #include <gnuradio/fft/fft_shift.h>
+#include <gnuradio/random.h>
 #include <iostream>
 
 #include <volk/volk.h>
@@ -398,6 +399,16 @@ namespace gr {
             return abs_squared(samples);
         }
 
+        std::vector<std::complex<float>> misc_utils::create_gaussian_noise(const uint32_t sample_count) {
+            std::vector<std::complex<float>> output(sample_count);
+            gr::random rand;
+            for (auto & sample : output) {
+                sample = rand.rayleigh_complex();
+            }
+
+            return output;
+        }
+
 
     } /* namespace droneid */
 } /* namespace gr */

gnuradio/gr-droneid/lib/qa_dot_prod.cc

@@ -0,0 +1,64 @@
+/* -*- c++ -*- */
+/*
+ * Copyright 2022 gr-droneid author.
+ *
+ * This is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3, or (at your option)
+ * any later version.
+ *
+ * This software is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this software; see the file COPYING.  If not, write to
+ * the Free Software Foundation, Inc., 51 Franklin Street,
+ * Boston, MA 02110-1301, USA.
+ */
+
+#include <iostream>
+
+#include <gnuradio/attributes.h>
+#include "qa_dot_prod.h"
+
+#include <droneid/dot_prod.h>
+#include <droneid/misc_utils.h>
+
+#include <boost/test/unit_test.hpp>
+
+#include <gnuradio/blocks/vector_source.h>
+#include <gnuradio/blocks/vector_sink.h>
+#include <gnuradio/blocks/file_source.h>
+#include <gnuradio/top_block.h>
+
+namespace gr {
+    namespace droneid {
+
+
+
+        BOOST_AUTO_TEST_CASE(moocow_test) {
+            auto tb = gr::make_top_block("top");
+
+            const auto noise = misc_utils::create_gaussian_noise(12000);
+            const auto taps_offset = 4562;
+            const auto taps_size = 1024;
+            const auto taps = std::vector<gr_complex>(noise.begin() + taps_offset, noise.begin() + taps_offset + taps_size);
+
+            auto source = gr::blocks::vector_source<gr_complex>::make(noise);
+            auto sink = gr::blocks::vector_sink<gr_complex>::make();
+
+            auto uut = droneid::dot_prod::make(taps);
+
+            tb->connect(source, 0, uut, 0);
+            tb->connect(uut, 0, sink, 0);
+
+            tb->run();
+
+            std::cout << "Sent in " << noise.size() << " samples, got back " << sink->data().size() << " samples\n";
+        }
+
+    } /* namespace droneid */
+} /* namespace gr */
+

gnuradio/gr-droneid/lib/qa_dot_prod.h

@@ -0,0 +1,45 @@
+/* -*- c++ -*- */
+/*
+ * Copyright 2022 gr-droneid author.
+ *
+ * This is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3, or (at your option)
+ * any later version.
+ *
+ * This software is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this software; see the file COPYING.  If not, write to
+ * the Free Software Foundation, Inc., 51 Franklin Street,
+ * Boston, MA 02110-1301, USA.
+ */
+
+#ifndef _QA_DOT_PROD_H_
+#define _QA_DOT_PROD_H_
+
+#include <cppunit/extensions/HelperMacros.h>
+#include <cppunit/TestCase.h>
+
+namespace gr {
+  namespace droneid {
+
+    class qa_dot_prod : public CppUnit::TestCase
+    {
+    public:
+      CPPUNIT_TEST_SUITE(qa_dot_prod);
+      CPPUNIT_TEST(t1);
+      CPPUNIT_TEST_SUITE_END();
+
+    private:
+      void t1();
+    };
+
+  } /* namespace droneid */
+} /* namespace gr */
+
+#endif /* _QA_DOT_PROD_H_ */
+

gnuradio/gr-droneid/lib/qa_variance.cc

@@ -0,0 +1,142 @@
+/* -*- c++ -*- */
+/*
+ * Copyright 2022 gr-droneid author.
+ *
+ * This is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3, or (at your option)
+ * any later version.
+ *
+ * This software is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this software; see the file COPYING.  If not, write to
+ * the Free Software Foundation, Inc., 51 Franklin Street,
+ * Boston, MA 02110-1301, USA.
+ */
+
+#include <gnuradio/attributes.h>
+#include <cppunit/TestAssert.h>
+#include "qa_variance.h"
+#include <droneid/variance.h>
+#include <boost/test/unit_test.hpp>
+
+#include <gnuradio/blocks/vector_source.h>
+#include <gnuradio/blocks/vector_sink.h>
+#include <gnuradio/top_block.h>
+#include <gnuradio/random.h>
+
+#include <MatlabEngine.hpp>
+#include <MatlabDataArray.hpp>
+
+namespace gr {
+    namespace droneid {
+
+        std::vector<gr_complex> create_noise(const uint32_t sample_count) {
+            std::vector<gr_complex> samples(sample_count);
+            gr::random rng;
+            for (auto & sample : samples) {
+                sample = rng.rayleigh_complex();
+            }
+
+            return samples;
+        }
+
+        std::vector<float> calculate_true_variance(const std::vector<gr_complex> & samples, const uint32_t window_size) {
+            using namespace matlab::engine;
+            using namespace matlab::data;
+
+            ArrayFactory factory;
+
+            auto matlab = startMATLAB();
+
+            auto samples_matlab = factory.createArray<gr_complex>(
+                    ArrayDimensions({samples.size(), 1}), &samples[0], &samples[samples.size() - 1]);
+
+            matlab->setVariable("samples", samples_matlab);
+            matlab->setVariable("window_size", factory.createScalar(window_size));
+
+            const auto program =
+                    u"scores = zeros(length(samples) - window_size, 1);"
+                    u"for idx = 1:length(scores)"
+                    u"  window = samples(idx:idx + window_size - 1);"
+                    u"  scores(idx) = var(window);"
+                    u"end";
+
+            matlab->eval(program);
+
+            const auto scores = matlab->getVariable("scores");
+
+            std::vector<float> output(scores.getNumberOfElements());
+            for (uint32_t idx = 0; idx < output.size(); idx++) {
+                output[idx] = static_cast<float>(scores[idx]);
+            }
+
+            return output;
+        }
+
+        void plot_results(const std::vector<float> & expected_values, const std::vector<float> & uut_output,
+                          const uint32_t window_size) {
+            auto matlab = matlab::engine::startMATLAB();
+            matlab::data::ArrayFactory factory;
+
+            const auto matlab_values = factory.createArray<float>(
+                    matlab::data::ArrayDimensions({expected_values.size(), 1}),
+                    &expected_values[0], &expected_values[expected_values.size() - 1]);
+
+            const auto uut_values = factory.createArray<float>(
+                    matlab::data::ArrayDimensions({uut_output.size(), 1}),
+                    &uut_output[0], &uut_output[uut_output.size() - 1]);
+
+            matlab->setVariable("matlab", matlab_values);
+            matlab->setVariable("uut", uut_values);
+            matlab->setVariable("window_size", factory.createScalar(window_size));
+
+            matlab->eval(
+                    u"figure(1);"
+                    u"subplot(3, 1, 1); plot(matlab); title('MATLAB');"
+                    u"subplot(3, 1, 2); plot(uut); title('UUT');"
+                    u"subplot(3, 1, 3); plot(uut - matlab); title('Delta');"
+                    u"pause");
+        }
+
+        BOOST_AUTO_TEST_CASE(test_variance) {
+            const float max_delta = 60e-3;
+            auto tb = gr::make_top_block("top");
+
+            const auto sample_count = 12000;
+            const auto window_size = 1024;
+
+            const auto noise = create_noise(sample_count);
+
+            auto source = gr::blocks::vector_source<gr_complex>::make(noise);
+            auto sink = gr::blocks::vector_sink<float>::make();
+
+            auto uut = gr::droneid::variance::make(window_size);
+
+            tb->connect(source, 0, uut, 0);
+            tb->connect(uut, 0, sink, 0);
+
+            tb->run();
+
+            const auto expected_values = calculate_true_variance(noise, window_size);
+            const auto & uut_output = sink->data();
+
+            for (uint32_t idx = 0; idx < expected_values.size(); idx++) {
+                const auto delta = abs(expected_values[idx] - uut_output[idx]);
+                BOOST_ASSERT(delta < max_delta);
+            }
+
+            BOOST_ASSERT_MSG(expected_values.size() == uut_output.size(),
+                             "Did not get expected number of samples back");
+
+            plot_results(expected_values, uut_output, window_size);
+
+        }
+
+    } /* namespace droneid */
+} /* namespace gr */
+

gnuradio/gr-droneid/lib/qa_variance.h

@@ -0,0 +1,45 @@
+/* -*- c++ -*- */
+/*
+ * Copyright 2022 gr-droneid author.
+ *
+ * This is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3, or (at your option)
+ * any later version.
+ *
+ * This software is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this software; see the file COPYING.  If not, write to
+ * the Free Software Foundation, Inc., 51 Franklin Street,
+ * Boston, MA 02110-1301, USA.
+ */
+
+#ifndef _QA_VARIANCE_H_
+#define _QA_VARIANCE_H_
+
+#include <cppunit/extensions/HelperMacros.h>
+#include <cppunit/TestCase.h>
+
+namespace gr {
+  namespace droneid {
+
+    class qa_variance : public CppUnit::TestCase
+    {
+    public:
+      CPPUNIT_TEST_SUITE(qa_variance);
+      CPPUNIT_TEST(t1);
+      CPPUNIT_TEST_SUITE_END();
+
+    private:
+      void t1();
+    };
+
+  } /* namespace droneid */
+} /* namespace gr */
+
+#endif /* _QA_VARIANCE_H_ */
+

gnuradio/gr-droneid/lib/variance_impl.cc

@@ -0,0 +1,96 @@
+/* -*- c++ -*- */
+/*
+ * Copyright 2022 gr-droneid author.
+ *
+ * This is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3, or (at your option)
+ * any later version.
+ *
+ * This software is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this software; see the file COPYING.  If not, write to
+ * the Free Software Foundation, Inc., 51 Franklin Street,
+ * Boston, MA 02110-1301, USA.
+ */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <gnuradio/io_signature.h>
+#include "variance_impl.h"
+
+#include <volk/volk.h>
+
+namespace gr {
+    namespace droneid {
+
+        variance::sptr
+        variance::make(uint32_t window_size) {
+            return gnuradio::get_initial_sptr
+                    (new variance_impl(window_size));
+        }
+
+        /*
+         * The private constructor
+         */
+        variance_impl::variance_impl(uint32_t window_size)
+                : gr::block("variance",
+                                 gr::io_signature::make(1, 1, sizeof(gr_complex)),
+                                 gr::io_signature::make(1, 1, sizeof(float))),
+                                 window_size_(window_size),
+                                 window_size_recip_(1.0f / static_cast<float>(window_size)) {
+
+        }
+
+        /*
+         * Our virtual destructor.
+         */
+        variance_impl::~variance_impl() {
+        }
+
+        int variance_impl::general_work(int noutput_items, gr_vector_int &ninput_items,
+                                        gr_vector_const_void_star &input_items, gr_vector_void_star &output_items) {
+            const auto * in = (gr_complex *) input_items[0];
+            auto * out = (float *) output_items[0];
+
+            buffer_.insert(buffer_.end(), in, in + noutput_items);
+
+            if (buffer_.size() < window_size_) {
+                consume_each(noutput_items);
+                return 0;
+            }
+
+            const auto num_steps = buffer_.size() - window_size_;
+
+            if (abs_squared_.size() < buffer_.size()) {
+                abs_squared_.resize(buffer_.size());
+                sum_.resize(num_steps);
+            }
+
+            volk_32fc_magnitude_squared_32f(&abs_squared_[0], &buffer_[0], buffer_.size());
+            float running_sum;
+            volk_32f_accumulator_s32f_a(&running_sum, &abs_squared_[0], window_size_);
+            sum_[0] = running_sum;
+
+            for (uint32_t idx = 1; idx < num_steps; idx++) {
+                running_sum = running_sum - abs_squared_[idx - 1] + abs_squared_[idx + window_size_];
+                sum_[idx] = running_sum;
+            }
+
+            volk_32f_s32f_multiply_32f(out, &sum_[0], window_size_recip_, num_steps);
+
+            buffer_.erase(buffer_.begin(), buffer_.begin() + num_steps);
+
+            consume_each(noutput_items);
+            return num_steps;
+        }
+
+    } /* namespace droneid */
+} /* namespace gr */
+

gnuradio/gr-droneid/lib/variance_impl.h

@@ -0,0 +1,52 @@
+/* -*- c++ -*- */
+/*
+ * Copyright 2022 gr-droneid author.
+ *
+ * This is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3, or (at your option)
+ * any later version.
+ *
+ * This software is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this software; see the file COPYING.  If not, write to
+ * the Free Software Foundation, Inc., 51 Franklin Street,
+ * Boston, MA 02110-1301, USA.
+ */
+
+#ifndef INCLUDED_DRONEID_VARIANCE_IMPL_H
+#define INCLUDED_DRONEID_VARIANCE_IMPL_H
+
+#include <droneid/variance.h>
+
+namespace gr {
+    namespace droneid {
+
+        class variance_impl : public variance {
+        private:
+            const uint32_t window_size_;
+            const float window_size_recip_;
+
+            std::vector<gr_complex> buffer_;
+            std::vector<float> abs_squared_;
+            std::vector<float> sum_;
+
+        public:
+            variance_impl(uint32_t window_size);
+
+            ~variance_impl();
+
+            int general_work(int noutput_items, gr_vector_int &ninput_items, gr_vector_const_void_star &input_items,
+                             gr_vector_void_star &output_items) override;
+
+        };
+
+    } // namespace droneid
+} // namespace gr
+
+#endif /* INCLUDED_DRONEID_VARIANCE_IMPL_H */
+

gnuradio/gr-droneid/python/CMakeLists.txt

@@ -46,3 +46,5 @@ GR_ADD_TEST(qa_extractor ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/qa_ext
 GR_ADD_TEST(qa_time_sync ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/qa_time_sync.py)
 GR_ADD_TEST(qa_demodulation ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/qa_demodulation.py)
 GR_ADD_TEST(qa_normalized_xcorr_estimate ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/qa_normalized_xcorr_estimate.py)
+GR_ADD_TEST(qa_variance ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/qa_variance.py)
+GR_ADD_TEST(qa_dot_prod ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/qa_dot_prod.py)

gnuradio/gr-droneid/python/qa_dot_prod.py

@@ -0,0 +1,41 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Copyright 2022 gr-droneid author.
+#
+# This is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3, or (at your option)
+# any later version.
+#
+# This software is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this software; see the file COPYING.  If not, write to
+# the Free Software Foundation, Inc., 51 Franklin Street,
+# Boston, MA 02110-1301, USA.
+#
+
+from gnuradio import gr, gr_unittest
+from gnuradio import blocks
+import droneid_swig as droneid
+
+class qa_dot_prod(gr_unittest.TestCase):
+
+    def setUp(self):
+        self.tb = gr.top_block()
+
+    def tearDown(self):
+        self.tb = None
+
+    def test_001_t(self):
+        # set up fg
+        self.tb.run()
+        # check data
+
+
+if __name__ == '__main__':
+    gr_unittest.run(qa_dot_prod)

gnuradio/gr-droneid/python/qa_variance.py

@@ -0,0 +1,41 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Copyright 2022 gr-droneid author.
+#
+# This is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3, or (at your option)
+# any later version.
+#
+# This software is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this software; see the file COPYING.  If not, write to
+# the Free Software Foundation, Inc., 51 Franklin Street,
+# Boston, MA 02110-1301, USA.
+#
+
+from gnuradio import gr, gr_unittest
+from gnuradio import blocks
+import droneid_swig as droneid
+
+class qa_variance(gr_unittest.TestCase):
+
+    def setUp(self):
+        self.tb = gr.top_block()
+
+    def tearDown(self):
+        self.tb = None
+
+    def test_001_t(self):
+        # set up fg
+        self.tb.run()
+        # check data
+
+
+if __name__ == '__main__':
+    gr_unittest.run(qa_variance)

gnuradio/gr-droneid/swig/droneid_swig.i

@@ -17,6 +17,8 @@
 #include "droneid/decode.h"
 #include "droneid/normalized_xcorr.h"
 #include "droneid/normalized_xcorr_estimate.h"
+#include "droneid/variance.h"
+#include "droneid/dot_prod.h"
 //#include "droneid/utils.h"
 %}
 
@@ -40,3 +42,7 @@ GR_SWIG_BLOCK_MAGIC2(droneid, decode);
 %include "droneid/normalized_xcorr.h"
 %include "droneid/normalized_xcorr_estimate.h"
 GR_SWIG_BLOCK_MAGIC2(droneid, normalized_xcorr_estimate);
+%include "droneid/variance.h"
+GR_SWIG_BLOCK_MAGIC2(droneid, variance);
+%include "droneid/dot_prod.h"
+GR_SWIG_BLOCK_MAGIC2(droneid, dot_prod);