common/alcomplex.cpp

   1
   2 #include "config.h"
   3
   4 #include "alcomplex.h"
   5
   6 #include <algorithm>
   7 #include <cassert>
   8 #include <cmath>
   9 #include <cstddef>
  10 #include <utility>
  11
  12 #include "albit.h"
  13 #include "alnumbers.h"
  14 #include "alnumeric.h"
  15 #include "opthelpers.h"
  16
  17
  18 namespace {
  19
  20 using ushort = unsigned short;
  21 using ushort2 = std::pair<ushort,ushort>;
  22
  23 /* Because std::array doesn't have constexpr non-const accessors in C++14. */
  24 template<typename T, size_t N>
  25 struct our_array {
  26     T mData[N];
  27 };
  28
  29 constexpr size_t BitReverseCounter(size_t log2_size) noexcept
  30 {
  31     /* Some magic math that calculates the number of swaps needed for a
  32      * sequence of bit-reversed indices when index < reversed_index.
  33      */
  34     return (1u<<(log2_size-1)) - (1u<<((log2_size-1u)/2u));
  35 }
  36
  37 template<size_t N>
  38 constexpr auto GetBitReverser() noexcept
  39 {
  40     static_assert(N <= sizeof(ushort)*8, "Too many bits for the bit-reversal table.");
  41
  42     our_array<ushort2, BitReverseCounter(N)> ret{};
  43     const size_t fftsize{1u << N};
  44     size_t ret_i{0};
  45
  46     /* Bit-reversal permutation applied to a sequence of fftsize items. */
  47     for(size_t idx{1u};idx < fftsize-1;++idx)
  48     {
  49         size_t revidx{0u}, imask{idx};
  50         for(size_t i{0};i < N;++i)
  51         {
  52             revidx = (revidx<<1) | (imask&1);
  53             imask >>= 1;
  54         }
  55
  56         if(idx < revidx)
  57         {
  58             ret.mData[ret_i].first  = static_cast<ushort>(idx);
  59             ret.mData[ret_i].second = static_cast<ushort>(revidx);
  60             ++ret_i;
  61         }
  62     }
  63     assert(ret_i == al::size(ret.mData));
  64     return ret;
  65 }
  66
  67 /* These bit-reversal swap tables support up to 10-bit indices (1024 elements),
  68  * which is the largest used by OpenAL Soft's filters and effects. Larger FFT
  69  * requests, used by some utilities where performance is less important, will
  70  * use a slower table-less path.
  71  */
  72 constexpr auto BitReverser2 = GetBitReverser<2>();
  73 constexpr auto BitReverser3 = GetBitReverser<3>();
  74 constexpr auto BitReverser4 = GetBitReverser<4>();
  75 constexpr auto BitReverser5 = GetBitReverser<5>();
  76 constexpr auto BitReverser6 = GetBitReverser<6>();
  77 constexpr auto BitReverser7 = GetBitReverser<7>();
  78 constexpr auto BitReverser8 = GetBitReverser<8>();
  79 constexpr auto BitReverser9 = GetBitReverser<9>();
  80 constexpr auto BitReverser10 = GetBitReverser<10>();
  81 constexpr al::span<const ushort2> gBitReverses[11]{
  82     {}, {},
  83     BitReverser2.mData,
  84     BitReverser3.mData,
  85     BitReverser4.mData,
  86     BitReverser5.mData,
  87     BitReverser6.mData,
  88     BitReverser7.mData,
  89     BitReverser8.mData,
  90     BitReverser9.mData,
  91     BitReverser10.mData
  92 };
  93
  94 } // namespace
  95
  96 void complex_fft(const al::span<std::complex<double>> buffer, const double sign)
  97 {
  98     const size_t fftsize{buffer.size()};
  99     /* Get the number of bits used for indexing. Simplifies bit-reversal and
 100      * the main loop count.
 101      */
 102     const size_t log2_size{static_cast<size_t>(al::countr_zero(fftsize))};
 103
 104     if(unlikely(log2_size >= al::size(gBitReverses)))
 105     {
 106         for(size_t idx{1u};idx < fftsize-1;++idx)
 107         {
 108             size_t revidx{0u}, imask{idx};
 109             for(size_t i{0};i < log2_size;++i)
 110             {
 111                 revidx = (revidx<<1) | (imask&1);
 112                 imask >>= 1;
 113             }
 114
 115             if(idx < revidx)
 116                 std::swap(buffer[idx], buffer[revidx]);
 117         }
 118     }
 119     else for(auto &rev : gBitReverses[log2_size])
 120         std::swap(buffer[rev.first], buffer[rev.second]);
 121
 122     /* Iterative form of Danielson-Lanczos lemma */
 123     const double pi{al::numbers::pi * sign};
 124     size_t step2{1u};
 125     for(size_t i{0};i < log2_size;++i)
 126     {
 127         const double arg{pi / static_cast<double>(step2)};
 128
 129         /* TODO: Would std::polar(1.0, arg) be any better? */
 130         const std::complex<double> w{std::cos(arg), std::sin(arg)};
 131         std::complex<double> u{1.0, 0.0};
 132         const size_t step{step2 << 1};
 133         for(size_t j{0};j < step2;j++)
 134         {
 135             for(size_t k{j};k < fftsize;k+=step)
 136             {
 137                 std::complex<double> temp{buffer[k+step2] * u};
 138                 buffer[k+step2] = buffer[k] - temp;
 139                 buffer[k] += temp;
 140             }
 141
 142             u *= w;
 143         }
 144
 145         step2 <<= 1;
 146     }
 147 }
 148
 149 void complex_hilbert(const al::span<std::complex<double>> buffer)
 150 {
 151     inverse_fft(buffer);
 152
 153     const double inverse_size = 1.0/static_cast<double>(buffer.size());
 154     auto bufiter = buffer.begin();
 155     const auto halfiter = bufiter + (buffer.size()>>1);
 156
 157     *bufiter *= inverse_size; ++bufiter;
 158     bufiter = std::transform(bufiter, halfiter, bufiter,
 159         [inverse_size](const std::complex<double> &c) -> std::complex<double>
 160         { return c * (2.0*inverse_size); });
 161     *bufiter *= inverse_size; ++bufiter;
 162
 163     std::fill(bufiter, buffer.end(), std::complex<double>{});
 164
 165     forward_fft(buffer);
 166 }