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 <functional>
  11 #include <utility>
  12
  13 #include "albit.h"
  14 #include "alnumbers.h"
  15 #include "alnumeric.h"
  16 #include "opthelpers.h"
  17
  18
  19 namespace {
  20
  21 using ushort = unsigned short;
  22 using ushort2 = std::pair<ushort,ushort>;
  23
  24 constexpr size_t BitReverseCounter(size_t log2_size) noexcept
  25 {
  26     /* Some magic math that calculates the number of swaps needed for a
  27      * sequence of bit-reversed indices when index < reversed_index.
  28      */
  29     return (1u<<(log2_size-1)) - (1u<<((log2_size-1u)/2u));
  30 }
  31
  32
  33 template<size_t N>
  34 struct BitReverser {
  35     static_assert(N <= sizeof(ushort)*8, "Too many bits for the bit-reversal table.");
  36
  37     ushort2 mData[BitReverseCounter(N)]{};
  38
  39     constexpr BitReverser()
  40     {
  41         const size_t fftsize{1u << N};
  42         size_t ret_i{0};
  43
  44         /* Bit-reversal permutation applied to a sequence of fftsize items. */
  45         for(size_t idx{1u};idx < fftsize-1;++idx)
  46         {
  47             size_t revidx{0u}, imask{idx};
  48             for(size_t i{0};i < N;++i)
  49             {
  50                 revidx = (revidx<<1) | (imask&1);
  51                 imask >>= 1;
  52             }
  53
  54             if(idx < revidx)
  55             {
  56                 mData[ret_i].first  = static_cast<ushort>(idx);
  57                 mData[ret_i].second = static_cast<ushort>(revidx);
  58                 ++ret_i;
  59             }
  60         }
  61         assert(ret_i == al::size(mData));
  62     }
  63 };
  64
  65 /* These bit-reversal swap tables support up to 10-bit indices (1024 elements),
  66  * which is the largest used by OpenAL Soft's filters and effects. Larger FFT
  67  * requests, used by some utilities where performance is less important, will
  68  * use a slower table-less path.
  69  */
  70 constexpr BitReverser<2> BitReverser2{};
  71 constexpr BitReverser<3> BitReverser3{};
  72 constexpr BitReverser<4> BitReverser4{};
  73 constexpr BitReverser<5> BitReverser5{};
  74 constexpr BitReverser<6> BitReverser6{};
  75 constexpr BitReverser<7> BitReverser7{};
  76 constexpr BitReverser<8> BitReverser8{};
  77 constexpr BitReverser<9> BitReverser9{};
  78 constexpr BitReverser<10> BitReverser10{};
  79 constexpr std::array<al::span<const ushort2>,11> gBitReverses{{
  80     {}, {},
  81     BitReverser2.mData,
  82     BitReverser3.mData,
  83     BitReverser4.mData,
  84     BitReverser5.mData,
  85     BitReverser6.mData,
  86     BitReverser7.mData,
  87     BitReverser8.mData,
  88     BitReverser9.mData,
  89     BitReverser10.mData
  90 }};
  91
  92 } // namespace
  93
  94 template<typename Real>
  95 std::enable_if_t<std::is_floating_point<Real>::value>
  96 complex_fft(const al::span<std::complex<Real>> buffer, const al::type_identity_t<Real> 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(log2_size >= gBitReverses.size()) UNLIKELY
 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 Real pi{al::numbers::pi_v<Real> * sign};
 124     size_t step2{1u};
 125     for(size_t i{0};i < log2_size;++i)
 126     {
 127         const Real arg{pi / static_cast<Real>(step2)};
 128
 129         /* TODO: Would std::polar(1.0, arg) be any better? */
 130         const std::complex<Real> w{std::cos(arg), std::sin(arg)};
 131         std::complex<Real> 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<Real> 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     using namespace std::placeholders;
 152
 153     inverse_fft(buffer);
 154
 155     const double inverse_size = 1.0/static_cast<double>(buffer.size());
 156     auto bufiter = buffer.begin();
 157     const auto halfiter = bufiter + (buffer.size()>>1);
 158
 159     *bufiter *= inverse_size; ++bufiter;
 160     bufiter = std::transform(bufiter, halfiter, bufiter,
 161         [scale=inverse_size*2.0](std::complex<double> d){ return d * scale; });
 162     *bufiter *= inverse_size; ++bufiter;
 163
 164     std::fill(bufiter, buffer.end(), std::complex<double>{});
 165
 166     forward_fft(buffer);
 167 }
 168
 169
 170 template void complex_fft<>(const al::span<std::complex<float>> buffer, const float sign);
 171 template void complex_fft<>(const al::span<std::complex<double>> buffer, const double sign);