Dft.h

Go to the documentation of this file.

// $Id: Dft.h,v 1.2 2010/02/23 17:19:57 eml Exp $

#ifndef GCP_UTIL_DFT_H
#define GCP_UTIL_DFT_H

#define COMPLEX_AV

#include <iostream>
#include <vector>

#include "gcp/util/common/CircBuf.h"
#include "gcp/util/common/Frequency.h"
#include "gcp/util/common/TimeVal.h"

#include <fftw3.h>

#define APOD_FN(fn) double (fn)(unsigned iSamp, unsigned n)

namespace gcp {
  namespace util {

    class Dft {
    public:

      enum Apodization {
    APOD_RECTANGLE = 0, // Rectangular (Dirichlet) window
    APOD_TRIANGLE  = 1, // Triangle apodization (goes to zero at edges)
    APOD_HAMMING   = 2, // Hamming window
    APOD_HANN      = 3, // Hann window
    APOD_COS       = 4, // Cosine window
    APOD_SINC      = 5, // Sinc window
      };

      Dft(int n, bool optimize=true, Apodization apod=APOD_RECTANGLE);

      Dft(const Dft& objToBeCopied);

      Dft(Dft& objToBeCopied);

      void operator=(const Dft& objToBeAssigned);

      void operator=(Dft& objToBeAssigned);

      friend std::ostream& operator<<(std::ostream& os, Dft& obj);

      virtual ~Dft();

      //-----------------------------------------------------------------------
      // Set methods
      //-----------------------------------------------------------------------

      void resize(unsigned int);

      void setTimeRes(TimeVal tVal);

      void setAverage(bool doAv);

      void setVectorAverage(bool vecAv);

      void setApodizationType(Apodization apod);

      //-----------------------------------------------------------------------
      // Query methods
      //-----------------------------------------------------------------------

      unsigned inputSize();

      unsigned transformSize();

      static unsigned transformSize(unsigned n);

      Frequency getFrequencyResolution();

      // Static method to return the same

      static Frequency getFrequencyResolution(unsigned npt, TimeVal timeRes);

      std::vector<double> powerSpectrum(double* data);

      std::vector<double> abs();

      void abs(float* outputArray, 
           bool& first, double& min, double& max, bool doLog=false);

      void updateMinMax(unsigned i, double& val, bool doLog, 
            bool& first, double& min, double& max);

      //-----------------------------------------------------------------------
      // Manual methods for using this class
      //-----------------------------------------------------------------------

      void fillInputArray(double* data);

      void computeTransform();

      double* getInputData();

      fftw_complex* getTransform();

      //-----------------------------------------------------------------------
      // Automatic methods for computing transforms
      //-----------------------------------------------------------------------

      void pushSample(double sample); 

      bool transformIsReady();

      //-----------------------------------------------------------------------
      // Apodization functions
      //-----------------------------------------------------------------------

      static APOD_FN(apodRectangle);
      static APOD_FN(apodTriangle);
      static APOD_FN(apodHamming);
      static APOD_FN(apodHann);
      static APOD_FN(apodCos);
      static APOD_FN(apodSinc);

    private:

      void computePlan(unsigned flag); 

      void addToAverage();

      // Remove the mean from the input array

      void removeMean();

      // Assert a pending apodization type

      void assertApodizationType(Apodization apod);

      // Return the apodization coefficient for the current sample

      double apodizationCoefficient(unsigned iSamp);

      void copySamples();

    private:

      // Fundamental variables

      int n_;                       // The size transform we will compute
      TimeVal tRes_;                // The time resolution of the samples
      bool haveRes_;                // True once the resolution is set
      Apodization currentApodType_; // The current apodization type
      Apodization pendingApodType_; // The pending apodization type
      bool apodTypePending_;        // An apodization command was received,
                        // that should be put into effect the next
                        // time the sample counter rolls over
      APOD_FN(*apodFn_);            // The current apodization function

      double mean_;
      bool subtractMean_;

      // FFTW specific variables

      double* in_;          // The input array to be transformed
      fftw_complex* out_;   // The output of the transformed
      fftw_plan plan_;      // Instructions to fftw for the best method
                // to FFT
      bool optimize_;       // True if we want fftw to perform expensive
                // tests to determine the optimal FFT
                // method.  Note that these tests will only
                // be done on calls to resize()

      // Averaging variables

      unsigned nAv_;        // The number of transforms currently in
                // the running average
      bool doAv_;           // If true, average the transforms

      bool vecAv_;          // If true, vector-average the transforms.
                // If false, rms-average the power spectra
      
      bool linear_;         // If true, output data in linear (rather
                // than log) units

      fftw_complex* outAv_; // The averaged output of the transform
      std::valarray<double> absBuf_; // The average output of the
                     // power spectra

      // Auto-transform variables

      std::valarray<double> sampBuf_; // Buffer into which we store
                      // samples
      std::valarray<double> apodBuf_; // Buffer into which we store
                      // apodization coefficients
#if 0
      CircBuf<double> buf_;     // A circular buffer for storing
                    // on-the-fly samples
      CircBuf<double> apodBuf_; // A circular buffer for storing
                    // on-the-fly samples
#endif
      unsigned nSinceLastTransform_;
      bool transformIsReady_;

      
    }; // End class Dft

  } // End namespace util
} // End namespace gcp


#endif // End #ifndef GCP_UTIL_DFT_H