---

support.h

Go to the documentation of this file.

#ifndef __SUPPORT_FFT__
#define __SUPPORT_FFT__

#include <errno.h>
#include <stdio.h>

#include <sstream>
#include <string>

#include "fftw3.h"

#include "msg_stream.h"


namespace FFT{

  /*
     #################
     # class Support #
     #################
  */


  
  class Support {
    
  public:

    typedef unsigned int  size_type;

    typedef double        real_type;
    typedef fftw_complex  complex_type;
    
    typedef real_type*    space_data_type;
    typedef complex_type* frequency_data_type;
    typedef void*         buffer_type;
    
    
    inline Support(const size_type w,
                   const size_type h);

    inline ~Support();



    inline size_type x_space_size() const;
    inline size_type y_space_size() const;


    inline size_type x_frequency_size() const;
    inline size_type y_frequency_size() const;
   
    inline buffer_type create_buffer() const;

    inline void destroy_buffer(buffer_type buffer) const;



    
    template <typename Real_data_2D_array>
    void save_space_data(Real_data_2D_array* const array) const;

    template <typename Real_data_2D_array>
    void load_space_data(const Real_data_2D_array& array);


    
    
    template <typename Complex_data_2D_array>
    void save_frequency_data(Complex_data_2D_array* const array) const;

    template <typename Complex_data_2D_array>
    void load_frequency_data(const Complex_data_2D_array& array);


    
    
    inline void save_space_data_into_buffer(buffer_type buffer) const;

    inline void load_space_data_from_buffer(buffer_type buffer);
    



    inline void save_frequency_data_into_buffer(buffer_type buffer) const;

    inline void load_frequency_data_from_buffer(buffer_type buffer);


    

    inline void space_to_frequency();

    inline void frequency_to_space();

    


    inline void multiply_frequency_data_by_buffer(buffer_type buffer,
                                                  const bool  scale = false);




    

    // ####################
    // # Static functions #
    // ####################


    static inline void set_fftw_flags(const unsigned flags);

    
    static inline void set_wisdom_file(const std::string file_name);

    
  private:


    fftw_plan forward, backward;

    space_data_type     space_data;

    frequency_data_type frequency_data;


    
    const size_type width;
    const size_type height;


    
    const size_type support_w;
    const size_type support_h;
    const size_type mem_size;
    
    static inline size_type index(const size_type x,
                                  const size_type y,
                                  const size_type w,
                                  const size_type h);
   

    // ####################
    // # Static functions #
    // ####################


    static unsigned fftw_flags;

    static std::string wisdom_file;

    static size_type support_number;

    static inline void new_support();

    static inline void delete_support();
  };



/*
  ###########################################################
  ###########################################################
  ###########################################################
  ##############                               ##############
  ##############  I M P L E M E N T A T I O N  ##############
  ##############                               ##############
  ###########################################################
  ###########################################################
  ###########################################################
*/


  

  // ##################
  // # Static section #
  // ##################

  // Static values initialization
  
  unsigned           Support::fftw_flags     = FFTW_EXHAUSTIVE;
  std::string        Support::wisdom_file    = std::string("");
  Support::size_type Support::support_number = 0;

  
  

  // ######################
  // # Non-static section #
  // ######################


  void Support::set_fftw_flags(const unsigned flags){
    fftw_flags = flags;
  }

  void Support::set_wisdom_file(const std::string file_name){
    wisdom_file = file_name;
  }



  void Support::new_support(){
    support_number++;

    // Job is done if not the first support or if not file support required.
    if ((support_number > 1)||(wisdom_file.size() == 0)){
      return;
    }

    
    // If it is the first support, try to load the wisdom data.

    
    // Opening the wisdom file.
    FILE* f = fopen(wisdom_file.c_str(),"r");
    
    if (f != NULL){ // System says "ok"
      if(fftw_import_wisdom_from_file(f) == 0){ // But FFTW says "error".
        Message::warning<<"error while reading FFTW wisdom file \""
                        <<wisdom_file<<"\""<<Message::done;
      }
      fclose(f);
    }
    else if (errno != ENOENT){ // System says "error".
      Message::warning<<"error while opening (read mode) FFTW wisdom file \""
                      <<wisdom_file<<"\""<<Message::done;
    }
    // Wisdom is now initialized.
  }


  
  void Support::delete_support(){
    support_number--;

    // Job is done if not the last support.
    if ((support_number > 0)||(wisdom_file.size() == 0)){
      return;
    }

    // If it is the last support, try to save the wisdom data.

    // Save wisdom.
    FILE* f = fopen(wisdom_file.c_str(),"w");

    if (f != NULL){ // System says "ok".
      fftw_export_wisdom_to_file(f);
        fclose(f);
    }
    else{ // System says "error".
      Message::warning<<"error while opening (write mode) FFTW wisdom file \""
                      <<wisdom_file<<"\"\n"
                      <<"no file saved"<<Message::done;
    }
  }

  Support::Support(const size_type w,
                   const size_type h):
    width(w),
    height(h),
    
    support_w(w),
    support_h(h/2+1),
    mem_size(w*(h/2+1)*sizeof(complex_type)) {

    new_support();
    
    // Pointer to the space data.
    space_data = reinterpret_cast<space_data_type>(fftw_malloc(mem_size));

    // Same pointer to perform "in place" computation.
    frequency_data = reinterpret_cast<frequency_data_type>(space_data);

    
    // We make the best measurement possible because we use the wisdom system.
    forward  = fftw_plan_dft_r2c_2d(width,height,
                                    space_data,frequency_data,
                                    fftw_flags);
    
    backward = fftw_plan_dft_c2r_2d(width,height,
                                    frequency_data,space_data,
                                    fftw_flags);
  }





  
  Support::~Support(){

    delete_support();
    
    // Free memory.
    fftw_destroy_plan(forward);
    fftw_destroy_plan(backward);

    fftw_free(space_data);
  }



  
  Support::size_type Support::x_space_size() const{
    return width;
  }


  
  Support::size_type Support::y_space_size() const{
    return height;
  }


  Support::size_type Support::x_frequency_size() const{
    return support_w;
  }


  Support::size_type Support::y_frequency_size() const{
    return support_h;
  }


  
  Support::size_type Support::index(const size_type x,
                                    const size_type y,
                                    const size_type w,
                                    const size_type h){
    return y + x*h;
  }


  
  
  Support::buffer_type Support::create_buffer() const{

    return fftw_malloc(mem_size);
  }

  void Support::destroy_buffer(buffer_type buffer) const{
    
    fftw_free(buffer);
  }

  

  
  template <typename Real_data_2D_array>
  void Support::save_space_data(Real_data_2D_array* const array) const{

    array->resize(width,height);

    for(size_type x=0;x<width;x++){
      for(size_type y=0;y<height;y++){
        const size_type i = index(x,y,support_w,2*support_h);

        (*array)(x,y) = space_data[i];
      }
    }
  }

  template <typename Real_data_2D_array>
  void Support::load_space_data(const Real_data_2D_array& array){

    if ((array.x_size()!=width)||(array.y_size()!=height)){
      Message::error<<"Support::load_space_data : size does not match"
                    <<Message::done;
    }

    for(size_type x=0;x<width;x++){
      for(size_type y=0;y<height;y++){
        const size_type i = index(x,y,support_w,2*support_h);

        space_data[i] = array(x,y);
      }
    }

  }



  
  template <typename Complex_data_2D_array>
  void Support::save_frequency_data(Complex_data_2D_array* const array) const{

    typedef typename Complex_data_2D_array::value_type cplx_type;
    
    array->resize(support_w,support_h);

    for(size_type x=0;x<support_w;x++){
      for(size_type y=0;y<support_h;y++){
        const size_type i = index(x,y,support_w,support_h);

        (*array)(x,y) = cplx_type(frequency_data[i][0],frequency_data[i][1]);
      }
    }
  }

  template <typename Complex_data_2D_array>
  void Support::load_frequency_data(const Complex_data_2D_array& array){

    if ((array.x_size()!=support_w)||(array.y_size()!=support_h)){
      Message::error<<"Support::load_frequency_data : size does not match"
                    <<Message::done;
    }

    for(size_type x=0;x<support_w;x++){
      for(size_type y=0;y<support_h;y++){
        const size_type i = index(x,y,support_w,support_h);

        frequency_data[i][0] = array(x,y).real();
        frequency_data[i][1] = array(x,y).imag();
      }
    }

  }



  
  
  void Support::save_space_data_into_buffer(buffer_type buffer) const{
    memcpy(buffer,space_data,mem_size);
  }

  void Support::load_space_data_from_buffer(buffer_type buffer){
    memcpy(space_data,buffer,mem_size);
  }

  

  
  void Support::save_frequency_data_into_buffer(buffer_type buffer) const{
    memcpy(buffer,frequency_data,mem_size);
  }

  void Support::load_frequency_data_from_buffer(buffer_type buffer){
    memcpy(frequency_data,buffer,mem_size);
  }





  void Support::space_to_frequency(){
    fftw_execute(forward);
  }


  void Support::frequency_to_space(){
    fftw_execute(backward);
  }




  void Support::multiply_frequency_data_by_buffer(buffer_type buffer,
                                                  const bool  scale){
    
    const size_type size = support_w*support_h;
    const frequency_data_type freq_buffer =
      reinterpret_cast<frequency_data_type>(buffer);

    const real_type s = scale ? 1.0/(width*height) : 1.0;
    
    for(size_type i=0;i<size;i++){

      const double a_re = frequency_data[i][0];
      const double a_im = frequency_data[i][1];
      
      const double b_re = freq_buffer[i][0];
      const double b_im = freq_buffer[i][1];
      
      frequency_data[i][0] = (a_re*b_re - a_im*b_im)*s;
      frequency_data[i][1] = (a_re*b_im + a_im*b_re)*s;
      
    }
    
  }

}


#endif

---