diff --git a/grc/CMakeLists.txt b/grc/CMakeLists.txt index 539de7d52657c7589a40a51a1d32a68ebfac11dd..ab1f4ef9a53ea89a3e3b4c1eebfce6cd76dee255 100644 --- a/grc/CMakeLists.txt +++ b/grc/CMakeLists.txt @@ -7,5 +7,5 @@ # install(FILES - DESTINATION share/gnuradio/grc/blocks + elen90089_corr_est_cc.block.yml DESTINATION share/gnuradio/grc/blocks ) diff --git a/grc/elen90089_corr_est_cc.block.yml b/grc/elen90089_corr_est_cc.block.yml new file mode 100644 index 0000000000000000000000000000000000000000..fa8de8cae93f145d501df25024527a491300fe1d --- /dev/null +++ b/grc/elen90089_corr_est_cc.block.yml @@ -0,0 +1,31 @@ +id: elen90089_corr_est_cc +label: CDC TS Correlator +category: '[elen90089]' + +templates: + imports: import elen90089 + make: elen90089.corr_est_cc(${sequence}, ${threshold}) + +parameters: +- id: sequence + label: Sequence + dtype: complex_vector +- id: threshold + label: Threshold + default: '0.5' + dtype: float + +inputs: +- label: in + domain: stream + dtype: complex +outputs: +- label: out + domain: stream + dtype: complex +- label: out + domain: stream + dtype: float + optional: True + +file_format: 1 diff --git a/include/elen90089/CMakeLists.txt b/include/elen90089/CMakeLists.txt index c546585f9092d2e0d5c98b97cc5343524298f554..1091fdada7bfb1c50398c98ebbfbee26a719471b 100644 --- a/include/elen90089/CMakeLists.txt +++ b/include/elen90089/CMakeLists.txt @@ -11,5 +11,6 @@ ######################################################################## install(FILES api.h + corr_est_cc.h DESTINATION include/elen90089 ) diff --git a/include/elen90089/corr_est_cc.h b/include/elen90089/corr_est_cc.h new file mode 100644 index 0000000000000000000000000000000000000000..95870f4af654ea1af6af7f706c5ac0e69b643885 --- /dev/null +++ b/include/elen90089/corr_est_cc.h @@ -0,0 +1,60 @@ +/* -*- c++ -*- */ +/* + * Copyright 2022 University of Melbourne. + * + * SPDX-License-Identifier: GPL-3.0-or-later + */ + +#ifndef INCLUDED_ELEN90089_CORR_EST_CC_H +#define INCLUDED_ELEN90089_CORR_EST_CC_H + +#include <elen90089/api.h> +#include <gnuradio/sync_block.h> + +namespace gr { +namespace elen90089 { + +/*! + * \brief Search for symbol sequence within sample stream via correlation + * \ingroup elen90089 + * + * \details + * Input: + * \li Stream of complex samples. + * + * Output: + * \li Output stream of input complex samples delayed by sequence length and + * tagged with correlation start and frequency offset estimate + * \li tag 'corr_start': start sample of identified sequence + * \li tag 'freq_est': estimated frequency offset + * + * \li Optional 2nd output stream providing advanced correlator output + */ +class ELEN90089_API corr_est_cc : virtual public gr::sync_block +{ +public: + typedef std::shared_ptr<corr_est_cc> sptr; + + /*! + * Make a block that searchers for the \p sequence vector within input + * sample stream and outputs stream tagged with correlation and frequency + * offset estimate. + * + * \param sequence Symbol sequence to correlate agains. + * \param threshold Correlation threshold to declare sequence found + * (0.0 to 1.0) + */ + static sptr make(const std::vector<gr_complex>& sequence, + float threshold=0.5); + + virtual std::vector<gr_complex> sequence() const = 0; + virtual void set_sequence(const std::vector<gr_complex>& sequence) = 0; + + virtual float threshold() const = 0; + virtual void set_threshold(float threshold) = 0; +}; + +} // namespace elen90089 +} // namespace gr + +#endif /* INCLUDED_ELEN90089_CORR_EST_CC_H */ diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index 2aaed005ef4afcd43413e5edf6b0e1c57212cd9d..4c6376a2a8f3d35df68aadcb8cee041ca7c9cc5d 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -12,6 +12,7 @@ include(GrPlatform) #define LIB_SUFFIX list(APPEND elen90089_sources + corr_est_cc_impl.cc ) set(elen90089_sources "${elen90089_sources}" PARENT_SCOPE) @@ -21,7 +22,12 @@ if(NOT elen90089_sources) endif(NOT elen90089_sources) add_library(gnuradio-elen90089 SHARED ${elen90089_sources}) -target_link_libraries(gnuradio-elen90089 gnuradio::gnuradio-runtime) + +target_link_libraries( + gnuradio-elen90089 + gnuradio::gnuradio-runtime + gnuradio-filter) + target_include_directories(gnuradio-elen90089 PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/../include> PUBLIC $<INSTALL_INTERFACE:include> diff --git a/lib/corr_est_cc_impl.cc b/lib/corr_est_cc_impl.cc new file mode 100644 index 0000000000000000000000000000000000000000..5cb25a7e5b63ecd591b01eedf81bc1d79e623a03 --- /dev/null +++ b/lib/corr_est_cc_impl.cc @@ -0,0 +1,186 @@ +/* -*- c++ -*- */ +/* + * Copyright 2022 University of Melbourne. + * + * SPDX-License-Identifier: GPL-3.0-or-later + */ + +#include <gnuradio/io_signature.h> +#include <gnuradio/math.h> +#include <volk/volk.h> +#include "corr_est_cc_impl.h" + +namespace gr { +namespace elen90089 { + +corr_est_cc::sptr corr_est_cc::make(const std::vector<gr_complex>& sequence, + float threshold) +{ + return gnuradio::make_block_sptr<corr_est_cc_impl>(sequence, threshold); +} + +corr_est_cc_impl::corr_est_cc_impl(const std::vector<gr_complex>& sequence, + float threshold) + : gr::sync_block("corr_est_cc", + gr::io_signature::make(1, 1, sizeof(gr_complex)), + gr::io_signature::make2(1, 2, sizeof(gr_complex), sizeof(float))), + d_src_id(pmt::intern(alias())), + d_sequence(sequence.size()), + d_threshold(threshold), + d_filter(1, sequence), + d_corr(s_nitems), + d_corr_mag(s_nitems), + d_y_mag(s_nitems + sequence.size()), + d_y_accum(0), + d_skip(0) +{ + set_max_noutput_items(s_nitems); + set_sequence(sequence); +} + +corr_est_cc_impl::~corr_est_cc_impl() {} + +std::vector<gr_complex> corr_est_cc_impl::sequence() const +{ + return d_sequence; +} + +void corr_est_cc_impl::set_sequence(const std::vector<gr_complex>& sequence) +{ + gr::thread::scoped_lock lock(d_setlock); + + d_sequence = sequence; + + // create time-reversed conjugate of training sequence for filtering + std::vector<gr_complex> taps(sequence.size()); + std::reverse_copy(sequence.begin(), sequence.end(), taps.begin()); + float scale = 0.0; + for (size_t i = 0; i < taps.size(); i++) { + taps[i] = std::conj(taps[i]); + scale += std::norm(taps[i]); + } + scale = 1.0 / sqrt(scale); + for (size_t i = 0; i < taps.size(); i++) { + taps[i] *= scale; + } + + // set taps and block output multiple to FFT kernel's internal nsamples + const int nsamples = d_filter.set_taps(taps); + set_output_multiple(nsamples); + + // keep a history of the length of the sync word to delay for tagging + set_history(d_sequence.size()); + declare_sample_delay(0, sequence.size()); + declare_sample_delay(1, 0); +} + +float corr_est_cc_impl::threshold() const +{ + return d_threshold; +} + +void corr_est_cc_impl::set_threshold(float threshold) +{ + gr::thread::scoped_lock lock(d_setlock); + d_threshold = threshold; +} + +double corr_est_cc_impl::estimate_freq_offset(const gr_complex* samples) +{ + double accum = 0; + + gr_complex y_p = samples[0] / d_sequence[0]; + for (int i = 1; i < d_sequence.size(); i++) { + gr_complex y = samples[i] / d_sequence[i]; + gr_complex p = y * conj(y_p); + accum += fast_atan2f(p.imag(), p.real()); + y_p = y; + } + + return (accum / (d_sequence.size() - 1)); +} + +int corr_est_cc_impl::work(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) +{ + auto in = static_cast<const gr_complex*>(input_items[0]); + auto out = static_cast<gr_complex*>(output_items[0]); + + float* corr_mag; + if (output_items.size() > 1) { + corr_mag = static_cast<float*>(output_items[1]); + } else { + corr_mag = d_corr_mag.data(); + } + + unsigned int hist_len = history() - 1; + + float threshold; + { + gr::thread::scoped_lock lock(d_setlock); + threshold = d_threshold; + } + + // correlate input samples with training sequence + d_filter.filter(noutput_items, &in[hist_len], d_corr.data()); + + // calculate squared magnitude of correlation result + volk_32fc_magnitude_squared_32f(corr_mag, d_corr.data(), noutput_items); + + // calculate squared magnitued of input samples for normalization + volk_32fc_magnitude_squared_32f(d_y_mag.data(), in, noutput_items + history()); + + // search for correlation peaks above threshold + for (int i = 0; i < noutput_items; i++) { + // add sample to moving average + d_y_accum += d_y_mag[hist_len + i]; + + // normalize cross-correlation + corr_mag[i] /= d_y_accum; + + // drop oldest sample from moving average + d_y_accum -= d_y_mag[i]; + + if (d_skip > 0) { + d_skip--; + continue; + } + + if (corr_mag[i] > d_threshold) { + // coarse frequency estimate + double freq_est = estimate_freq_offset(&in[i]); + + // tag output + add_item_tag(0, nitems_written(0) + i, + pmt::intern("corr_start"), + pmt::from_double(corr_mag[i]), + d_src_id); + add_item_tag(0, nitems_written(0) + i, + pmt::intern("freq_est"), + pmt::from_double(freq_est), + d_src_id); + if (output_items.size() > 1) { + add_item_tag(1, nitems_written(0) + i, + pmt::intern("corr_start"), + pmt::from_double(corr_mag[i]), + d_src_id); + add_item_tag(1, nitems_written(0) + i, + pmt::intern("freq_est"), + pmt::from_double(freq_est), + d_src_id); + } + + // skip remaining symbols in sequence + d_skip = d_sequence.size() - 1; + } + } + + // pass through input samples but with delay of history() + memcpy(out, &in[0], sizeof(gr_complex)*noutput_items); + + return noutput_items; +} + +} /* namespace elen90089 */ +} /* namespace gr */ diff --git a/lib/corr_est_cc_impl.h b/lib/corr_est_cc_impl.h new file mode 100644 index 0000000000000000000000000000000000000000..9fc44e640ea25a0f498340ced0f26d79080c3a2b --- /dev/null +++ b/lib/corr_est_cc_impl.h @@ -0,0 +1,58 @@ +/* -*- c++ -*- */ +/* + * Copyright 2022 University of Melbourne. + * + * SPDX-License-Identifier: GPL-3.0-or-later + */ + +#ifndef INCLUDED_ELEN90089_CORR_EST_CC_IMPL_H +#define INCLUDED_ELEN90089_CORR_EST_CC_IMPL_H + +#include <elen90089/corr_est_cc.h> +#include <gnuradio/filter/fft_filter.h> + +using namespace gr::filter; + +namespace gr { +namespace elen90089 { + +class corr_est_cc_impl : public corr_est_cc +{ +private: + pmt::pmt_t d_src_id; + + std::vector<gr_complex> d_sequence; + float d_threshold; + + kernel::fft_filter_ccc d_filter; + volk::vector<gr_complex> d_corr; + volk::vector<float> d_corr_mag; + volk::vector<float> d_y_mag; + float d_y_accum; + int d_skip; + + static constexpr int s_nitems = 24*1024; + + double estimate_freq_offset(const gr_complex* samples); + +public: + corr_est_cc_impl(const std::vector<gr_complex>& sequence, + float threshold = 0.5); + + ~corr_est_cc_impl() override; + + std::vector<gr_complex> sequence() const override; + void set_sequence(const std::vector<gr_complex>& sequence) override; + + float threshold() const override; + void set_threshold(float threshold) override; + + int work(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) override; +}; + +} // namespace elen90089 +} // namespace gr + +#endif /* INCLUDED_ELEN90089_CORR_EST_CC_IMPL_H */ diff --git a/python/bindings/CMakeLists.txt b/python/bindings/CMakeLists.txt index 4398c7122a04df231e10163d4d949bf64b4ab26e..3c360f53dc91f633be59d0ff3af5286ab08d87f6 100644 --- a/python/bindings/CMakeLists.txt +++ b/python/bindings/CMakeLists.txt @@ -29,7 +29,7 @@ include(GrPybind) ######################################################################## list(APPEND elen90089_python_files - python_bindings.cc) + corr_est_cc_python.cc python_bindings.cc) GR_PYBIND_MAKE_OOT(elen90089 ../.. diff --git a/python/bindings/corr_est_cc_python.cc b/python/bindings/corr_est_cc_python.cc new file mode 100644 index 0000000000000000000000000000000000000000..e602ccef753146f0647d14ea8b151868f41d7d09 --- /dev/null +++ b/python/bindings/corr_est_cc_python.cc @@ -0,0 +1,64 @@ +/* + * Copyright 2022 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * SPDX-License-Identifier: GPL-3.0-or-later + * + */ + +/***********************************************************************************/ +/* This file is automatically generated using bindtool and can be manually edited */ +/* The following lines can be configured to regenerate this file during cmake */ +/* If manual edits are made, the following tags should be modified accordingly. */ +/* BINDTOOL_GEN_AUTOMATIC(0) */ +/* BINDTOOL_USE_PYGCCXML(0) */ +/* BINDTOOL_HEADER_FILE(corr_est_cc.h) */ +/* BINDTOOL_HEADER_FILE_HASH(53737d29a95e3a8c008b00324c083e88) */ +/***********************************************************************************/ + +#include <pybind11/complex.h> +#include <pybind11/pybind11.h> +#include <pybind11/stl.h> + +namespace py = pybind11; + +#include <elen90089/corr_est_cc.h> +// pydoc.h is automatically generated in the build directory +#include <corr_est_cc_pydoc.h> + +void bind_corr_est_cc(py::module& m) +{ + using corr_est_cc = ::gr::elen90089::corr_est_cc; + + py::class_<corr_est_cc, + gr::sync_block, + gr::block, + gr::basic_block, + std::shared_ptr<corr_est_cc>>(m, "corr_est_cc", D(corr_est_cc)) + + .def(py::init(&corr_est_cc::make), + py::arg("sequence"), + py::arg("threshold") = 0.5, + D(corr_est_cc, make)) + + .def("sequence", + &corr_est_cc::sequence, + D(corr_est_cc, sequence)) + + .def("set_sequence", + &corr_est_cc::set_sequence, + py::arg("sequence"), + D(corr_est_cc, set_sequence)) + + .def("threshold", + &corr_est_cc::threshold, + D(corr_est_cc, threshold)) + + .def("set_threshold", + &corr_est_cc::set_threshold, + py::arg("threshold"), + D(corr_est_cc, set_threshold)); + +} + diff --git a/python/bindings/docstrings/corr_est_cc_pydoc_template.h b/python/bindings/docstrings/corr_est_cc_pydoc_template.h new file mode 100644 index 0000000000000000000000000000000000000000..94f560f8f1e0b96938824787b9d41777c6cd18ce --- /dev/null +++ b/python/bindings/docstrings/corr_est_cc_pydoc_template.h @@ -0,0 +1,42 @@ +/* + * Copyright 2022 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * SPDX-License-Identifier: GPL-3.0-or-later + * + */ +#include "pydoc_macros.h" +#define D(...) DOC(gr,elen90089, __VA_ARGS__ ) +/* + This file contains placeholders for docstrings for the Python bindings. + Do not edit! These were automatically extracted during the binding process + and will be overwritten during the build process + */ + + + + static const char *__doc_gr_elen90089_corr_est_cc = R"doc()doc"; + + + static const char *__doc_gr_elen90089_corr_est_cc_corr_est_cc_0 = R"doc()doc"; + + + static const char *__doc_gr_elen90089_corr_est_cc_corr_est_cc_1 = R"doc()doc"; + + + static const char *__doc_gr_elen90089_corr_est_cc_make = R"doc()doc"; + + + static const char *__doc_gr_elen90089_corr_est_cc_sequence = R"doc()doc"; + + + static const char *__doc_gr_elen90089_corr_est_cc_set_sequence = R"doc()doc"; + + + static const char *__doc_gr_elen90089_corr_est_cc_threshold = R"doc()doc"; + + + static const char *__doc_gr_elen90089_corr_est_cc_set_threshold = R"doc()doc"; + + diff --git a/python/bindings/python_bindings.cc b/python/bindings/python_bindings.cc index ea20a8e31e850105b46bc0d5dded22a9a2d87bd9..f3152253fbc8e008943fb9c9d006597da61a4cce 100644 --- a/python/bindings/python_bindings.cc +++ b/python/bindings/python_bindings.cc @@ -21,6 +21,7 @@ namespace py = pybind11; // Please do not delete /**************************************/ // BINDING_FUNCTION_PROTOTYPES( + void bind_corr_est_cc(py::module& m); // ) END BINDING_FUNCTION_PROTOTYPES @@ -49,5 +50,6 @@ PYBIND11_MODULE(elen90089_python, m) // Please do not delete /**************************************/ // BINDING_FUNCTION_CALLS( + bind_corr_est_cc(m); // ) END BINDING_FUNCTION_CALLS -} +} \ No newline at end of file diff --git a/python/qa_corr_est_cc.py b/python/qa_corr_est_cc.py new file mode 100644 index 0000000000000000000000000000000000000000..d9bfa1fdef74197705265c988acddc5604ad8560 --- /dev/null +++ b/python/qa_corr_est_cc.py @@ -0,0 +1,63 @@ +from gnuradio import gr, gr_unittest +from gnuradio import blocks +import pmt +import elen90089 +import numpy as np + + +class qa_corr_est_cc(gr_unittest.TestCase): + + def setUp(self): + self.tb = gr.top_block() + + def tearDown(self): + self.tb = None + + def test_001(self): + """Test correlation result""" + pad = (0,) * 16 + code = (1, 0, 1, 1, 1, 0, 1, 1) + code_sym = [(2.0*x - 1.0) for x in code] + src_bits = pad + code + pad + src_data = [(2.0*x - 1.0) for x in src_bits] + # build flowgraph + src = blocks.vector_source_c(src_data) + cor = elen90089.corr_est_cc(code_sym, threshold=0.5) + dst = blocks.tag_debug(gr.sizeof_gr_complex, "", "corr_start") + dst.set_display(False) + self.tb.connect(src, cor, dst) + self.tb.run() + # check test results + result_data = dst.current_tags() + index = len(pad) + len(code) - 1 + self.assertEqual(len(result_data), 1) + self.assertEqual(result_data[0].offset, index) + corr_mag = pmt.to_python(result_data[0].value) + self.assertAlmostEqual(corr_mag, 1.0, places=6) + + def test_002(self): + """Test coarse frequency offset estimate.""" + phase_inc = 2*np.pi*0.02 # freq offset + pad = (0,) * 16 + code = (1, 0, 1, 1, 1, 0, 1, 1) + code_sym = [(2.0*x - 1.0) for x in code] + src_bits = pad + code + pad + src_data = [(2.0*x - 1.0) for x in src_bits] + # build flowgraph + src = blocks.vector_source_c(src_data) + rot = blocks.rotator_cc(phase_inc) + cor = elen90089.corr_est_cc(code_sym, threshold=0.5) + dst = blocks.tag_debug(gr.sizeof_gr_complex, "", "freq_est") + dst.set_display(False) + self.tb.connect(src, rot, cor, dst) + self.tb.run() + # check test results + result_data = dst.current_tags() + index = len(pad) + len(code) - 1 + self.assertEqual(len(result_data), 1) + self.assertEqual(result_data[0].offset, index) + freq_est = pmt.to_python(result_data[0].value) + self.assertAlmostEqual(freq_est, phase_inc, places=6) + +if __name__ == '__main__': + gr_unittest.run(qa_corr_est_cc)