#ifndef _WAVETHRESH_H_
#define _WAVETHRESH_H_


#include <math.h>
#include <stdio.h>

/**

  Support for wavelet coefficient thresholding.

  In wavelet coefficient thresholding a selected set of
  the wavelet coefficients that result from the wavelet
  transform are modified.

  The algorithm used here uses the "Universal Threshold" that Donoho
  and Johnstone have discussed in their papers.  Thresholds and the
  "Universal Threshold" are also discussed in <i>Noise Reduction by
  Wavelet Thresholding</i> by Jansen, Springer Verlag, 2001.

  The "Universal Threshold is:

  <pre>
     UnivThresh = sqrt( 2 * log( N ) ) * stdDev
  </pre>

  Here stdDev is the standard deviation of N wavelet coefficients.
  Log is the natural log.

  Thresholding is used to remove noise.  Another way to state this
  is that the signal is smoothed.  This code uses "soft thresholding"
  where a wavelet coefficient that is less than the threshold is set
  to zero and a coefficient whose absolute value is greater than the
  threshold is moved toward zero by the threshold amount.

  This class provides globally accessible static functions.  It
  is not meant to be declared as an object and has no state.

 */
template<class T>
class wavethresh
{
private:
  wavethresh( const wavethresh &rhs );
  wavethresh() {}
  ~wavethresh() {}

public:

  /**
    A container for the mean and standard deviation
   */
  class bellInfo {
  private:
    double mean_;
    double stdDev_;
  public:
    bellInfo()
    {
      mean_ = 0.0;
      stdDev_ = 0.0;
    }
    ~bellInfo() {}
    bellInfo( const bellInfo &rhs )
    {
      mean_ = rhs.mean_;
      stdDev_ = rhs.stdDev_;
    }

    void mean( const double m ) { mean_ = m; }
    double mean() { return mean_; }

    void stdDev( const double s ) { stdDev_ = s; }
    double stdDev() { return stdDev_; }
  }; // bellInfo


  /**
    Calculate teh mean and standard deviation for the section
    of <i>vec</i> from start to (start+size)-1.

    \param vec An object or array where the [] operator returns a double.

    \param start the index of the first element in the data section
    
    \param size the size of the data section

    \param stats The result (mean and standard deviation) is returned in
                 this argument.

   */
  static void bellStats(const T &vec,
			const size_t start,
			const size_t size,
			bellInfo &stats )
  {
    stats.mean( 0.0 );
    stats.stdDev( 0.0 );
    int i;
    
    // calculate the mean (a.k.a average)
    double sum = 0.0;
    for (i = start; i < start + size; i++) {
      sum = sum + vec[i];
    }
    double mean = sum / static_cast<double>(size);
    
    // calculate the standard deviation sum
    double stdDevSum = 0;
    double x;
    for (i = start; i < start + size; i++) {
      x = vec[i] - mean;
      stdDevSum = stdDevSum + (x * x);
    }
    double variance = stdDevSum / static_cast<double>(size-1);
    double stdDev = sqrt( variance );
    
    stats.mean( mean );
    stats.stdDev( stdDev );
    
  } // bellStats


  static double soft_thresh( const double val, const double tau )
  {
    double sign = 1.0;
    double absval = fabs( val );

    if (val < 0)
      sign = -1.0;

    double new_val;

    if (absval < tau) {
      new_val = 0.0;
    }
    else {
      new_val = sign * (absval - tau);
    }
    
    return new_val;
  }


  /**
    Calculate a wavelet threshold and apply it to the wavelet 
    coefficients using soft thresholding.
   */
  static void thresh(T &vec, 
		     const size_t start, 
		     const size_t size )
  {
    bellInfo info;

    bellStats( vec, start, size, info );

    double stdDev = info.stdDev();
    double Tau = sqrt( 2 * log( size ) ) * stdDev;

    size_t i;
    for (i = start; i < start + size; i++) {
      vec[ i ] = soft_thresh( vec[i], Tau );
    }
  } // thresh


  /**
    Print the mean and standard deviation for the wavelet coefficient
    bands whose length is greater than 32.

    \param vec An object or array where the [] operator returns a double
    
    \param N the number of data elements (doubles) in <i>vec</i>
   */
  static void printHarmonicStdDev( const T &vec, const size_t N)
  {
    const size_t minSize = 32;
    size_t band = N / 2;

    bellInfo stats;

    while (band >= minSize ) {
      bellStats( vec, band, band, stats );
      printf("coefficient band %3d: mean = %7.6f, stddev = %7.6f\n",
	     band, stats.mean(), stats.stdDev() );
      band = band / 2;
    } // while
  } // printHarmonicStdDev

}; // wavethresh


#endif