Added variance block

gnuradio/gr-droneid/grc/CMakeLists.txt

@@ -25,5 +25,6 @@ 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 DESTINATION share/gnuradio/grc/blocks
 )

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,6 @@ install(FILES
     lte_decode.h
     decode.h
     normalized_xcorr.h
-    normalized_xcorr_estimate.h DESTINATION include/droneid
+    normalized_xcorr_estimate.h
+    variance.h DESTINATION include/droneid
 )

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,7 @@ list(APPEND droneid_sources
     decode_impl.cc
     normalized_xcorr.cc
     normalized_xcorr_estimate_impl.cc
+    variance_impl.cc
 )
 
 set(droneid_sources "${droneid_sources}" PARENT_SCOPE)
@@ -82,6 +83,7 @@ include(GrTest)
 list(APPEND test_droneid_sources
         qa_normalized_xcorr.cc
         qa_normalized_xcorr_estimate.cc
+        qa_variance.cc
 )
 # Anything we need to link to for the unit tests go here
 list(APPEND GR_TEST_TARGET_DEPS gnuradio-droneid)

gnuradio/gr-droneid/lib/qa_variance.cc

@@ -0,0 +1,151 @@
+/* -*- 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;
+        }
+
+        BOOST_AUTO_TEST_CASE(foo2) {
+            auto tb = gr::make_top_block("top");
+
+            const auto sample_count = 12000;
+            const auto window_size = 1024;
+
+            const auto noise = create_noise(sample_count);
+
+//            std::vector<gr_complex> kk(sample_count);
+//            for (uint32_t idx = 0; idx < sample_count; idx++) {
+//                kk[idx] = {(float)idx, (float)idx};
+//            }
+
+            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();
+
+//            tb->start();
+//
+//            uint32_t last_sample_count = 0;
+//            std::this_thread::sleep_for(std::chrono::milliseconds(10));
+//            while(true) {
+//                if (last_sample_count != sink->data().size()) {
+//                    last_sample_count = sink->data().size();
+//                    std::this_thread::sleep_for(std::chrono::milliseconds(10));
+//                } else {
+//                    break;
+//                }
+//            }
+//
+//            tb->stop();
+
+            const auto expected_values = calculate_true_variance(noise, window_size);
+            const auto & uut_output = sink->data();
+
+            std::cout << "Got " << expected_values.size() << " elements from MATLAB and ";
+            std::cout << uut_output.size() << " elements from GNU Radio\n";
+
+            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");
+
+        }
+
+    } /* 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,4 @@ 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)

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,7 @@
 #include "droneid/decode.h"
 #include "droneid/normalized_xcorr.h"
 #include "droneid/normalized_xcorr_estimate.h"
+#include "droneid/variance.h"
 //#include "droneid/utils.h"
 %}
 
@@ -40,3 +41,5 @@ 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);