FFT.hh

/* 
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS 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 Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/


#ifndef __FFT_hh__
#define __FFT_hh__

#include <bias_config.h>
#include <Base/Common/W32Compat.hh>
#include <math.h>

#include <Base/Debug/Debug.hh>
#include <Base/Math/Vector.hh>

#include <complex>
#include <iostream>
#include <fftw3.h>

namespace BIAS {


/**
 * @class FFT
     @ingroup g_mathalgo
     @brief Wrapper to the fftw3 library.
     @note Use fast as follows:
    - call Init() with appropriate size
    - call GetInput() and fill it with data
    - call Compute()
    - get Results with GetPhase(), GetMagnitude() or GetRes()
    
    a slow interface for Vector<double> also exists
    an example exists in Examples/ExampleFFT, see ExampleFFT.cpp
    @author woelk 03/2004 
*/
  class BIASMathAlgo_EXPORT FFT : public Debug
  {
   
  public:
    FFT();
    ~FFT();

    /** initializes the internal variables. This may take some time 
        (i.e. several seconds) since some measurements are done */
    void Init(int size);

    /** cal this before a second call to Init() */
    void Release();
    
    /** returns the internal data. Do _not_ delete */
    inline double *GetInput()
    { return _in; }

    /** slow, memcopys in to internal input data */
    void Compute(double *in);

    /** slow, memcopys in to internal input data */
    void Compute(Vector<double> &in);

    /** input must be filled before with GetInput */
    inline void Compute()
    { fftw_execute(_p); }

    /** returns the resulting phase, phase must be of length 
        size/2+1 (rounded downwards) */
    void GetPhase(double *phase);
    
    /** returns the resulting magnitude, mag must be of length 
       size/2+1 (rounded downwards) */
    void GetMagnitude(double *mag);

    /** returns the internal data. Do _not_ delete 
        res is of length size/2+1 */
    inline std::complex<double>* GetResult()
    { return reinterpret_cast<std::complex<double> *>(_out); }

    /** returns the magnitude */
    void GetPhase(Vector<double>& phase);

    /** returns the phase */
    void GetMagnitude(Vector<double>& mag);

    void DumpResult(std::ostream& os=std::cout);

  protected:
    fftw_complex *_out;
    double *_in;
    fftw_plan _p;
    size_t _InputSize; // in bytes
    int _Size; // length of input
    int _OutNum; // length of output = _Size/2 + 1
  };

} // namespace

#endif // __FFT_hh__