WienerInterpolator.cpp

Go to the documentation of this file.

#include "WienerInterpolator.h"
#include <iostream>
#include <iterator>
#include <algorithm>

namespace gplib
  {

    WienerInterpolator::WienerInterpolator(const int filterlength) :
      LSSOFilter(filterlength), xms(0.0)
      {
      }

    WienerInterpolator::~WienerInterpolator()
      {
      }

    void WienerInterpolator::AdaptFilter(const gplib::rvec &Input,
        const gplib::rvec &Desired)
      {
        const int size = Input.size();
        const int filterlength = GetWeights().size();
        gplib::rvec wk1(size);
        gplib::rvec wk2(size);
        gplib::rvec wkm(filterlength);
        double power = inner_prod(Input, Input);
        xms = power / size;
        std::fill_n(wkm.begin(), filterlength, 0.0);

        wk1(0) = Input(0);
        wk2(size - 2) = Input(size - 1);
        for (int i = 1; i < size - 1; ++i)
          {
            wk1(i) = Input(i);
            wk2(i - 1) = Input(i);
          }
        for (int j = 0; j < filterlength; ++j)
          {
            double num = 0.0;
            double denom = 0.0;
            for (int k = 0; k < (size - j - 1); ++k)
              {
                num += wk1(k) * wk2(k);
                denom += wk1(k) * wk1(k) + wk2(k) * wk2(k);
              }
            SetWeights()(j) = 2.0 * num / denom;
            xms *= (1.0 - GetWeights()(j) * GetWeights()(j));
            for (int k = 0; k < j; ++k)
              {
                SetWeights()(k) = wkm(k) - GetWeights()(j) * wkm(j - k - 1);
              }
            if (j == filterlength - 1)
              {
                reverse(SetWeights().begin(), SetWeights().end());
                return;
              }
            for (int k = 0; k <= j; ++k)
              wkm(k) = GetWeights()(k);
            for (int k = 0; k < size - j - 2; ++k)
              {
                wk1(k) -= wkm(j) * wk2(k);
                wk2(k) = wk2(k + 1) - wkm(j) * wk1(k + 1);
              }
          }
      }
  }