recentering of filter taps

CMakeLists.txt

@@ -158,6 +158,9 @@ if(NOT CPPUNIT_FOUND)
     message(FATAL_ERROR "CppUnit required to compile digitalhf")
 endif()
 
+include(GrMiscUtils)
+GR_LOGGING()
+
 ########################################################################
 # Setup doxygen option
 ########################################################################

examples/test_cc.grc

@@ -317,7 +317,7 @@
     </param>
     <param>
       <key>value</key>
-      <value>0.002</value>
+      <value>0.001</value>
     </param>
     <param>
       <key>_enabled</key>
@@ -353,15 +353,15 @@
     </param>
     <param>
       <key>start</key>
-      <value>0</value>
+      <value>0.0001</value>
     </param>
     <param>
       <key>step</key>
-      <value>0.001</value>
+      <value>0.0001</value>
     </param>
     <param>
       <key>stop</key>
-      <value>0.02</value>
+      <value>0.01</value>
     </param>
     <param>
       <key>rangeType</key>
@@ -396,7 +396,7 @@
     </param>
     <param>
       <key>value</key>
-      <value>15</value>
+      <value>9</value>
     </param>
   </block>
   <block>
@@ -423,7 +423,7 @@
     </param>
     <param>
       <key>value</key>
-      <value>15</value>
+      <value>9</value>
     </param>
   </block>
   <block>
@@ -453,6 +453,49 @@
       <value>4</value>
     </param>
   </block>
+  <block>
+    <key>variable_rrc_filter_taps</key>
+    <param>
+      <key>comment</key>
+      <value></value>
+    </param>
+    <param>
+      <key>_enabled</key>
+      <value>True</value>
+    </param>
+    <param>
+      <key>alpha</key>
+      <value>0.25</value>
+    </param>
+    <param>
+      <key>_coordinate</key>
+      <value>(469, 352)</value>
+    </param>
+    <param>
+      <key>_rotation</key>
+      <value>0</value>
+    </param>
+    <param>
+      <key>gain</key>
+      <value>1.0</value>
+    </param>
+    <param>
+      <key>id</key>
+      <value>rrc_taps</value>
+    </param>
+    <param>
+      <key>ntaps</key>
+      <value>11*sps</value>
+    </param>
+    <param>
+      <key>samp_rate</key>
+      <value>samp_rate</value>
+    </param>
+    <param>
+      <key>sym_rate</key>
+      <value>samp_rate/sps</value>
+    </param>
+  </block>
   <block>
     <key>variable</key>
     <param>
@@ -477,7 +520,7 @@
     </param>
     <param>
       <key>value</key>
-      <value>12000</value>
+      <value>12001</value>
     </param>
   </block>
   <block>
@@ -774,7 +817,7 @@
     </param>
     <param>
       <key>file</key>
-      <value>/Users/chm/Downloads/julusdalen.proxy.kiwisdr.com_2018-10-28T17_36_28Z_4644.80_iq.wav</value>
+      <value>/Users/chm/Downloads/kiwi-sm2gct.mooo.com_2018-10-30T10_35_39Z_8700.00_iq.wav</value>
     </param>
     <param>
       <key>_coordinate</key>
@@ -876,6 +919,112 @@
       <value>STANAG_4285</value>
     </param>
   </block>
+  <block>
+    <key>fir_filter_xxx</key>
+    <param>
+      <key>alias</key>
+      <value></value>
+    </param>
+    <param>
+      <key>comment</key>
+      <value></value>
+    </param>
+    <param>
+      <key>affinity</key>
+      <value></value>
+    </param>
+    <param>
+      <key>decim</key>
+      <value>1</value>
+    </param>
+    <param>
+      <key>_enabled</key>
+      <value>True</value>
+    </param>
+    <param>
+      <key>_coordinate</key>
+      <value>(330, 293)</value>
+    </param>
+    <param>
+      <key>_rotation</key>
+      <value>0</value>
+    </param>
+    <param>
+      <key>id</key>
+      <value>fir_filter_xxx_0</value>
+    </param>
+    <param>
+      <key>maxoutbuf</key>
+      <value>0</value>
+    </param>
+    <param>
+      <key>minoutbuf</key>
+      <value>0</value>
+    </param>
+    <param>
+      <key>samp_delay</key>
+      <value>0</value>
+    </param>
+    <param>
+      <key>taps</key>
+      <value>rrc_taps</value>
+    </param>
+    <param>
+      <key>type</key>
+      <value>ccc</value>
+    </param>
+  </block>
+  <block>
+    <key>fractional_resampler_xx</key>
+    <param>
+      <key>alias</key>
+      <value></value>
+    </param>
+    <param>
+      <key>comment</key>
+      <value></value>
+    </param>
+    <param>
+      <key>affinity</key>
+      <value></value>
+    </param>
+    <param>
+      <key>_enabled</key>
+      <value>1</value>
+    </param>
+    <param>
+      <key>_coordinate</key>
+      <value>(96, 293)</value>
+    </param>
+    <param>
+      <key>_rotation</key>
+      <value>0</value>
+    </param>
+    <param>
+      <key>id</key>
+      <value>fractional_resampler_xx_0</value>
+    </param>
+    <param>
+      <key>maxoutbuf</key>
+      <value>0</value>
+    </param>
+    <param>
+      <key>minoutbuf</key>
+      <value>0</value>
+    </param>
+    <param>
+      <key>phase_shift</key>
+      <value>0</value>
+    </param>
+    <param>
+      <key>resamp_ratio</key>
+      <value>12001./11999.*0+1</value>
+    </param>
+    <param>
+      <key>type</key>
+      <value>complex</value>
+    </param>
+  </block>
   <block>
     <key>import</key>
     <param>
@@ -2315,13 +2464,13 @@
   </connection>
   <connection>
     <source_block_id>blocks_multiply_const_vxx_0</source_block_id>
-    <sink_block_id>blocks_throttle_0</sink_block_id>
+    <sink_block_id>fractional_resampler_xx_0</sink_block_id>
     <source_key>0</source_key>
     <sink_key>0</sink_key>
   </connection>
   <connection>
     <source_block_id>blocks_throttle_0</source_block_id>
-    <sink_block_id>preamble</sink_block_id>
+    <sink_block_id>fir_filter_xxx_0</sink_block_id>
     <source_key>0</source_key>
     <sink_key>0</sink_key>
   </connection>
@@ -2349,6 +2498,18 @@
     <source_key>soft_dec</source_key>
     <sink_key>in</sink_key>
   </connection>
+  <connection>
+    <source_block_id>fir_filter_xxx_0</source_block_id>
+    <sink_block_id>preamble</sink_block_id>
+    <source_key>0</source_key>
+    <sink_key>0</sink_key>
+  </connection>
+  <connection>
+    <source_block_id>fractional_resampler_xx_0</source_block_id>
+    <sink_block_id>blocks_throttle_0</sink_block_id>
+    <source_key>0</source_key>
+    <sink_key>0</sink_key>
+  </connection>
   <connection>
     <source_block_id>preamble</source_block_id>
     <sink_block_id>blocks_complex_to_mag_0</sink_block_id>

lib/adaptive_dfe_impl.cc

@@ -22,9 +22,14 @@
 #include "config.h"
 #endif
 
-#include <gnuradio/io_signature.h>
+#include <boost/format.hpp>
+
 #include <gnuradio/expj.h>
+#include <gnuradio/io_signature.h>
+#include <gnuradio/logger.h>
+
 #include <volk/volk.h>
+
 #include "adaptive_dfe_impl.h"
 
 #define VOLK_SAFE_DELETE(x) \
@@ -96,6 +101,8 @@ adaptive_dfe_impl::adaptive_dfe_impl(int sps, // samples per symbol
   , _hist_symbols(nullptr)
   , _hist_sample_index(0)
   , _hist_symbol_index(0)
+  , _ignore_filter_updates(0)
+  , _saved_samples()
   , _sample_counter(0)
   , _constellations()
   , _npwr()
@@ -118,6 +125,8 @@ adaptive_dfe_impl::adaptive_dfe_impl(int sps, // samples per symbol
   , _ud(0)
   , _state(WAIT_FOR_PREAMBLE)
 {
+  GR_LOG_DECLARE_LOGPTR(d_logger);
+  GR_LOG_ASSIGN_LOGPTR(d_logger, "adaptive_dfe");
   message_port_register_out(_msg_port_name);
 }
 
@@ -152,11 +161,11 @@ adaptive_dfe_impl::general_work(int noutput_items,
 
   int nout = 0;
   int i = 0;
-  for (; i<ninput_items[0] && nout < noutput_items; ++i) {
+  for (; i<ninput_items[0] && nout < noutput_items;) {
     assert(nout < noutput_items);
 
     if (_state == WAIT_FOR_PREAMBLE) {
-      insert_sample(in[i]);
+      insert_sample(in[i++]);
       uint64_t offset = 0;
       float phase_est = 0;
       if (get_correlation_tag(i, offset, phase_est)) {
@@ -168,12 +177,14 @@ adaptive_dfe_impl::general_work(int noutput_items,
         _descrambled_symbols.clear();
         // _hist_sample_index = 0;
         _hist_symbol_index = 0;
+        _ignore_filter_updates = 0;
+        _saved_samples.clear();
         std::fill_n(_hist_symbols, 2*_nW, gr_complex(0));
         std::fill_n(_taps_samples, _nB+_nF+1, gr_complex(0));
         std::fill_n(_taps_symbols, _nW, gr_complex(0));
         _samples.clear();
         _phase = -phase_est;
-        _taps_samples[_nB+1] = 1;
+        _taps_samples[_nB+1] = 0.01;
         _taps_symbols[0] = 1;
         GILLock gil_lock;
         try {
@@ -224,11 +235,11 @@ adaptive_dfe_impl::general_work(int noutput_items,
         }
         continue;
       }
-      _samples.push_back(in[i]);
+      _samples.push_back(in[i++]);
     } // INITIAL_DOPPLER_ESTIMATE_CONTINUE
 
     if (_state == INITIAL_DOPPLER_ESTIMATE_CONTINUE) {
-      std::cout << "INITIAL_DOPPLER_ESTIMATE_CONTINUE\n";
+      GR_LOG_DEBUG(d_logger, "INITIAL_DOPPLER_ESTIMATE_CONTINUE");
       while (!_samples.empty() && nout < noutput_items) {
         insert_sample(_samples.back());
         _sample_counter += 1;
@@ -246,13 +257,15 @@ adaptive_dfe_impl::general_work(int noutput_items,
 
     if (_state == DO_FILTER) {
       if ((_sample_counter%_sps) == 0) {
-        if (_symbol_counter == _symbols.size()) {
+        if (_symbol_counter == _symbols.size()) { // frame is ready
           _symbol_counter = 0;
           GILLock gil_lock;
           try {
+            // update doppler estimate
             update_doppler_information(_physicalLayer.attr("get_doppler")
                                        (complex_vector_to_ndarray(_descrambled_symbols),
                                         complex_vector_to_ndarray(_samples)));
+            // publish soft decisions
             if (!_vec_soft_decisions.empty()) {
               unsigned int const bits_per_symbol = _constellations[_constellation_index]->bits_per_symbol();
               _msg_metadata = pmt::dict_add(_msg_metadata, pmt::mp("bits_per_symbol"), pmt::from_long(bits_per_symbol));
@@ -262,16 +275,31 @@ adaptive_dfe_impl::general_work(int noutput_items,
               _vec_soft_decisions.clear();
             }
             _samples.clear();
+            // get information about the following frame
             update_frame_information(_physicalLayer.attr("get_frame")());
           } catch (boost::python::error_already_set const&) {
             PyErr_Print();
           }
-        }
+        } // frame is ready
+        if (_ignore_filter_updates == 0) {
           out[nout++] = filter();
+          if (_symbol_counter+1 == _symbols.size())
+            recenter_filter_taps();
+        } else {
+          _ignore_filter_updates -= 1;
         }
-      insert_sample(in[i]);
+      } // (_sample_counter%_sps) == 0
-      if (_need_samples)
+
+
+      if (_need_samples) {
         _samples.push_back(_hist_samples[_hist_sample_index+_nB+1]);
+      }
+      if (_saved_samples.empty()) {
+        insert_sample(in[i++]);
+      } else {
+        insert_sample(_saved_samples.back());
+        _saved_samples.pop_back();
+      }
       _sample_counter += 1;
     } // DO_FILTER
   } // next input sample
@@ -306,10 +334,11 @@ bool adaptive_dfe_impl::start()
   std::fill_n(_taps_samples,   (_nB+_nF+1), gr_complex(0));
   std::fill_n(_taps_symbols,           _nW, gr_complex(0));
 
-  _taps_samples[_nB+1] = 1;
+  _taps_samples[_nB+1] = 0.01;
   _taps_symbols[0]     = 1;
 
-  std::cout << "adaptive_dfe_impl::start() " << _nB << " " << _nF << " " << _mu << " " << _alpha << std::endl;
+  GR_LOG_DEBUG(d_logger,str(boost::format("adaptive_dfe_impl::start() nB=%d nF=%d mu=%f alpha=%f")
+                             % _nB % _nF % _mu % _alpha));
   GILLock gil_lock;
   try {
     boost::python::object module        = boost::python::import(boost::python::str("digitalhf.physical_layer." + _py_module_name));
@@ -325,7 +354,7 @@ bool adaptive_dfe_impl::start()
 bool adaptive_dfe_impl::stop()
 {
   gr::thread::scoped_lock lock(d_setlock);
-  std::cout << "adaptive_dfe_impl::stop()" << std::endl;
+  GR_LOG_DEBUG(d_logger, "adaptive_dfe_impl::stop()");
   GILLock gil_lock;
   _physicalLayer = boost::python::object();
   VOLK_SAFE_DELETE(_taps_samples);
@@ -334,6 +363,7 @@ bool adaptive_dfe_impl::stop()
   VOLK_SAFE_DELETE(_hist_symbols);
   return true;
 }
+
 gr_complex adaptive_dfe_impl::filter() {
   gr_complex filter_output = 0;
   volk_32fc_x2_dot_prod_32fc(&filter_output,
@@ -371,17 +401,10 @@ gr_complex adaptive_dfe_impl::filter() {
     known_symbol = _scramble[_symbol_counter] * descrambled_symbol;
   }
   gr_complex err =  filter_output - known_symbol;
-  int jMax=0;
-  float tMax=0;
   for (int j=0; j<_nB+_nF+1; ++j) {
     assert(_hist_sample_index+j < 2*(_nB+_nF+1));
     _taps_samples[j] -= _mu*err*std::conj(_hist_samples[_hist_sample_index+j]);
-    // if (std::abs(_taps_samples[j]) > tMax) {
-    //   tMax = std::abs(_taps_samples[j]);
-    //   jMax = j;
-    // }
   }
-  // std::cout << "taps_max: " << jMax << " " << tMax << std::endl;
   for (int j=0; j<_nW; ++j) {
     assert(_hist_symbol_index+j < 2*_nW);
     _taps_symbols[j] -= _mu*err*std::conj(_hist_symbols[_hist_symbol_index+j]) + _alpha*_taps_symbols[j];
@@ -397,9 +420,49 @@ gr_complex adaptive_dfe_impl::filter() {
   return filter_output*std::conj(_scramble[_symbol_counter++]);
 }
 
+void adaptive_dfe_impl::recenter_filter_taps() {
+  // get max(abs(taps))
+  ssize_t const idx_max = std::distance(_taps_samples,
+                                        std::max_element(_taps_samples+_nB+1-3*_sps, _taps_samples+_nB+1+3*_sps,
+                                                         [](gr_complex a, gr_complex b) {
+                                                           return std::norm(a) < std::norm(b);
+                                                         }));
+  GR_LOG_DEBUG(d_logger, str(boost::format("idx_max=%2d abs(tap_max)=%f") % idx_max % std::abs(_taps_samples[idx_max])));
+
+  if (idx_max-_nB-1 >= 2*_sps && _saved_samples.empty() && _ignore_filter_updates==0) {
+    // maximum is right of the center tap
+    //   -> shift taps to the left left
+    GR_LOG_DEBUG(d_logger, "shift left");
+    std::copy(_taps_samples+2*_sps, _taps_samples+_nB+_nF+1, _taps_samples);
+    std::fill_n(_taps_samples+_nB+_nF+1-2*_sps, 2*_sps, gr_complex(0));
+    //      and omit the next two calls to filter in order to keep the alignment between samples and taps
+    _ignore_filter_updates = 2;
+
+  } else if (idx_max-_nB-1 <= -2*_sps && _saved_samples.empty() && _ignore_filter_updates==0) {
+    // maximum is left of the center tap
+    //   -> shift taps to the right
+    GR_LOG_DEBUG(d_logger, "shift right");
+    std::copy_backward(_taps_samples, _taps_samples+_nB+_nF+1-2*_sps,
+                       _taps_samples+_nB+_nF+1);
+    std::fill_n(_taps_samples, 2*_sps, gr_complex(0));
+    //      save the last 2*_sps samples (will be reinserted)
+    _saved_samples.resize(2*_sps);
+    std::reverse_copy(_hist_samples+_hist_sample_index+(_nB+_nF+1)-2*_sps,
+                      _hist_samples+_hist_sample_index+(_nB+_nF+1),
+                      _saved_samples.begin());
+    //      shift samples index
+    _hist_sample_index += (_nB+_nF+1)-2*_sps;
+    _hist_sample_index %= (_nB+_nF+1);
+    //      set the 1st 2*_sps unknown old samples to zero
+    for (int l=_hist_sample_index; l<_hist_sample_index+2*_sps; ++l) {
+      int const k = (l+_nB+_nF+1)%(2*(_nB+_nF+1));
+      _hist_samples[l] = _hist_samples[k] = gr_complex(0);
+    }
+  }
+}
 void adaptive_dfe_impl::set_mode(std::string mode) {
   gr::thread::scoped_lock lock(d_setlock);
-  std::cout << "adaptive_dfe_impl::set_mode " << mode << std::endl;
+  GR_LOG_DEBUG(d_logger, "adaptive_dfe_impl::set_mode "+ mode);
   GILLock gil_lock;
   try {
     _physicalLayer.attr("set_mode")(mode);
@@ -483,13 +546,14 @@ void adaptive_dfe_impl::update_pll(float doppler) {
     _ud  = delta_f;
     _df +=_b[0]*_ud + _b[1]*ud_old;
   }
-  std::cout << "PLL: " << _df << " " << delta_f << std::endl;
+  GR_LOG_DEBUG(d_logger, str(boost::format("PLL: df=%f delta_f=%f (rad/symb)") % _df % delta_f));
 }
 void adaptive_dfe_impl::insert_sample(gr_complex z) {
   // insert sample into the circular buffer
   _hist_samples[_hist_sample_index] = _hist_samples[_hist_sample_index+_nB+_nF+1] = z * gr_expj(-_phase);
   if (++_hist_sample_index == _nB+_nF+1)
     _hist_sample_index = 0;
+  if (z != gr_complex(0))
     update_local_oscillator();
 }
 void adaptive_dfe_impl:: update_local_oscillator() {
@@ -503,7 +567,7 @@ bool adaptive_dfe_impl::get_correlation_tag(uint64_t i, uint64_t& offset, float&
   std::vector<tag_t> v;
   get_tags_in_window(v, 0, i,i+1);
   for (int j=0; j<v.size(); ++j) {
-    std::cout << "tag " << v[j].key << " " << v[j].offset-nitems_read(0) << std::endl;
+    std::cout << "tag " << v[j].key << " " << v[j].offset-nitems_read(0) << " " << v[j].value << std::endl;
     if (v[j].key == pmt::mp("phase_est")) {
       phase_est = pmt::to_double(v[j].value);
       std::cout << "phase_est " << v[j].offset <<" " << nitems_read(0) << " " << phase_est << std::endl;

lib/adaptive_dfe_impl.h

@@ -50,7 +50,10 @@ private:
   int _hist_sample_index;
   int _hist_symbol_index;
 
-  std::size_t _sample_counter;
+  int _ignore_filter_updates;
+  std::vector<gr_complex> _saved_samples;
+
+  uint64_t _sample_counter;
 
   std::vector<gr::digital::constellation_sptr> _constellations;
   std::vector<float> _npwr;
@@ -88,6 +91,7 @@ private:
 
   void update_local_oscillator();
   gr_complex filter();
+  void recenter_filter_taps();
 
   void insert_sample(gr_complex z);
   void update_pll(float doppler);

python/physical_layer/STANAG_4285.py

@@ -68,12 +68,14 @@ class PhysicalLayer(object):
         if len(symbols) != 0:
             idx = range(30,80) if self._is_first_frame else range(80)
             z = symbols[idx]*np.conj(self._preamble['symb'][idx])
+            print('quality_preamble',np.sum(np.real(z)<0))
             success = np.sum(np.real(z)<0) < 30
         return success,doppler
 
     def quality_data(self, s):
         """quality check for the data frame"""
         known_symbols = np.mod(range(176),48)>=32
+        print('quality_data',np.sum(np.real(s[known_symbols])<0))
         success = np.sum(np.real(s[known_symbols])<0) < 20
         return success,0 ## no doppler estimate for data frames