attempt to regularize kalman_hsu

lib/adaptive_dfe_impl.cc

@@ -232,10 +232,10 @@ bool adaptive_dfe_impl::start()
   GR_LOG_DEBUG(d_logger,str(boost::format("adaptive_dfe_impl::start() nB=%d nF=%d mu=%f alpha=%f")
                              % _nB % _nF % _mu % _alpha));
   _filter_update = lms::make(_mu);
-  // _filter_update = nlms::make(0.5);
+  //_filter_update = nlms::make(0.5f, 1.0f);
-  // _filter_update = rls::make(0.001f, 0.999f); // not stable
+  //_filter_update = rls::make(0.001f, 0.9f); // not stable
-  // _filter_update = kalman_exp::make(1.0f, 0.999f); // too slow
+  //_filter_update = kalman_exp::make(1.0f, 0.9f); // too slow
-  // _filter_update = kalman_hsu::make(0.008f, 0.1f); // not stable
+  // _filter_update = kalman_hsu::make(0.02f, 0.1f); // not stable
   return true;
 }
 bool adaptive_dfe_impl::stop()

lib/kalman_exp.hpp

@@ -16,7 +16,6 @@ namespace digitalhf {
 
 class kalman_exp : public filter_update {
 public:
-  kalman_exp(float r, float lambda);
   virtual ~kalman_exp();
 
   static sptr make(float r, float lambda);
@@ -29,6 +28,8 @@ protected:
   void resize(size_t);
 
 private:
+  kalman_exp(float r, float lambda);
+
   typedef std::vector<gr_complex, volk_allocator<gr_complex> > vec_type;
   float    _lambda;
   float    _r;

lib/kalman_hsu.cc

@@ -3,7 +3,7 @@
 #include <cassert>
 #include "kalman_hsu.hpp"
 #include <volk/volk.h>
-
+#include <iostream>
 namespace gr {
 namespace digitalhf {
 
@@ -92,8 +92,10 @@ gr_complex const* kalman_hsu::update(gr_complex const* beg,
     }
   }
 
-  for (unsigned i=0; i<n; ++i)
+  for (unsigned i=0; i<n; ++i) {
+    _d[i] = 1.0f/(10.0f+1.0f/_d[i]);// regularization
     _g[i] *= y;
+  }
 
   return &_g[0];
 }

lib/kalman_hsu.hpp

@@ -20,7 +20,6 @@ namespace digitalhf {
 
 class kalman_hsu : public filter_update {
 public:
-  kalman_hsu(float q, float e);
   virtual ~kalman_hsu();
 
   static sptr make(float q, float e);
@@ -39,6 +38,8 @@ protected:
     return i+j*(j-1)/2; // lower-triangular matrix index -> linear index
   }
 private:
+  kalman_hsu(float q, float e);
+
   typedef std::vector<gr_complex, volk_allocator<gr_complex > > complex_vec_type;
   typedef std::vector<float, volk_allocator<float> > real_vec_type;
   float _q;

lib/lms.hpp

@@ -14,7 +14,6 @@ namespace digitalhf {
 
 class lms : public filter_update {
 public:
-  lms(float mu);
   virtual ~lms();
 
   static sptr make(float mu);
@@ -27,6 +26,8 @@ protected:
   void resize(size_t);
 
 private:
+  lms(float mu);
+
   typedef std::vector<gr_complex, volk_allocator<gr_complex> > vec_type;
   float    _mu;
   vec_type _gain;

lib/nlms.cc

@@ -7,12 +7,13 @@
 namespace gr {
 namespace digitalhf {
 
-filter_update::sptr nlms::make(float mu) {
+filter_update::sptr nlms::make(float mu, float delta) {
-  return filter_update::sptr(new nlms(mu));
+  return filter_update::sptr(new nlms(mu, delta));
 }
 
-nlms::nlms(float mu)
+nlms::nlms(float mu, float delta)
   : _mu(mu)
+  , _delta(delta)
   , _gain()
   , _tmp()
   , _t1() {
@@ -35,21 +36,22 @@ void nlms::reset() {
 }
 
 gr_complex const* nlms::update(gr_complex const* beg,
-                              gr_complex const* end) {
+                               gr_complex const* end) {
   assert(end-beg > 0);
   size_t const n = end - beg;
   resize(n);
   volk_32fc_conjugate_32fc(&_tmp[0], beg, n);
   volk_32fc_magnitude_squared_32f(&_t1[0], beg, n);
-  float norm = 0.0f;
+  float norm = _delta;
   volk_32f_accumulator_s32f(&norm, &_t1[0], n);
-  volk_32f_s32f_multiply_32f((float*)&_gain[0], (float const*)&_tmp[0], _mu/norm, 2*n);
+  volk_32f_s32f_multiply_32f((float*)&_gain[0], (float const*)&_tmp[0], std::max(0.006f, _mu/norm), 2*n);
 
   return &_gain.front();
 }
 
 void nlms::set_parameters(std::map<std::string, float>const & p) {
-  _mu = p.at("mu");
+  _mu    = p.at("mu");
+  _delta = p.at("delta");
 }
 
 } // namespace digitalhf

lib/nlms.hpp

@@ -14,10 +14,9 @@ namespace digitalhf {
 
 class nlms : public filter_update {
 public:
-  nlms(float mu);
   virtual ~nlms();
 
-  static sptr make(float mu);
+  static sptr make(float mu, float delta);
 
   virtual void reset();
   virtual gr_complex const* update(gr_complex const*, gr_complex const*);
@@ -27,9 +26,12 @@ protected:
   void resize(size_t);
 
 private:
+  nlms(float mu, float delta);
+
   typedef std::vector<gr_complex, volk_allocator<gr_complex> > complex_vec_type;
   typedef std::vector<float, volk_allocator<float> > real_vec_type;
   float _mu;
+  float _delta;
   complex_vec_type _gain;
   complex_vec_type _tmp;
   real_vec_type    _t1;

lib/rls.cc

@@ -53,7 +53,7 @@ gr_complex const* rls::update(gr_complex const* beg,
   volk_32fc_x2_conjugate_dot_prod_32fc(&uPu, &_pu[0], beg, n);
 
   volk_32fc_conjugate_32fc(&_tmp[0], &_pu[0], n);
-  volk_32f_s32f_multiply_32f((float*)&_gain[0], (float const*)&_tmp[0], 1.0f/(_lambda + std::real(uPu)), 2*n);
+  volk_32f_s32f_multiply_32f((float*)&_gain[0], (float const*)&_tmp[0], 1.0f/(_lambda + std::real(uPu)/_lambda), 2*n);
 //  volk_32fc_s32fc_multiply_32fc(&_gain[0], &_tmp[0], 1/(_lambda + uPu), n);
 
   for (unsigned i=0; i<n; ++i) {
@@ -61,7 +61,6 @@ gr_complex const* rls::update(gr_complex const* beg,
     volk_32fc_s32fc_multiply_32fc(&_tmp[0], &_gain[0], -_pu[i], n);
     volk_32fc_x2_add_32fc(&_inv_corr[k], &_inv_corr[k], &_tmp[0], n);
     volk_32f_s32f_multiply_32f((float*)&_inv_corr[k], (float const*)&_inv_corr[k], 1.0f/_lambda, 2*n);
-    _inv_corr[k+i] = std::real(_inv_corr[k+i]);
   }
 
   return &_gain.front();

lib/rls.hpp

@@ -14,7 +14,6 @@ namespace digitalhf {
 
 class rls : public filter_update {
 public:
-  rls(float delta, float lambda);
   virtual ~rls();
 
   static sptr make(float delta, float lambda);
@@ -27,6 +26,8 @@ protected:
   void resize(size_t);
 
 private:
+  rls(float delta, float lambda);
+
   typedef std::vector<gr_complex, volk_allocator<gr_complex> > vec_type;
   float    _delta;
   float    _lambda;