Line data    Source code

       1             : /**
       2             :  * \file FFT.cpp
       3             :  */
       4             : 
       5             : #include "FFT.h"
       6             : 
       7             : #include <ATK/config.h>
       8             : #if ATK_USE_FFTW == 1
       9             : #include <fftw3.h>
      10             : #endif
      11             : #if ATK_USE_IPP == 1
      12             : #include <ipp.h>
      13             : #endif
      14             : 
      15             : #include <gsl/gsl>
      16             : 
      17             : #include <algorithm>
      18             : #include <complex>
      19             : #include <cstdlib>
      20             : 
      21             : namespace ATK
      22             : {
      23             :   template<class DataType_>
      24             :   class FFT<DataType_>::FFTImpl
      25             :   {
      26             : #if ATK_USE_FFTW == 1
      27             :     fftw_plan fft_plan;
      28             :     fftw_plan fft_reverse_plan;
      29             :     fftw_complex* input_data;
      30             :     fftw_complex* output_freqs;
      31             : #endif
      32             : #if ATK_USE_IPP == 1
      33             :     Ipp64fc *pSrc;
      34             :     Ipp64fc *pDst;
      35             :     IppsFFTSpec_C_64fc* pFFTSpec;
      36             :     IppsDFTSpec_C_64fc* pDFTSpec;
      37             :     Ipp8u* pFFTSpecBuf;
      38             :     Ipp8u* pFFTInitBuf;
      39             :     Ipp8u* pFFTWorkBuf;
      40             :     bool power_of_two;
      41             : #endif
      42             :   public:
      43         332 :     FFTImpl()
      44             : #if ATK_USE_FFTW == 1
      45         332 :     :fft_plan(nullptr), fft_reverse_plan(nullptr), input_data(nullptr), output_freqs(nullptr)
      46             : #endif
      47             : #if ATK_USE_IPP == 1
      48             :     :pSrc(nullptr), pDst(nullptr), pFFTSpec(nullptr), pDFTSpec(nullptr), pFFTSpecBuf(nullptr), pFFTInitBuf(nullptr), pFFTWorkBuf(nullptr), power_of_two(false)
      49             : #endif
      50             :     {
      51         332 :     }
      52             :     
      53         332 :     ~FFTImpl()
      54             :     {
      55             : #if ATK_USE_FFTW == 1
      56         332 :       fftw_free(input_data);
      57         332 :       fftw_free(output_freqs);
      58         332 :       fftw_destroy_plan(fft_plan);
      59         332 :       fftw_destroy_plan(fft_reverse_plan);
      60             : #endif
      61             : #if ATK_USE_IPP == 1
      62             :       if (pSrc)
      63             :       {
      64             :         ippFree(pSrc);
      65             :       }
      66             :       if (pDst)
      67             :       {
      68             :         ippFree(pDst);
      69             :       }
      70             :       if (pFFTSpec)
      71             :       {
      72             :         ippFree(pFFTSpec);
      73             :       }
      74             :       if (pFFTInitBuf)
      75             :       {
      76             :         ippFree(pFFTInitBuf);
      77             :       }
      78             :       if (pFFTWorkBuf)
      79             :       {
      80             :         ippFree(pFFTWorkBuf);
      81             :       }
      82             : #endif
      83         332 :     }
      84             : 
      85         331 :     void set_size(gsl::index size)
      86             :     {
      87             : #if ATK_USE_FFTW == 1
      88         331 :       fftw_free(input_data);
      89         331 :       fftw_free(output_freqs);
      90         331 :       fftw_destroy_plan(fft_plan);
      91         331 :       fftw_destroy_plan(fft_reverse_plan);
      92             :       
      93         331 :       input_data = fftw_alloc_complex(size);
      94         331 :       output_freqs = fftw_alloc_complex(size);
      95             :       
      96         331 :       fft_plan = fftw_plan_dft_1d(size, input_data, output_freqs, FFTW_FORWARD, FFTW_ESTIMATE);
      97         331 :       fft_reverse_plan = fftw_plan_dft_1d(size, output_freqs, input_data, FFTW_BACKWARD, FFTW_ESTIMATE);
      98             : #endif
      99             : #if ATK_USE_IPP == 1
     100             :       power_of_two = ((size & (size - 1)) == 0);
     101             :       if (pSrc)
     102             :       {
     103             :         ippFree(pSrc);
     104             :       }
     105             :       if (pDst)
     106             :       {
     107             :         ippFree(pDst);
     108             :       }
     109             :       if (pFFTSpec)
     110             :       {
     111             :         ippFree(pFFTSpec);
     112             :       }
     113             :       pFFTSpec = nullptr;
     114             :       if (pDFTSpec)
     115             :       {
     116             :         ippFree(pDFTSpec);
     117             :       }
     118             :       pDFTSpec = nullptr;
     119             :       if (pFFTSpecBuf)
     120             :       {
     121             :         ippFree(pFFTSpecBuf);
     122             :       }
     123             :       if (pFFTInitBuf)
     124             :       {
     125             :         ippFree(pFFTInitBuf);
     126             :       }
     127             :       if (pFFTWorkBuf)
     128             :       {
     129             :         ippFree(pFFTWorkBuf);
     130             :       }
     131             :       int sizeFFTSpec, sizeFFTInitBuf, sizeFFTWorkBuf;
     132             :       if (power_of_two)
     133             :       {
     134             :         ippsFFTGetSize_C_64fc(static_cast<int>(std::lround(std::log(size) / std::log(2))), IPP_FFT_DIV_FWD_BY_N, ippAlgHintAccurate, &sizeFFTSpec, &sizeFFTInitBuf, &sizeFFTWorkBuf);
     135             :       }
     136             :       else
     137             :       {
     138             :         ippsDFTGetSize_C_64fc(static_cast<int>(size), IPP_FFT_DIV_FWD_BY_N, ippAlgHintAccurate, &sizeFFTSpec, &sizeFFTInitBuf, &sizeFFTWorkBuf);
     139             :       }
     140             :       
     141             :       pFFTSpecBuf = ippsMalloc_8u(sizeFFTSpec);
     142             :       pFFTInitBuf = ippsMalloc_8u(sizeFFTInitBuf);
     143             :       pFFTWorkBuf = ippsMalloc_8u(sizeFFTWorkBuf);
     144             :       pSrc = ippsMalloc_64fc(static_cast<int>(size));
     145             :       pDst = ippsMalloc_64fc(static_cast<int>(size));
     146             :       if (power_of_two)
     147             :       {
     148             :         ippsFFTInit_C_64fc(&pFFTSpec, static_cast<int>(std::lround(std::log(size) / std::log(2))), IPP_FFT_DIV_FWD_BY_N, ippAlgHintAccurate, pFFTSpecBuf, pFFTInitBuf);
     149             :       }
     150             :       else
     151             :       {
     152             :         pDFTSpec = (IppsDFTSpec_C_64fc*)ippsMalloc_8u(sizeFFTSpec);
     153             :         ippsDFTInit_C_64fc(static_cast<int>(size), IPP_FFT_DIV_FWD_BY_N, ippAlgHintAccurate, pDFTSpec, pFFTInitBuf);
     154             :       }
     155             :       if (pFFTInitBuf)
     156             :         ippFree(pFFTInitBuf);
     157             :       pFFTInitBuf = nullptr;
     158             : #endif
     159         331 :     }
     160             : 
     161       17106 :     void process(const DataType_* input, gsl::index input_size, gsl::index size) const
     162             :     {
     163             : #if ATK_USE_FFTW == 1
     164       17106 :       double factor = static_cast<double>(size);
     165    41648000 :       for(gsl::index j = 0; j < std::min(input_size, size); ++j)
     166             :       {
     167    41630900 :         input_data[j][0] = input[j] / factor;
     168    41630900 :         input_data[j][1] = 0;
     169             :       }
     170     2134690 :       for(gsl::index j = input_size; j < size; ++j)
     171             :       {
     172     2117590 :         input_data[j][0] = 0;
     173     2117590 :         input_data[j][1] = 0;
     174             :       }
     175       17106 :       fftw_execute(fft_plan);
     176             : #endif
     177             : #if ATK_USE_IPP == 1
     178             :       for (gsl::index j = 0; j < std::min(input_size, size); ++j)
     179             :       {
     180             :         pSrc[j].re = input[j];
     181             :         pSrc[j].im = 0;
     182             :       }
     183             :       for (gsl::index j = input_size; j < size; ++j)
     184             :       {
     185             :         pSrc[j].re = 0;
     186             :         pSrc[j].im = 0;
     187             :       }
     188             :       if (power_of_two)
     189             :       {
     190             :         ippsFFTFwd_CToC_64fc(pSrc, pDst, pFFTSpec, pFFTWorkBuf);
     191             :       }
     192             :       else
     193             :       {
     194             :         ippsDFTFwd_CToC_64fc(pSrc, pDst, pDFTSpec, pFFTWorkBuf);
     195             :       }
     196             : #endif
     197       17106 :     }
     198             :     
     199       16498 :     void process_forward(std::complex<DataType_>* output, gsl::index size) const
     200             :     {
     201     4246760 :       for(gsl::index j = 0; j < size; ++j)
     202             :       {
     203             : #if ATK_USE_FFTW == 1
     204     4230260 :         output[j] = std::complex<DataType_>(output_freqs[j][0], output_freqs[j][1]);
     205             : #endif
     206             : #if ATK_USE_IPP == 1
     207             :         output[j] = std::complex<DataType_>(pDst[j].re, pDst[j].im);
     208             : #endif
     209             :       }
     210       16498 :     }
     211             :     
     212       16402 :     void process_backward(const std::complex<DataType_>* input, DataType_* output, gsl::index input_size, gsl::index size) const
     213             :     {
     214             : #if ATK_USE_FFTW == 1
     215     4213860 :       for(gsl::index j = 0; j < size; ++j)
     216             :       {
     217     4197460 :         output_freqs[j][0] = input[j].real();
     218     4197460 :         output_freqs[j][1] = input[j].imag();
     219             :       }
     220       16402 :       fftw_execute(fft_reverse_plan);
     221     4213860 :       for(gsl::index j = 0; j < input_size; ++j)
     222             :       {
     223     4197460 :         output[j] = input_data[j][0];
     224             :       }
     225             : #endif
     226             : #if ATK_USE_IPP == 1
     227             :       for (gsl::index j = 0; j < size; ++j)
     228             :       {
     229             :         pDst[j].re = input[j].real();
     230             :         pDst[j].im = input[j].imag();
     231             :       }
     232             :       if (power_of_two)
     233             :       {
     234             :         ippsFFTInv_CToC_64fc(pDst, pDst, pFFTSpec, pFFTWorkBuf);
     235             :       }
     236             :       else
     237             :       {
     238             :         ippsDFTInv_CToC_64fc(pDst, pDst, pDFTSpec, pFFTWorkBuf);
     239             :       }
     240             :       for (gsl::index j = 0; j < input_size; ++j)
     241             :       {
     242             :         output[j] = pDst[j].re;
     243             :       }
     244             : #endif
     245       16402 :     }
     246             : 
     247           0 :     void process(const std::complex<DataType_>* input, gsl::index input_size, gsl::index size) const
     248             :     {
     249             : #if ATK_USE_FFTW == 1
     250           0 :       double factor = static_cast<double>(size);
     251           0 :       for (gsl::index j = 0; j < std::min(input_size, size); ++j)
     252             :       {
     253           0 :         input_data[j][0] = std::real(input[j]) / factor;
     254           0 :         input_data[j][1] = std::imag(input[j]) / factor;
     255             :       }
     256           0 :       for (gsl::index j = input_size; j < size; ++j)
     257             :       {
     258           0 :         input_data[j][0] = 0;
     259           0 :         input_data[j][1] = 0;
     260             :       }
     261           0 :       fftw_execute(fft_plan);
     262             : #endif
     263             : #if ATK_USE_IPP == 1
     264             :       for (gsl::index j = 0; j < std::min(input_size, size); ++j)
     265             :       {
     266             :         pSrc[j].re = std::real(input[j]);
     267             :         pSrc[j].im = std::imag(input[j]);
     268             :       }
     269             :       for (gsl::index j = input_size; j < size; ++j)
     270             :       {
     271             :         pSrc[j].re = 0;
     272             :         pSrc[j].im = 0;
     273             :       }
     274             :       if (power_of_two)
     275             :       {
     276             :         ippsFFTFwd_CToC_64fc(pSrc, pDst, pFFTSpec, pFFTWorkBuf);
     277             :       }
     278             :       else
     279             :       {
     280             :         ippsDFTFwd_CToC_64fc(pSrc, pDst, pDFTSpec, pFFTWorkBuf);
     281             :       }
     282             : #endif
     283           0 :     }
     284             :     
     285           0 :     void process_backward(const std::complex<DataType_>* input, std::complex<DataType_>* output, gsl::index input_size, gsl::index size) const
     286             :     {
     287             : #if ATK_USE_FFTW == 1
     288           0 :       for (gsl::index j = 0; j < size; ++j)
     289             :       {
     290           0 :         output_freqs[j][0] = input[j].real();
     291           0 :         output_freqs[j][1] = input[j].imag();
     292             :       }
     293           0 :       fftw_execute(fft_reverse_plan);
     294           0 :       for (gsl::index j = 0; j < input_size; ++j)
     295             :       {
     296           0 :         output[j] = std::complex<DataType_>(input_data[j][0], input_data[j][1]);
     297             :       }
     298             : #endif
     299             : #if ATK_USE_IPP == 1
     300             :       for (gsl::index j = 0; j < size; ++j)
     301             :       {
     302             :         pDst[j].re = input[j].real();
     303             :         pDst[j].im = input[j].imag();
     304             :       }
     305             :       if (power_of_two)
     306             :       {
     307             :         ippsFFTInv_CToC_64fc(pDst, pDst, pFFTSpec, pFFTWorkBuf);
     308             :       }
     309             :       else
     310             :       {
     311             :         ippsDFTInv_CToC_64fc(pDst, pDst, pDFTSpec, pFFTWorkBuf);
     312             :       }
     313             :       for (gsl::index j = 0; j < input_size; ++j)
     314             :       {
     315             :         output[j] = std::complex<DataType_>(pDst[j].re, pDst[j].im);
     316             :       }
     317             : #endif
     318           0 :     }
     319             :     
     320         608 :     void get_amp(std::vector<DataType_>& amp) const
     321             :     {
     322    19760300 :       for(size_t i = 0; i < amp.size(); ++i)
     323             :       {
     324             : #if ATK_USE_FFTW == 1
     325    19759700 :         amp[i] = output_freqs[i][0] * output_freqs[i][0];
     326    19759700 :         amp[i] += output_freqs[i][1] * output_freqs[i][1];
     327    19759700 :         amp[i] = std::sqrt(amp[i]);
     328             : #endif
     329             : #if ATK_USE_IPP == 1
     330             :         amp[i] = pDst[i].re * pDst[i].re;
     331             :         amp[i] += pDst[i].im * pDst[i].im;
     332             :         amp[i] = std::sqrt(amp[i]);
     333             : #endif
     334             :       }
     335         608 :     }
     336             :     
     337           0 :     void get_angle(std::vector<DataType_>& angle) const
     338             :     {
     339           0 :       for(size_t i = 0; i < angle.size(); ++i)
     340             :       {
     341             : #if ATK_USE_FFTW == 1
     342           0 :         angle[i] = std::arg(std::complex<DataType_>(output_freqs[i][0], output_freqs[i][1]));
     343             : #endif
     344             : #if ATK_USE_IPP == 1
     345             :         angle[i] = std::arg(std::complex<DataType_>(pDst[i].re, pDst[i].im));
     346             : #endif
     347             :       }
     348           0 :     }
     349             :   };
     350             : 
     351             :   template<class DataType_>
     352         332 :   FFT<DataType_>::FFT()
     353         332 :   :size(0), impl(std::make_unique<FFTImpl>())
     354             :   {
     355         332 :   }
     356             :   
     357             : 
     358             :   template<class DataType_>
     359         332 :   FFT<DataType_>::~FFT()
     360             :   {
     361         332 :   }
     362             : 
     363             :   template<class DataType_>
     364         332 :   void FFT<DataType_>::set_size(gsl::index size_)
     365             :   {
     366         332 :     if(size == size_)
     367             :     {
     368           1 :       return;
     369             :     }
     370         331 :     size = size_;
     371         331 :     impl->set_size(size);
     372             :   }
     373             : 
     374             :   template<class DataType_>
     375           0 :   gsl::index FFT<DataType_>::get_size() const
     376             :   {
     377           0 :     return size;
     378             :   }
     379             : 
     380             :   template<class DataType_>
     381       17106 :   void FFT<DataType_>::process(const DataType_* input, gsl::index input_size) const
     382             :   {
     383       17106 :     impl->process(input, input_size, size);
     384       17106 :   }
     385             :   
     386             :   template<class DataType_>
     387       16498 :   void FFT<DataType_>::process_forward(const DataType_* input, std::complex<DataType_>* output, gsl::index input_size) const
     388             :   {
     389       16498 :     process(input, input_size);
     390       16498 :     impl->process_forward(output, size);
     391       16498 :   }
     392             : 
     393             :   template<class DataType_>
     394       16402 :   void FFT<DataType_>::process_backward(const std::complex<DataType_>* input, DataType_* output, gsl::index input_size) const
     395             :   {
     396       16402 :     impl->process_backward(input, output, input_size, size);
     397       16402 :   }
     398             : 
     399             :   template<class DataType_>
     400           0 :   void FFT<DataType_>::process(const std::complex<DataType_>* input, gsl::index input_size) const
     401             :   {
     402           0 :     impl->process(input, input_size, size);
     403           0 :   }
     404             : 
     405             :   template<class DataType_>
     406           0 :   void FFT<DataType_>::process_forward(const std::complex<DataType_>* input, std::complex<DataType_>* output, gsl::index input_size) const
     407             :   {
     408           0 :     process(input, input_size);
     409           0 :     impl->process_forward(output, size);
     410           0 :   }
     411             : 
     412             :   template<class DataType_>
     413           0 :   void FFT<DataType_>::process_backward(const std::complex<DataType_>* input, std::complex<DataType_>* output, gsl::index input_size) const
     414             :   {
     415           0 :     impl->process_backward(input, output, input_size, size);
     416           0 :   }
     417             : 
     418             :   template<class DataType_>
     419         608 :   void FFT<DataType_>::get_amp(std::vector<DataType_>& amp) const
     420             :   {
     421         608 :     amp.assign(size / 2 + 1, 0);
     422         608 :     impl->get_amp(amp);
     423         608 :   }
     424             :   
     425             :   template<class DataType_>
     426           0 :   void FFT<DataType_>::get_angle(std::vector<DataType_>& angle) const
     427             :   {
     428           0 :     angle.assign(size / 2 + 1, 0);
     429           0 :     impl->get_angle(angle);
     430           0 :   }
     431             : 
     432             :   template class FFT<double>;
     433             : }