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 "math_defs.h"
  12
  13
  14 void complex_fft(const al::span<std::complex<double>> buffer, const double sign)
  15 {
  16     const size_t fftsize{buffer.size()};
  17     /* Bit-reversal permutation applied to a sequence of FFTSize items */
  18     for(size_t i{1u};i < fftsize-1;i++)
  19     {
  20         size_t j{0u};
  21         for(size_t mask{1u};mask < fftsize;mask <<= 1)
  22         {
  23             if((i&mask) != 0)
  24                 j++;
  25             j <<= 1;
  26         }
  27         j >>= 1;
  28
  29         if(i < j)
  30             std::swap(buffer[i], buffer[j]);
  31     }
  32
  33     /* Iterative form of DanielsonLanczos lemma */
  34     size_t step{2u};
  35     for(size_t i{1u};i < fftsize;i<<=1, step<<=1)
  36     {
  37         const size_t step2{step >> 1};
  38         double arg{al::MathDefs<double>::Pi() / static_cast<double>(step2)};
  39
  40         std::complex<double> w{std::cos(arg), std::sin(arg)*sign};
  41         std::complex<double> u{1.0, 0.0};
  42         for(size_t j{0};j < step2;j++)
  43         {
  44             for(size_t k{j};k < fftsize;k+=step)
  45             {
  46                 std::complex<double> temp{buffer[k+step2] * u};
  47                 buffer[k+step2] = buffer[k] - temp;
  48                 buffer[k] += temp;
  49             }
  50
  51             u *= w;
  52         }
  53     }
  54 }
  55
  56 void complex_hilbert(const al::span<std::complex<double>> buffer)
  57 {
  58     complex_fft(buffer, 1.0);
  59
  60     const double inverse_size = 1.0/static_cast<double>(buffer.size());
  61     auto bufiter = buffer.begin();
  62     const auto halfiter = bufiter + (buffer.size()>>1);
  63
  64     *bufiter *= inverse_size; ++bufiter;
  65     bufiter = std::transform(bufiter, halfiter, bufiter,
  66         [inverse_size](const std::complex<double> &c) -> std::complex<double>
  67         { return c * (2.0*inverse_size); });
  68     *bufiter *= inverse_size; ++bufiter;
  69
  70     std::fill(bufiter, buffer.end(), std::complex<double>{});
  71
  72     complex_fft(buffer, -1.0);
  73 }