/*!
 * \file
 * \brief Definition of frequency domain filter class
 * \author Simon Wood
 *
 * -------------------------------------------------------------------------
 *
 * IT++ - C++ library of mathematical, signal processing, speech processing,
 *        and communications classes and functions
 *
 * Copyright (C) 1995-2008  (see AUTHORS file for a list of contributors)
 *
 * This program 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.
 *
 * This program 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 should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 *
 * -------------------------------------------------------------------------
 */

#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>


namespace itpp {

  /*!
    \brief Freq_Filt Frequency domain filtering using the overlap-add technique
    \ingroup filters

    The Freq_Filt class implements an FFT based filter using the overlap-add
    technique. The data is filtered by first transforming the input sequence
    into the frequency domain with an efficient FFT implementation (i.e. FFTW)
    and then multiplied with a Fourier transformed version of the impulse
    response. The resulting data is then inversed Fourier transformed to return
    a filtered time domain signal.

    Freq_Filt is a templated class. The template paramter \c Num_T defines the
    data type for the impulse response \c b, input data \c x and output data
    \c y.

    The class constructor chooses an optimal FFT length and data block
    size that minimizes execution time.

    For example,
    \code
    vec b = "1 2 3 4";
    Freq_Filt<double> FF(b,8000);
    \endcode

    where 8000 corresponds to the expected vector length of the data
    to be filtered. The actual number of elements does not have to be
    exact, but it should be close.

    Here is a complete example:
    \code
    vec b = "1 2 3 4";
    vec x(20);
    x(0) = 1;

    // Define a filter object for doubles
    Freq_Filt<double> FF(b,x.length());

    // Filter the data
    vec y = FF.filter(x);

    // Check the FFT and block sizes that were used
    int fftsize = FF.getFFTSize();
    int blk = FF.getBlkSize();
    \endcode

    To facilitate large data sets the Freq_Filt class has a streaming
    option. In this mode of operation data history is internally
    stored. This allows the class to be used for large data sets that
    are read from disk or streamed in real-time.

    \code
    bool stream = true;
    vec b = "1 2 3 4";
    Freq_Filt<double> FF(b,1000);

    vec x,y;
    while(!EOF) {
    x = "read buffer of data";
    y = FF.filter(x,stream);
    cout << << endl;
    }
    \endcode
  */


  template<class Num_T>
    class Freq_Filt {
    public:
    /*!
      \brief Constructor

      The empty constructor makes it possible to have other container objects
      of the Freq_Filt class
    */
    Freq_Filt() {}

    /*!
      \brief Constructor with initialization of the impulse response \b b.

      Create a filter object with impulse response \b b. The FFT size and
      data block size are also initialized.
      \code
      vec b = "1 2 3 4";
      vec x(20);
      Freq_Filt FF(b,x.length());
      \endcode
    */
    Freq_Filt(const Vec<Num_T> &b, const int xlength) {init(b,xlength);}

    /*!
      \brief Filter data in the input vector \b x

      Filters data in batch mode or streaming mode
      \code
      FF.filter(x);	// Filters data in batch mode
      FF.filter(x,1);	// Filters data in streaming mode
      \endcode
    */
    Vec<Num_T> filter(const Vec<Num_T> &x, const int strm = 0);

    //! Return FFT size
    int get_fft_size() { return fftsize; }

    //! Return the data block size
    int get_blk_size() { return blksize; }

    //! Destructor
    ~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 aare 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;
    }

} // namespace itpp

#endif // #ifndef FREQ_FILT_H