src/fft_fftw.cpp

   1 #include "fft_fftw.h"
   2
   3 #include "fft.h"
   4
   5 #ifdef OPENMP
   6   #include <omp.h>
   7 #endif
   8
   9 #ifdef ASYNC
  10 #define FFTW_PLAN_WITH_THREADS() fftw_plan_with_nthreads(m_num_threads);
  11 #elif OPENMP
  12 #define FFTW_PLAN_WITH_THREADS() fftw_plan_with_nthreads(omp_get_max_threads());
  13 #else
  14 #define FFTW_PLAN_WITH_THREADS()
  15 #endif
  16
  17 Fftw::Fftw()
  18     : m_num_threads(1)
  19 {
  20 }
  21
  22 Fftw::Fftw(int num_threads)
  23     : m_num_threads(num_threads)
  24 {
  25 }
  26
  27 void Fftw::init(unsigned width, unsigned height, unsigned num_of_feats, unsigned num_of_scales)
  28 {
  29     m_width = width;
  30     m_height = height;
  31     m_num_of_feats = num_of_feats;
  32     m_num_of_scales = num_of_scales;
  33
  34 #if defined(ASYNC) || defined(OPENMP)
  35     fftw_init_threads();
  36 #endif //OPENMP
  37
  38 #ifndef CUFFTW
  39     std::cout << "FFT: FFTW" << std::endl;
  40 #else
  41     std::cout << "FFT: cuFFTW" << std::endl;
  42 #endif
  43
  44     {
  45     cv::Mat in_f = cv::Mat::zeros(m_height, m_width, CV_32FC1);
  46     ComplexMat out_f(m_height, m_width / 2 + 1, 1);
  47     plan_f = fftwf_plan_dft_r2c_2d(m_height, m_width,
  48                                    reinterpret_cast<float*>(in_f.data),
  49                                    reinterpret_cast<fftwf_complex*>(out_f.get_p_data()),
  50                                    FFTW_MEASURE);
  51     }
  52
  53     {
  54         cv::Mat in_fw = cv::Mat::zeros(m_height * m_num_of_feats, m_width, CV_32F);
  55         ComplexMat out_fw(m_height, m_width / 2 + 1, m_num_of_feats);
  56         float *in = reinterpret_cast<float*>(in_fw.data);
  57         fftwf_complex *out = reinterpret_cast<fftwf_complex*>(out_fw.get_p_data());
  58         int rank = 2;
  59         int n[] = {(int)m_height, (int)m_width};
  60         int howmany = m_num_of_feats;
  61         int idist = m_height*m_width, odist = m_height*(m_width/2+1);
  62         int istride = 1, ostride = 1;
  63         int *inembed = NULL, *onembed = NULL;
  64
  65         FFTW_PLAN_WITH_THREADS();
  66         plan_fw = fftwf_plan_many_dft_r2c(rank, n, howmany,
  67                                           in, inembed, istride, idist,
  68                                           out, onembed, ostride, odist,
  69                                           FFTW_MEASURE);
  70     }
  71     if(num_of_scales > 1)
  72     {
  73         cv::Mat in_all = cv::Mat::zeros(m_height * (num_of_scales*m_num_of_feats), m_width, CV_32F);
  74         ComplexMat out_all(m_height, m_width / 2 + 1, num_of_scales*m_num_of_feats);
  75         float *in = reinterpret_cast<float*>(in_all.data);
  76         fftwf_complex *out = reinterpret_cast<fftwf_complex*>(out_all.get_p_data());
  77         int rank = 2;
  78         int n[] = {(int)m_height, (int)m_width};
  79         int howmany = num_of_scales*m_num_of_feats;
  80         int idist = m_height*m_width, odist = m_height*(m_width/2+1);
  81         int istride = 1, ostride = 1;
  82         int *inembed = NULL, *onembed = NULL;
  83
  84         FFTW_PLAN_WITH_THREADS();
  85         plan_fw_all_scales = fftwf_plan_many_dft_r2c(rank, n, howmany,
  86                                                      in,  inembed, istride, idist,
  87                                                      out, onembed, ostride, odist,
  88                                                      FFTW_MEASURE);
  89     }
  90
  91     {
  92         ComplexMat in_i(m_height,m_width,m_num_of_feats);
  93         cv::Mat out_i = cv::Mat::zeros(m_height, m_width, CV_32FC(m_num_of_feats));
  94         fftwf_complex *in = reinterpret_cast<fftwf_complex*>(in_i.get_p_data());
  95         float *out = reinterpret_cast<float*>(out_i.data);
  96         int rank = 2;
  97         int n[] = {(int)m_height, (int)m_width};
  98         int howmany = m_num_of_feats;
  99         int idist = m_height*(m_width/2+1), odist = 1;
 100         int istride = 1, ostride = m_num_of_feats;
 101         int inembed[] = {(int)m_height, (int)m_width/2+1}, *onembed = n;
 102
 103         FFTW_PLAN_WITH_THREADS();
 104         plan_i_features = fftwf_plan_many_dft_c2r(rank, n, howmany,
 105                                                   in,  inembed, istride, idist,
 106                                                   out, onembed, ostride, odist,
 107                                                   FFTW_MEASURE);
 108     }
 109
 110     {
 111         ComplexMat in_i1(m_height,m_width,1);
 112         cv::Mat out_i1 = cv::Mat::zeros(m_height, m_width, CV_32FC1);
 113         fftwf_complex *in = reinterpret_cast<fftwf_complex*>(in_i1.get_p_data());
 114         float *out = reinterpret_cast<float*>(out_i1.data);
 115         int rank = 2;
 116         int n[] = {(int)m_height, (int)m_width};
 117         int howmany = 1;
 118         int idist = m_height*(m_width/2+1), odist = 1;
 119         int istride = 1, ostride = 1;
 120         int inembed[] = {(int)m_height, (int)m_width/2+1}, *onembed = n;
 121
 122         FFTW_PLAN_WITH_THREADS();
 123         plan_i_1ch = fftwf_plan_many_dft_c2r(rank, n, howmany,
 124                                              in,  inembed, istride, idist,
 125                                              out, onembed, ostride, odist,
 126                                              FFTW_MEASURE);
 127     }
 128 }
 129
 130 void Fftw::set_window(const cv::Mat &window)
 131 {
 132     m_window = window;
 133 }
 134
 135 ComplexMat Fftw::forward(const cv::Mat &input)
 136 {
 137     cv::Mat complex_result(m_height, m_width / 2 + 1, CV_32FC2);
 138
 139     fftwf_execute_dft_r2c(plan_f,reinterpret_cast<float*>(input.data),reinterpret_cast<fftwf_complex*>(complex_result.data));
 140
 141     return ComplexMat(complex_result);
 142 }
 143
 144 ComplexMat Fftw::forward_window(const std::vector<cv::Mat> &input)
 145 {
 146     int n_channels = input.size();
 147     cv::Mat in_all(m_height * n_channels, m_width, CV_32F);
 148     for (int i = 0; i < n_channels; ++i) {
 149         cv::Mat in_roi(in_all, cv::Rect(0, i*m_height, m_width, m_height));
 150         in_roi = input[i].mul(m_window);
 151     }
 152     ComplexMat result(m_height, m_width/2 + 1, n_channels);
 153
 154     float *in = reinterpret_cast<float*>(in_all.data);
 155     fftwf_complex *out = reinterpret_cast<fftwf_complex*>(result.get_p_data());
 156
 157     if (n_channels <= 44)
 158         fftwf_execute_dft_r2c(plan_fw, in, out);
 159     else
 160         fftwf_execute_dft_r2c(plan_fw_all_scales, in, out);
 161
 162     return result;
 163 }
 164
 165 cv::Mat Fftw::inverse(const ComplexMat &inputf)
 166 {
 167     int n_channels = inputf.n_channels;
 168     cv::Mat real_result(m_height, m_width, CV_32FC(n_channels));
 169     fftwf_complex *in = reinterpret_cast<fftwf_complex*>(inputf.get_p_data());
 170     float *out = reinterpret_cast<float*>(real_result.data);
 171
 172     if(n_channels != 1)
 173         fftwf_execute_dft_c2r(plan_i_features, in, out);
 174     else
 175         fftwf_execute_dft_c2r(plan_i_1ch, in, out);
 176
 177     return real_result/(m_width*m_height);
 178 }
 179
 180 Fftw::~Fftw()
 181 {
 182   fftwf_destroy_plan(plan_f);
 183   fftwf_destroy_plan(plan_fw);
 184   fftwf_destroy_plan(plan_fw_all_scales);
 185   fftwf_destroy_plan(plan_i_features);
 186   fftwf_destroy_plan(plan_i_1ch);
 187 }