common/alcomplex.cpp

   1
   2 #include "config.h"
   3
   4 #include "alcomplex.h"
   5
   6 #include <cmath>
   7
   8 namespace {
   9
  10 constexpr double Pi{3.141592653589793238462643383279502884};
  11
  12 } // namespace
  13
  14 void complex_fft(std::complex<double> *FFTBuffer, int FFTSize, double Sign)
  15 {
  16     /* Bit-reversal permutation applied to a sequence of FFTSize items */
  17     for(int i{1};i < FFTSize-1;i++)
  18     {
  19         int j{0};
  20         for(int mask{1};mask < FFTSize;mask <<= 1)
  21         {
  22             if((i&mask) != 0)
  23                 j++;
  24             j <<= 1;
  25         }
  26         j >>= 1;
  27
  28         if(i < j)
  29             std::swap(FFTBuffer[i], FFTBuffer[j]);
  30     }
  31
  32     /* Iterative form of DanielsonLanczos lemma */
  33     int step{2};
  34     for(int i{1};i < FFTSize;i<<=1, step<<=1)
  35     {
  36         int step2{step >> 1};
  37         double arg{Pi / step2};
  38
  39         std::complex<double> w{std::cos(arg), std::sin(arg)*Sign};
  40         std::complex<double> u{1.0, 0.0};
  41         for(int j{0};j < step2;j++)
  42         {
  43             for(int k{j};k < FFTSize;k+=step)
  44             {
  45                 std::complex<double> temp{FFTBuffer[k+step2] * u};
  46                 FFTBuffer[k+step2] = FFTBuffer[k] - temp;
  47                 FFTBuffer[k] += temp;
  48             }
  49
  50             u *= w;
  51         }
  52     }
  53 }
  54
  55 void complex_hilbert(std::complex<double> *Buffer, int size)
  56 {
  57     const double inverse_size = 1.0/static_cast<double>(size);
  58
  59     for(int i{0};i < size;i++)
  60         Buffer[i].imag(0.0);
  61
  62     complex_fft(Buffer, size, 1.0);
  63
  64     int todo{size>>1};
  65     int i{0};
  66
  67     Buffer[i++] *= inverse_size;
  68     while(i < todo)
  69         Buffer[i++] *= 2.0*inverse_size;
  70     Buffer[i++] *= inverse_size;
  71
  72     for(;i < size;i++)
  73         Buffer[i] = std::complex<double>{};
  74
  75     complex_fft(Buffer, size, -1.0);
  76 }