wavelettransform.hh

#ifndef __WAVELETTRANSFORM_HH__
#define __WAVELETTRANSFORM_HH__

#include "waveletanalysis.hh"
#include "convolution.hh"
#include "movingaverage.hh"
#include <vector>


namespace wt {

class WaveletTransform: public WaveletAnalysis
{
public:
  WaveletTransform(const Wavelet &wavelet, const RVector &scales, bool subSample=false);
  WaveletTransform(const Wavelet &wavelet, double *scales, int Nscales, bool subSample=false);
  WaveletTransform(const WaveletAnalysis &other, bool subSample=false);
  virtual ~WaveletTransform();

  template <class iDerived, class oDerived>
  void operator() (const Eigen::DenseBase<iDerived> &signal, Eigen::DenseBase<oDerived> &out)
  {
    // signal length
    int N = signal.size();
    // Get start indices for each transformation block
    std::vector<size_t> blockIdxs(_filterBank.size());
    for (size_t j=0; j<_filterBank.size(); j++) {
      if (0 == j) {
        blockIdxs[0] = 0;
      } else {
        blockIdxs[j] = blockIdxs[j-1] + _filterBank[j-1]->numKernels();
      }
    }

    // Iterate over all convolution filters grouping wavelets with the same size
    for (size_t j=0; j<_filterBank.size(); j++)
    {
      // Get convolution filters
      Convolution *filters = _filterBank[j];
      // Get subsampling
      int M = filters->subSampling();
      // Get start column in output matrix
      size_t outCol = blockIdxs[j];
      // Get number of kernels in group / # columns in out
      int K = filters->numKernels();

      if (1 == M) {
        // w/o sub-sampling -> direct overlap-add convolution
        filters->apply(signal, out.block(0, outCol, N, K).derived());
        // continue with next block
        continue;
      }

      /*
       * Perform convolution with sub-sampling
       */
      // Number of samples in the sub-sampled signal
      int n = WT_IDIV_CEIL(N,M);
      // allocate some working arrays
      CVector subsig(n);     // <- Will hold the sub-sampled signal
      CMatrix subres(n, K);  // <- Will hold the sub-sampled result
      // For every shift within the M-fold sub-sampling:
      for (int m=0; m<M; m++) {
        // subsample input signal
        for (int i=0; i<n; i++) {
          subsig[i] = ((i*M+m) < N) ? CScalar(signal[i*M+m]) : 0;
        }
        // Apply overlap-add convolution
        filters->apply(subsig, subres.derived());
        // Store results into output buffer (interleaved rows)
        for (int i=0; i<n; i++) {
          if ( (i*M+m) < N ) {
            out.block(i*M+m, outCol, 1, K) = subres.row(i);
          }
        }
      }
    }
  }

protected:
  void init_trafo();

protected:
  bool _subSample;
  std::vector<Convolution *> _filterBank;
};

}

#endif // __WAVELETTRANSFORM_HH__