IterDecon.cpp

Go to the documentation of this file.

#include "IterDecon.h"
#include "miscfunc.h"
#include <boost/algorithm/minmax_element.hpp>
#include <iostream>
#include <cmath>

namespace gplib
  {

    IterDecon::IterDecon(const int inputsize, TsSpectrum &Spec) :
      AdaptiveFilter(inputsize, inputsize), Weights(inputsize), Spectrum(Spec)
      {
        for (int i = 0; i < inputsize; ++i)
          Weights(i) = 0;
      }

    IterDecon::~IterDecon()
      {
      }

    void IterDecon::PrintWeights(std::ostream &output)
      {
        std::copy(Weights.begin(), Weights.end(),
            std::ostream_iterator<double>(output, "\n"));
      }

    void IterDecon::AdaptFilter(const gplib::rvec &Input,
        const gplib::rvec &Desired)
      {
        typedef gplib::rvec::const_iterator iterator;
        gplib::rvec Corr(Desired.size());
        Correl(Desired, Input, Corr, Spectrum); //calculate the correlation between desired and input
        const double power = ublas::prec_inner_prod(Input, Input);
        std::pair<iterator, iterator> MinMaxPos(boost::minmax_element(
            Corr.begin(), Corr.begin() + Corr.size() / 2)); //find minimum and maximum
        int MaxIndex = MinMaxPos.second - Corr.begin(); //index of maximum
        int MinIndex = MinMaxPos.first - Corr.begin(); //index of minimum

        if (std::abs(*MinMaxPos.second) > std::abs(*MinMaxPos.first)) //if max has higher absolute value than min
          Weights(MaxIndex) += *MinMaxPos.second / power; //update weight for max
        else
          Weights(MinIndex) += *MinMaxPos.first / power; //update weight for min
      }

    void IterDecon::CalcOutput(const gplib::rvec &Input, gplib::rvec &Output)
      {
        Convolve(Input, Weights, Output, Spectrum);//multiplication with a toeplitz matrix is the same as cyclic convolution
        SetOutput(Output);
      }
  }