common/alcomplex.cpp

   1
   2 #include "config.h"
   3
   4 #include "alcomplex.h"
   5
   6 #include <algorithm>
   7 #include <cmath>
   8 #include <cstddef>
   9 #include <utility>
  10
  11 #include "albit.h"
  12 #include "alnumeric.h"
  13 #include "math_defs.h"
  14
  15
  16 void complex_fft(const al::span<std::complex<double>> buffer, const double sign)
  17 {
  18     const size_t fftsize{buffer.size()};
  19     /* Get the number of bits used for indexing. Simplifies bit-reversal and
  20      * the main loop count.
  21      */
  22     const size_t log2_size{static_cast<size_t>(al::countr_zero(fftsize))};
  23
  24     /* Bit-reversal permutation applied to a sequence of fftsize items. */
  25     for(size_t idx{1u};idx < fftsize-1;++idx)
  26     {
  27         size_t revidx{0u}, imask{idx};
  28         for(size_t i{0};i < log2_size;++i)
  29         {
  30             revidx = (revidx<<1) | (imask&1);
  31             imask >>= 1;
  32         }
  33
  34         if(idx < revidx)
  35             std::swap(buffer[idx], buffer[revidx]);
  36     }
  37
  38     /* Iterative form of Danielson-Lanczos lemma */
  39     size_t step2{1u};
  40     for(size_t i{0};i < log2_size;++i)
  41     {
  42         const double arg{al::MathDefs<double>::Pi() / static_cast<double>(step2)};
  43
  44         const std::complex<double> w{std::cos(arg), std::sin(arg)*sign};
  45         std::complex<double> u{1.0, 0.0};
  46         const size_t step{step2 << 1};
  47         for(size_t j{0};j < step2;j++)
  48         {
  49             for(size_t k{j};k < fftsize;k+=step)
  50             {
  51                 std::complex<double> temp{buffer[k+step2] * u};
  52                 buffer[k+step2] = buffer[k] - temp;
  53                 buffer[k] += temp;
  54             }
  55
  56             u *= w;
  57         }
  58
  59         step2 <<= 1;
  60     }
  61 }
  62
  63 void complex_hilbert(const al::span<std::complex<double>> buffer)
  64 {
  65     inverse_fft(buffer);
  66
  67     const double inverse_size = 1.0/static_cast<double>(buffer.size());
  68     auto bufiter = buffer.begin();
  69     const auto halfiter = bufiter + (buffer.size()>>1);
  70
  71     *bufiter *= inverse_size; ++bufiter;
  72     bufiter = std::transform(bufiter, halfiter, bufiter,
  73         [inverse_size](const std::complex<double> &c) -> std::complex<double>
  74         { return c * (2.0*inverse_size); });
  75     *bufiter *= inverse_size; ++bufiter;
  76
  77     std::fill(bufiter, buffer.end(), std::complex<double>{});
  78
  79     forward_fft(buffer);
  80 }