FIR.H

Go to the documentation of this file.

/* Copyright 2000-2018 Matt Flax <flatmax@flatmax.org>
   This file is part of GTK+ IOStream class set

   GTK+ IOStream 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 2 of the License, or
   (at your option) any later version.

   GTK+ IOStream 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 have received a copy of the GNU General Public License
   along with GTK+ IOStream
*/

#ifndef FIR_H
#define FIR_H

// Debug
#include "Debug.H"
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wignored-attributes"
#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
#include <Eigen/Dense>
#include <unsupported/Eigen/FFT>
#pragma GCC diagnostic pop

#define FIR_BLOCKSIZE_MISMATCH_ERROR FIR_ERROR_OFFSET-1
#define FIR_H_EMPTY_ERROR FIR_ERROR_OFFSET-2
#define FIR_CHANNEL_MISMATCH_ERROR FIR_ERROR_OFFSET-3

class FIRDebug :  virtual public Debug  {
public:
    FIRDebug(){
#ifndef NDEBUG
errors[FIR_BLOCKSIZE_MISMATCH_ERROR]=std::string("The input data was not of the same length you used as the variable for the method init. ");
errors[FIR_H_EMPTY_ERROR]=std::string("The fileter h is empty, please load using loadTimeDomainCoefficients. ");
errors[FIR_CHANNEL_MISMATCH_ERROR]=std::string("The input, output and h columns (channels) are not the same count. ");

#endif // NDEBUG
    }
};

#ifndef FP_TYPE
typedef double FP_TYPE; 
#endif

class FIR {
  Eigen::FFT<FP_TYPE> fft; 
  Eigen::Array<Eigen::FFT<FP_TYPE>::Complex, Eigen::Dynamic, Eigen::Dynamic> H; 
  Eigen::Matrix<FP_TYPE, Eigen::Dynamic, Eigen::Dynamic> h; 
  Eigen::Matrix<FP_TYPE, Eigen::Dynamic, Eigen::Dynamic> x; 
  Eigen::Matrix<FP_TYPE, Eigen::Dynamic, 1> yTemp; 
  Eigen::Array<Eigen::FFT<FP_TYPE>::Complex, Eigen::Dynamic, 1> Y; 

  void resetDFT();
protected:
  unsigned int N; 
  Eigen::Matrix<FP_TYPE, Eigen::Dynamic, Eigen::Dynamic> y; 
public:
    FIR(){N=0;} 

    void init(unsigned int blockSize);

    int loadTimeDomainCoefficients(const std::string fileName);

    void loadTimeDomainCoefficients(const Eigen::Matrix<FP_TYPE, Eigen::Dynamic, Eigen::Dynamic> hIn);

    template<typename Derived, typename DerivedOther>
    void filter(const Eigen::MatrixBase<Derived> &input, Eigen::DenseBase<DerivedOther> const &output) {
      if (input.rows()!=N){
        FIRDebug().evaluateError(FIR_BLOCKSIZE_MISMATCH_ERROR);
        return;
      }
      if (input.cols()!=h.cols() || output.cols() != h.cols()){
        FIRDebug().evaluateError(FIR_CHANNEL_MISMATCH_ERROR);
        return;
      }
      if (h.rows()==0) {
        FIRDebug().evaluateError(FIR_H_EMPTY_ERROR);
        return;
      }

      if (x.cols() != input.cols()){ // resize if necessary
        x.setZero(h.rows(), input.cols());
        y.setZero(h.rows(), input.cols());
      }
      y.topRows(y.rows()-N)=y.bottomRows(y.rows()-N); // keep the residual
      y.bottomRows(N).setZero();

      x.topRows(N)=input;

      for (int i=0; i<x.cols(); i++){ // perform the filter on each column
        fft.fwd(Y.data(), x.col(i).data(), x.rows()); // find the DFT of X=Z=dft(x) (store in Y)
        Y*=H.col(i); // convolve X (which is Y) with H
        fft.inv(yTemp.data(), Y.data(), Y.rows()); // take back to the time domain
        y.col(i)+=yTemp; // add to the residual
      }
      const_cast< Eigen::DenseBase<DerivedOther>& >(output)=y.topRows(N);
    }
};
#endif // FIR_H