mirror
/
dji_droneid
kopia lustrzana https://github.com/proto17/dji_droneid
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność

Moved dot_prod to normalized_xcorr_estimate 

The dot_prod module ended up being a full on normalized cross correlation module, so it made sense to replace the old normalized xcorr module
gr-droneid-update
David Protzman 2022-06-26 10:24:05 -04:00
rodzic 85f39daf0b
commit df985ba877
19 zmienionych plików z 337 dodań i 581 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

3
gnuradio/gr-droneid/grc/CMakeLists.txt
Wyświetl plik

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

42
gnuradio/gr-droneid/grc/droneid_dot_prod.block.yml
Wyświetl plik

@ -1,42 +0,0 @@
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

6
gnuradio/gr-droneid/grc/droneid_normalized_xcorr_estimate.block.yml
Wyświetl plik

@ -4,7 +4,7 @@ category: '[droneid]'
templates:
imports: import droneid
make: droneid.normalized_xcorr_estimate(${filter_taps})
make: droneid.normalized_xcorr_estimate(${taps})
# Make one 'parameters' list entry for every parameter you want settable from the GUI.
# Keys include:
@ -12,7 +12,7 @@ templates:
# * label (label shown in the GUI)
# * dtype (e.g. int, float, complex, byte, short, xxx_vector, ...)
parameters:
- id: filter_taps
- id: taps
label: Filter Taps
dtype: complex_vector
@ -33,7 +33,7 @@ inputs:
outputs:
- label: out
domain: stream
dtype: float
dtype: complex
vlen: 1
optional: 0

3
gnuradio/gr-droneid/include/droneid/CMakeLists.txt
Wyświetl plik

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

54
gnuradio/gr-droneid/include/droneid/dot_prod.h
Wyświetl plik

@ -1,54 +0,0 @@
/* -*- 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 */

43
gnuradio/gr-droneid/include/droneid/normalized_xcorr_estimate.h
Wyświetl plik

@ -22,33 +22,32 @@
#define INCLUDED_DRONEID_NORMALIZED_XCORR_ESTIMATE_H
#include <droneid/api.h>
#include <gnuradio/sync_block.h>
#include <gnuradio/block.h>
namespace gr {
namespace droneid {
namespace droneid {
/*!
* \brief <+description of block+>
* \ingroup droneid
*
*/
class DRONEID_API normalized_xcorr_estimate : virtual public gr::sync_block
{
public:
typedef boost::shared_ptr<normalized_xcorr_estimate> sptr;
/*!
* \brief <+description of block+>
* \ingroup droneid
*
*/
class DRONEID_API normalized_xcorr_estimate : virtual public gr::block {
public:
typedef boost::shared_ptr<normalized_xcorr_estimate> sptr;
/*!
* \brief Return a shared_ptr to a new instance of droneid::normalized_xcorr_estimate.
*
* To avoid accidental use of raw pointers, droneid::normalized_xcorr_estimate's
* constructor is in a private implementation
* class. droneid::normalized_xcorr_estimate::make is the public interface for
* creating new instances.
*/
static sptr make(std::vector<std::complex<float>> filter_taps);
};
/*!
* \brief Return a shared_ptr to a new instance of droneid::normalized_xcorr_estimate.
*
* To avoid accidental use of raw pointers, droneid::normalized_xcorr_estimate's
* constructor is in a private implementation
* class. droneid::normalized_xcorr_estimate::make is the public interface for
* creating new instances.
*/
static sptr make(const std::vector<gr_complex> & /*taps*/);
};
} // namespace droneid
} // namespace droneid
} // namespace gr
#endif /* INCLUDED_DRONEID_NORMALIZED_XCORR_ESTIMATE_H */

2
gnuradio/gr-droneid/lib/CMakeLists.txt
Wyświetl plik

@ -34,7 +34,6 @@ list(APPEND droneid_sources
normalized_xcorr.cc
normalized_xcorr_estimate_impl.cc
variance_impl.cc
dot_prod_impl.cc
)
set(droneid_sources "${droneid_sources}" PARENT_SCOPE)
@ -102,7 +101,6 @@ if (MATLAB_PATH)
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")

179
gnuradio/gr-droneid/lib/dot_prod_impl.cc
Wyświetl plik

@ -1,179 +0,0 @@
/* -*- 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 */

61
gnuradio/gr-droneid/lib/dot_prod_impl.h
Wyświetl plik

@ -1,61 +0,0 @@
/* -*- 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 */

137
gnuradio/gr-droneid/lib/normalized_xcorr_estimate_impl.cc
Wyświetl plik

@ -24,30 +24,42 @@
#include <gnuradio/io_signature.h>
#include "normalized_xcorr_estimate_impl.h"
#include <volk/volk.h>
#include <numeric>
#include <droneid/misc_utils.h>
namespace gr {
namespace droneid {
normalized_xcorr_estimate::sptr
normalized_xcorr_estimate::make(std::vector<std::complex<float>> filter_taps) {
normalized_xcorr_estimate::make(const std::vector<gr_complex> &taps) {
return gnuradio::get_initial_sptr
(new normalized_xcorr_estimate_impl(filter_taps));
(new normalized_xcorr_estimate_impl(taps));
}
/*
* The private constructor
*/
normalized_xcorr_estimate_impl::normalized_xcorr_estimate_impl(std::vector<std::complex<float>> filter_taps)
: gr::sync_block("normalized_xcorr_estimate",
gr::io_signature::make(1, 1, sizeof(gr_complex)),
gr::io_signature::make(1, 1, sizeof(float))) {
xcorr_ = std::unique_ptr<normalized_xcorr>(new normalized_xcorr(filter_taps));
normalized_xcorr_estimate_impl::normalized_xcorr_estimate_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()) {
set_history(filter_taps.size());
set_alignment(std::max(1, static_cast<int32_t>(volk_get_alignment() / sizeof(float))));
// 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};
}
/*
@ -57,18 +69,109 @@ namespace gr {
}
int
normalized_xcorr_estimate_impl::work(int noutput_items,
gr_vector_const_void_star &input_items,
gr_vector_void_star &output_items) {
normalized_xcorr_estimate_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 = (float *) output_items[0];
auto *out = (gr_complex *) output_items[0];
xcorr_->run(in, noutput_items, out);
consume_each(noutput_items);
// Do <+signal processing+>
// 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 noutput_items;
return num_steps;
}
} /* namespace droneid */

25
gnuradio/gr-droneid/lib/normalized_xcorr_estimate_impl.h
Wyświetl plik

@ -22,26 +22,33 @@
#define INCLUDED_DRONEID_NORMALIZED_XCORR_ESTIMATE_IMPL_H
#include <droneid/normalized_xcorr_estimate.h>
#include <droneid/normalized_xcorr.h>
namespace gr {
namespace droneid {
class normalized_xcorr_estimate_impl : public normalized_xcorr_estimate {
private:
std::unique_ptr<gr::droneid::normalized_xcorr> xcorr_;
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:
normalized_xcorr_estimate_impl(std::vector<std::complex<float>> filter_taps);
normalized_xcorr_estimate_impl(const std::vector<gr_complex> & taps);
~normalized_xcorr_estimate_impl();
// Where all the action really happens
int work(
int noutput_items,
gr_vector_const_void_star &input_items,
gr_vector_void_star &output_items
);
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

64
gnuradio/gr-droneid/lib/qa_dot_prod.cc
Wyświetl plik

@ -1,64 +0,0 @@
/* -*- 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 */

45
gnuradio/gr-droneid/lib/qa_dot_prod.h
Wyświetl plik

@ -1,45 +0,0 @@
/* -*- 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_ */

174
gnuradio/gr-droneid/lib/qa_normalized_xcorr.cc
Wyświetl plik

@ -26,49 +26,167 @@
#include <boost/test/unit_test.hpp>
#include <gnuradio/random.h>
#include <chrono>
#include <type_traits>
#include <droneid/misc_utils.h>
namespace gr {
namespace droneid {
using sample_t = float;
using complex_t = std::complex<sample_t>;
using complex_vec_t = std::vector<complex_t>;
#include <MatlabEngine.hpp>
#include <MatlabDataArray.hpp>
BOOST_AUTO_TEST_CASE(foo) {
auto rng = gr::random();
complex_vec_t noise_vector(1e6, {0, 0});
for (auto & sample : noise_vector) {
sample = rng.rayleigh_complex();
using namespace matlab::engine;
namespace gr {
namespace droneid {
using sample_t = float;
using complex_t = std::complex<sample_t>;
using complex_vec_t = std::vector<complex_t>;
template <typename OUTPUT_T>
std::vector<OUTPUT_T> run(const std::u16string & cmd, std::vector<matlab::data::Array> & inputs, MATLABEngine * engine) {
std::vector<OUTPUT_T> output;
matlab::data::Array matlab_output = engine->feval(cmd, inputs);
output.resize(matlab_output.getNumberOfElements());
for (uint32_t idx = 0; idx < matlab_output.getNumberOfElements(); idx++) {
output[idx] = static_cast<OUTPUT_T>(matlab_output[idx]);
}
return output;
}
const auto start_offset = 102313; // Just a random starting offset that's not on an even boundary
const auto filter_length = 1023;
template <typename OUTPUT_T>
std::vector<OUTPUT_T> get_vec(const std::u16string & variable_name, MATLABEngine * const engine_ptr) {
std::vector<OUTPUT_T> output;
const auto filter_taps = complex_vec_t(noise_vector.begin() + start_offset, noise_vector.begin() + start_offset + filter_length);
auto variable_value = engine_ptr->getVariable(variable_name);
output.reserve(variable_value.getNumberOfElements());
for (const auto & sample : matlab::data::getReadOnlyElements<OUTPUT_T>(variable_value)) {
output.push_back(sample);
}
return output;
}
template <typename OUTPUT_T>
std::vector<std::complex<OUTPUT_T>> get_complex_vec(const std::u16string & variable_name, MATLABEngine * const engine_ptr) {
std::vector<std::complex<OUTPUT_T>> output;
auto variable_value = engine_ptr->getVariable(variable_name);
output.reserve(variable_value.getNumberOfElements());
for (const std::complex<double> & sample : matlab::data::getReadOnlyElements<std::complex<double>>(variable_value)) {
output.push_back({
static_cast<OUTPUT_T>(sample.real()), static_cast<OUTPUT_T>(sample.imag())
});
}
return output;
}
uint64_t get_matrix_length(const std::u16string & variable_name, MATLABEngine * const engine_ptr) {
const auto cmd = u"length(" + variable_name + u");";
engine_ptr->eval(cmd);
const auto ret = engine_ptr->getVariable("ans");
return static_cast<uint64_t>(ret[0]);
}
BOOST_AUTO_TEST_CASE(foo) {
auto matlab_ptr = startMATLAB();
matlab::data::ArrayFactory factory;
normalized_xcorr xcorr(filter_taps);
gr::random random;
const auto burst = misc_utils::read_samples("/tmp/droneid_debug/burst_25", 0, 0);
std::vector<sample_t> mags(noise_vector.size() - filter_length, 0);
const auto full_sample_count = static_cast<uint32_t>(.1e6);
const auto start = std::chrono::high_resolution_clock::now();
xcorr.run(&noise_vector[0], noise_vector.size(), &mags[0]);
const auto end = std::chrono::high_resolution_clock::now();
const auto duration = std::chrono::duration<float>(end - start).count();
const auto rate = noise_vector.size() / duration;
std::cout << "Took " << duration << " seconds to run through " << noise_vector.size() << " samples\n";
std::cout << "Average of " << rate << " samples per second\n";
complex_vec_t samples(full_sample_count);
const int64_t padding = floor((static_cast<int64_t>(full_sample_count) - static_cast<int64_t>(burst.size())) / 2);
std::cout << "Padding: " << padding << "\n";
if (padding > 0) {
for (uint64_t idx = 0; idx < padding; idx++) {
samples[idx] = random.rayleigh_complex();
}
std::copy(burst.begin(), burst.end(), samples.begin() + padding);
for (uint64_t idx = padding + burst.size(); idx < samples.size(); idx++) {
samples[idx] = random.rayleigh_complex();
}
} else {
samples = burst;
}
auto filter_taps = misc_utils::create_zc_sequence(15.36e6, 600);
filter_taps.resize(filter_taps.size() / 1);
const auto filter_length = filter_taps.size();
misc_utils::write("/tmp/mags", &mags[0], sizeof(mags[0]), mags.size());
normalized_xcorr xcorr(filter_taps);
std::vector<sample_t> mags(samples.size() - filter_length, 0);
xcorr.run(&samples[0], samples.size(), &mags[0]);
const int iters = 10;
const auto start = std::chrono::high_resolution_clock::now();
for (int iter = 0; iter < iters; iter++) {
xcorr.run(&samples[0], samples.size(), &mags[0]);
}
const auto end = std::chrono::high_resolution_clock::now();
const auto duration = std::chrono::duration<float>(end - start).count();
const auto rate = (samples.size() * iters) / duration;
std::cout << "Took " << duration << " seconds to run through " << samples.size() << " samples\n";
std::cout << "Average of " << rate << " samples per second\n";
const auto max_element_iter = std::max_element(mags.begin(), mags.end());
const auto distance = std::distance(mags.begin(), max_element_iter);
std::cout << "Distance: " << distance << "\n";
}
const auto max_element_iter = std::max_element(mags.begin(), mags.end());
const auto distance = std::distance(mags.begin(), max_element_iter);
std::cout << "c++ Distance: " << distance << "\n";
} /* namespace droneid */
complex_vec_t matlab_golden_mags;
{
auto samples_vec = factory.createArray<complex_t>(
matlab::data::ArrayDimensions({samples.size(), 1}),
&samples[0],
&samples[samples.size() - 1]);
auto filter_taps_vec = factory.createArray<complex_t>(
matlab::data::ArrayDimensions({filter_taps.size(), 1}),
&filter_taps[0],
&filter_taps[filter_taps.size() - 1]);
matlab_ptr->setVariable("samples", samples_vec);
matlab_ptr->setVariable("filter", filter_taps_vec);
const std::u16string correlation_script =
u"scores1 = zeros(length(samples) - length(filter), 1);\n"
u"scores2 = zeros(length(samples) - length(filter), 1);\n"
u"for idx = 1:length(scores1)\n"
u" window = samples(idx:idx + length(filter) - 1);"
u" scores1(idx) = xcorr(window, filter, 'normalized', 0);"
u" scores2(idx) = sum(window .* conj(filter)) / length(filter);"
u"end";
matlab_ptr->eval(correlation_script);
matlab_golden_mags = get_complex_vec<sample_t>(u"scores1", matlab_ptr.get());
std::cout << "Got back " << matlab_golden_mags.size() << " correlation scores\n";
}
const auto matlab_golden_mags_sqrd = misc_utils::abs_squared(matlab_golden_mags);
matlab_ptr->setVariable(u"cpp_mags", factory.createArray<sample_t>(
matlab::data::ArrayDimensions({mags.size(), 1}), &mags[0], &mags[mags.size() - 1]));
matlab_ptr->eval(u"delta = cpp_mags - (abs(scores1).^2);");
const auto deltas = get_vec<sample_t>(u"delta", matlab_ptr.get());
matlab_ptr->eval(u"figure(1); plot(delta);");
matlab_ptr->eval(u"figure(2); subplot(3, 1, 1); plot(abs(cpp_mags)); title('CPP mags');");
matlab_ptr->eval(u"figure(2); subplot(3, 1, 2); plot(abs(scores1).^2); title('MATLAB mags');");
matlab_ptr->eval(u"figure(2); subplot(3, 1, 3); plot(10 * log10(abs(samples).^2)); title('Raw Samples');");
matlab_ptr->eval(u"pause;");
}
} /* namespace droneid */
} /* namespace gr */

33
gnuradio/gr-droneid/lib/qa_normalized_xcorr_estimate.cc
Wyświetl plik

@ -18,20 +18,43 @@
* Boston, MA 02110-1301, USA.
*/
#include <iostream>
#include <gnuradio/attributes.h>
#include <cppunit/TestAssert.h>
#include "qa_normalized_xcorr_estimate.h"
#include <droneid/normalized_xcorr_estimate.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 {
namespace droneid {
BOOST_AUTO_TEST_CASE(normalized_xcorr_estimate_test) {
auto tb = gr::make_top_block("top");
BOOST_AUTO_TEST_CASE(foo) {
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();
} /* namespace droneid */
auto uut = droneid::normalized_xcorr_estimate::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 */

1
gnuradio/gr-droneid/python/CMakeLists.txt
Wyświetl plik

@ -47,4 +47,3 @@ GR_ADD_TEST(qa_time_sync ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/qa_tim
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)

41
gnuradio/gr-droneid/python/qa_dot_prod.py
Wyświetl plik

@ -1,41 +0,0 @@
#!/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)

2
gnuradio/gr-droneid/python/qa_normalized_xcorr_estimate.py
Wyświetl plik

@ -21,7 +21,7 @@
from gnuradio import gr, gr_unittest
from gnuradio import blocks
import droneid_swig as droneid
import normalized_xcorr_estimate_swig as droneid
class qa_normalized_xcorr_estimate(gr_unittest.TestCase):

3
gnuradio/gr-droneid/swig/droneid_swig.i
Wyświetl plik

@ -18,7 +18,6 @@
#include "droneid/normalized_xcorr.h"
#include "droneid/normalized_xcorr_estimate.h"
#include "droneid/variance.h"
#include "droneid/dot_prod.h"
//#include "droneid/utils.h"
%}
@ -44,5 +43,3 @@ GR_SWIG_BLOCK_MAGIC2(droneid, decode);
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);