---

convolution.h

Go to the documentation of this file.

#ifndef __CONVOLUTION_FFT__
#define __CONVOLUTION_FFT__

#include <algorithm>

#include "fftw3.h"

#include "fft/support.h"



namespace FFT{

  
  

  template<typename Function>
  inline void convolve(const Function& f,
                       const Function& g,
                       Function* const h);

  
  template<typename Function>
  inline void convolve(const Function& f,
                       const Function& g,
                       Function* const h,
                       Support&        support);

  
  template<typename Function>
  void convolve(const Support::buffer_type& F,
                const Function& g,
                Function* const h,
                Support&        support);

  
  template<typename Function>
  void move_center(const Function&    f,
                   const unsigned int x,
                   const unsigned int y,
                   Function* const    result);



  


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

  
  template<typename Function>
  void convolve(const Support::buffer_type& F,
                const Function& g,
                Function* const h,
                Support&        support){

    // ################
    // ### FFT of g ###
    // ################

    // Store g data in the fftw structures.    
    support.load_space_data(g);
    
    // Compute the FFT.
    support.space_to_frequency();



    
    // Multiplication "in place" in the frequency space.    
    support.multiply_frequency_data_by_buffer(F,true);


    

    // ###################
    // ### Inverse FFT ###
    // ###################
    
    // Compute inverse FFT transform.
    support.frequency_to_space();

    // Return data into the original structures    
    support.save_space_data(h);
  }


  template<typename Function>
  void convolve(const Function& f,
                const Function& g,
                Function* const h,
                Support&  support){

    // ################
    // ### FFT of f ###
    // ################

    // We need a temporary buffer.
    Support::buffer_type buffer = support.create_buffer();

    // Store f data in the fftw structures.
    support.load_space_data(f);
    
    // Compute the FFT.
    support.space_to_frequency();
    
    // Copy the result in a safe place.
    support.save_frequency_data_into_buffer(buffer);


    convolve(buffer,g,h,support);
    
    support.destroy_buffer(buffer);
    
//     // ################
//     // ### FFT of g ###
//     // ################

//     // Store g data in the fftw structures.    
//     support.load_space_data(g);
    
//     // Compute the FFT.
//     support.space_to_frequency();



    
//     // Multiplication "in place" in the frequency space.    
//     support.multiply_frequency_data_by_buffer(buffer,true);


    

//     // ###################
//     // ### Inverse FFT ###
//     // ###################
    
//     // Compute inverse FFT transform.
//     support.frequency_to_space();

//     // Return data into the original structures    
//     support.save_space_data(h);

  }

  

  template<typename Function>
  void convolve(const Function& f,
                const Function& g,
                Function* const h){

    Support s(f.x_size(),f.y_size());

    convolve(f,g,h,s);
  }

  
  template<typename Function>
  void move_center(const Function&    f,
                   const unsigned int x_center,
                   const unsigned int y_center,
                   Function* const    result){

    typedef unsigned int size_type;

    const size_type width  = f.x_size();
    const size_type height = f.y_size();

    Function g(width,height);
    
    for(size_type x=0;x<width;x++){
      for(size_type y=0;y<height;y++){
        (*result)(x,y) = f((x + width - x_center) % width,
                           (y + height - y_center) % height);
      }
    }
  }

}

#endif

---