freq_filt.h

Go to the documentation of this file.

#ifndef FREQ_FILT_H
#define FREQ_FILT_H

#include <itpp/base/vec.h>
#include <itpp/base/converters.h>
#include <itpp/base/math/elem_math.h>
#include <itpp/base/matfunc.h>
#include <itpp/base/specmat.h>
#include <itpp/base/math/min_max.h>
#include <itpp/signal/transforms.h>
#include <itpp/base/math/elem_math.h>
#include <itpp/itexports.h>


namespace itpp
{

template<class Num_T>
class Freq_Filt
{
public:
  Freq_Filt() {}

  Freq_Filt(const Vec<Num_T> &b, const int xlength) {init(b, xlength);}

  Vec<Num_T> filter(const Vec<Num_T> &x, const int strm = 0);

  int get_fft_size() { return fftsize; }

  int get_blk_size() { return blksize; }

  ~Freq_Filt() {}

private:
  int fftsize, blksize;
  cvec B; // FFT of impulse vector
  Vec<Num_T> impulse;
  Vec<Num_T> old_data;
  cvec zfinal;

  void init(const Vec<Num_T> &b, const int xlength);
  vec overlap_add(const vec &x);
  svec overlap_add(const svec &x);
  ivec overlap_add(const ivec &x);
  cvec overlap_add(const cvec &x);
  void overlap_add(const cvec &x, cvec &y);
};


// Initialize impulse rsponse, FFT size and block size
template <class Num_T>
void Freq_Filt<Num_T>::init(const Vec<Num_T> &b, const int xlength)
{
  // Set the impulse response
  impulse = b;

  // Get the length of the impulse response
  int num_elements = impulse.length();

  // Initizlize old data
  old_data.set_size(0, false);

  // Initialize the final state
  zfinal.set_size(num_elements - 1, false);
  zfinal.zeros();

  vec Lvec;

  /*
   * Compute the FFT size and the data block size to use.
   * The FFT size is N = L + M -1, where L corresponds to
   * the block size (to be determined) and M corresponds
   * to the impulse length. The method used here is designed
   * to minimize the (number of blocks) * (number of flops per FFT)
   * Use the FFTW flop computation of 5*N*log2(N).
   */
  vec n = linspace(1, 20, 20);
  n = pow(2.0, n);
  ivec fftflops = to_ivec(elem_mult(5.0 * n, log2(n)));

  // Find where the FFT sizes are > (num_elements-1)
  //ivec index = find(n,">", static_cast<double>(num_elements-1));
  ivec index(n.length());
  int cnt = 0;
  for (int ii = 0; ii < n.length(); ii++) {
    if (n(ii) > (num_elements - 1)) {
      index(cnt) = ii;
      cnt += 1;
    }
  }
  index.set_size(cnt, true);

  fftflops = fftflops(index);
  n = n(index);

  // Minimize number of blocks * number of flops per FFT
  Lvec = n - (double)(num_elements - 1);
  int min_ind = min_index(elem_mult(ceil((double)xlength / Lvec), to_vec(fftflops)));
  fftsize = static_cast<int>(n(min_ind));
  blksize = static_cast<int>(Lvec(min_ind));

  // Finally, compute the FFT of the impulse response
  B = fft(to_cvec(impulse), fftsize);
}

// Filter data
template <class Num_T>
Vec<Num_T> Freq_Filt<Num_T>::filter(const Vec<Num_T> &input, const int strm)
{
  Vec<Num_T> x, tempv;
  Vec<Num_T> output;

  /*
   * If we are not in streaming mode we want to process all of the data.
   * However, if we are in streaming mode we want to process the data in
   * blocks that are commensurate with the designed 'blksize' parameter.
   * So, we do a little book keeping to divide the iinput vector into the
   * correct size.
   */
  if (!strm) {
    x = input;
    zfinal.zeros();
    old_data.set_size(0, false);
  }
  else { // we are in streaming mode
    tempv = concat(old_data, input);
    if (tempv.length() <= blksize) {
      x = tempv;
      old_data.set_size(0, false);
    }
    else {
      int end = tempv.length();
      int numblks = end / blksize;
      if ((end % blksize)) {
        x = tempv(0, blksize * numblks - 1);
        old_data = tempv(blksize * numblks, end - 1);
      }
      else {
        x = tempv(0, blksize * numblks - 1);
        old_data.set_size(0, false);
      }
    }
  }
  output = overlap_add(x);

  return output;
}

// Overlap-add routine
template<class Num_T>
void Freq_Filt<Num_T>::overlap_add(const cvec&x, cvec &y)
{
  int nb = impulse.length();
  int nx = x.length();

  y.set_size(nx, false);
  y.zeros();
  cvec X, Y;
  int istart = 0;
  int L = blksize;
  while (istart < nx) {
    int iend = std::min(istart + L - 1, nx - 1);

    X = fft(x(istart, iend), fftsize);
    Y = ifft(elem_mult(X, B));
    Y.set_subvector(0, Y(0, nb - 2) + zfinal);
    int yend = std::min(nx - 1, istart + fftsize - 1);
    y.set_subvector(istart, Y(0, yend - istart));
    zfinal = Y(fftsize - (nb - 1), fftsize - 1);
    istart += L;
  }
}

template<class Num_T>
vec Freq_Filt<Num_T>::overlap_add(const vec &x)
{
  cvec y; // Size of y is set later
  overlap_add(to_cvec(x), y);
  return real(y);
}

template<class Num_T>
svec Freq_Filt<Num_T>::overlap_add(const svec &x)
{
  cvec y; // Size of y is set later
  overlap_add(to_cvec(x), y);
  return to_svec(real(y));
}

template<class Num_T>
ivec Freq_Filt<Num_T>::overlap_add(const ivec &x)
{
  cvec y; // Size of y is set later
  overlap_add(to_cvec(x), y);
  return to_ivec(real(y));
}

template<class Num_T>
cvec Freq_Filt<Num_T>::overlap_add(const cvec &x)
{
  cvec y; // Size of y is set later
  overlap_add(x, y);
  return y;
}


// ----------------------------------------------------------------------
// Instantiations
// ----------------------------------------------------------------------

ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Freq_Filt<double>;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Freq_Filt<std::complex<double> >;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Freq_Filt<short>;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Freq_Filt<int>;


} // namespace itpp

#endif // #ifndef FREQ_FILT_H