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 imask{i + fftsize};imask;imask >>= 1)
  22             j = (j<<1) + (imask&1);
  23         j >>= 1;
  24
  25         if(i < j)
  26             std::swap(buffer[i], buffer[j]);
  27     }
  28
  29     /* Iterative form of Danielson-Lanczos lemma */
  30     size_t step{2u};
  31     for(size_t i{1u};i < fftsize;i<<=1, step<<=1)
  32     {
  33         const size_t step2{step >> 1};
  34         const double arg{al::MathDefs<double>::Pi() / static_cast<double>(step2)};
  35
  36         const std::complex<double> w{std::cos(arg), std::sin(arg)*sign};
  37         std::complex<double> u{1.0, 0.0};
  38         for(size_t j{0};j < step2;j++)
  39         {
  40             for(size_t k{j};k < fftsize;k+=step)
  41             {
  42                 std::complex<double> temp{buffer[k+step2] * u};
  43                 buffer[k+step2] = buffer[k] - temp;
  44                 buffer[k] += temp;
  45             }
  46
  47             u *= w;
  48         }
  49     }
  50 }
  51
  52 void complex_hilbert(const al::span<std::complex<double>> buffer)
  53 {
  54     inverse_fft(buffer);
  55
  56     const double inverse_size = 1.0/static_cast<double>(buffer.size());
  57     auto bufiter = buffer.begin();
  58     const auto halfiter = bufiter + (buffer.size()>>1);
  59
  60     *bufiter *= inverse_size; ++bufiter;
  61     bufiter = std::transform(bufiter, halfiter, bufiter,
  62         [inverse_size](const std::complex<double> &c) -> std::complex<double>
  63         { return c * (2.0*inverse_size); });
  64     *bufiter *= inverse_size; ++bufiter;
  65
  66     std::fill(bufiter, buffer.end(), std::complex<double>{});
  67
  68     forward_fft(buffer);
  69 }