common/alcomplex.cpp

   1
   2 #include "config.h"
   3
   4 #include "alcomplex.h"
   5
   6 #include <algorithm>
   7 #include <array>
   8 #include <cassert>
   9 #include <cstddef>
  10 #include <functional>
  11 #include <iterator>
  12 #include <utility>
  13
  14 #include "albit.h"
  15 #include "alnumbers.h"
  16 #include "alnumeric.h"
  17 #include "opthelpers.h"
  18
  19
  20 namespace {
  21
  22 using ushort = unsigned short;
  23 using ushort2 = std::array<ushort,2>;
  24 using complex_d = std::complex<double>;
  25
  26 constexpr std::size_t BitReverseCounter(std::size_t log2_size) noexcept
  27 {
  28     /* Some magic math that calculates the number of swaps needed for a
  29      * sequence of bit-reversed indices when index < reversed_index.
  30      */
  31     return (1_zu<<(log2_size-1)) - (1_zu<<((log2_size-1_zu)/2_zu));
  32 }
  33
  34
  35 template<std::size_t N>
  36 struct BitReverser {
  37     static_assert(N <= sizeof(ushort)*8, "Too many bits for the bit-reversal table.");
  38
  39     std::array<ushort2,BitReverseCounter(N)> mData{};
  40
  41     constexpr BitReverser()
  42     {
  43         const std::size_t fftsize{1u << N};
  44         std::size_t ret_i{0};
  45
  46         /* Bit-reversal permutation applied to a sequence of fftsize items. */
  47         for(std::size_t idx{1u};idx < fftsize-1;++idx)
  48         {
  49             std::size_t revidx{idx};
  50             revidx = ((revidx&0xaaaaaaaa) >> 1) | ((revidx&0x55555555) << 1);
  51             revidx = ((revidx&0xcccccccc) >> 2) | ((revidx&0x33333333) << 2);
  52             revidx = ((revidx&0xf0f0f0f0) >> 4) | ((revidx&0x0f0f0f0f) << 4);
  53             revidx = ((revidx&0xff00ff00) >> 8) | ((revidx&0x00ff00ff) << 8);
  54             revidx = (revidx >> 16) | ((revidx&0x0000ffff) << 16);
  55             revidx >>= 32-N;
  56
  57             if(idx < revidx)
  58             {
  59                 mData[ret_i][0]  = static_cast<ushort>(idx);
  60                 mData[ret_i][1] = static_cast<ushort>(revidx);
  61                 ++ret_i;
  62             }
  63         }
  64         assert(ret_i == std::size(mData));
  65     }
  66 };
  67
  68 /* These bit-reversal swap tables support up to 11-bit indices (2048 elements),
  69  * which is large enough for the filters and effects in OpenAL Soft. Larger FFT
  70  * requests will use a slower table-less path.
  71  */
  72 constexpr BitReverser<2> BitReverser2{};
  73 constexpr BitReverser<3> BitReverser3{};
  74 constexpr BitReverser<4> BitReverser4{};
  75 constexpr BitReverser<5> BitReverser5{};
  76 constexpr BitReverser<6> BitReverser6{};
  77 constexpr BitReverser<7> BitReverser7{};
  78 constexpr BitReverser<8> BitReverser8{};
  79 constexpr BitReverser<9> BitReverser9{};
  80 constexpr BitReverser<10> BitReverser10{};
  81 constexpr BitReverser<11> BitReverser11{};
  82 constexpr std::array<al::span<const ushort2>,12> gBitReverses{{
  83     {}, {},
  84     BitReverser2.mData,
  85     BitReverser3.mData,
  86     BitReverser4.mData,
  87     BitReverser5.mData,
  88     BitReverser6.mData,
  89     BitReverser7.mData,
  90     BitReverser8.mData,
  91     BitReverser9.mData,
  92     BitReverser10.mData,
  93     BitReverser11.mData
  94 }};
  95
  96 /* Lookup table for std::polar(1, pi / (1<<index)); */
  97 template<typename T>
  98 constexpr std::array<std::complex<T>,gBitReverses.size()-1> gArgAngle{{
  99     {static_cast<T>(-1.00000000000000000e+00), static_cast<T>(0.00000000000000000e+00)},
 100     {static_cast<T>( 0.00000000000000000e+00), static_cast<T>(1.00000000000000000e+00)},
 101     {static_cast<T>( 7.07106781186547524e-01), static_cast<T>(7.07106781186547524e-01)},
 102     {static_cast<T>( 9.23879532511286756e-01), static_cast<T>(3.82683432365089772e-01)},
 103     {static_cast<T>( 9.80785280403230449e-01), static_cast<T>(1.95090322016128268e-01)},
 104     {static_cast<T>( 9.95184726672196886e-01), static_cast<T>(9.80171403295606020e-02)},
 105     {static_cast<T>( 9.98795456205172393e-01), static_cast<T>(4.90676743274180143e-02)},
 106     {static_cast<T>( 9.99698818696204220e-01), static_cast<T>(2.45412285229122880e-02)},
 107     {static_cast<T>( 9.99924701839144541e-01), static_cast<T>(1.22715382857199261e-02)},
 108     {static_cast<T>( 9.99981175282601143e-01), static_cast<T>(6.13588464915447536e-03)},
 109     {static_cast<T>( 9.99995293809576172e-01), static_cast<T>(3.06795676296597627e-03)}
 110 }};
 111
 112 } // namespace
 113
 114 void complex_fft(const al::span<std::complex<double>> buffer, const double sign)
 115 {
 116     const std::size_t fftsize{buffer.size()};
 117     /* Get the number of bits used for indexing. Simplifies bit-reversal and
 118      * the main loop count.
 119      */
 120     const std::size_t log2_size{static_cast<std::size_t>(al::countr_zero(fftsize))};
 121
 122     if(log2_size < gBitReverses.size()) LIKELY
 123     {
 124         for(auto &rev : gBitReverses[log2_size])
 125             std::swap(buffer[rev[0]], buffer[rev[1]]);
 126
 127         /* Iterative form of Danielson-Lanczos lemma */
 128         for(std::size_t i{0};i < log2_size;++i)
 129         {
 130             const std::size_t step2{1_uz << i};
 131             const std::size_t step{2_uz << i};
 132             /* The first iteration of the inner loop would have u=1, which we
 133              * can simplify to remove a number of complex multiplies.
 134              */
 135             for(std::size_t k{0};k < fftsize;k+=step)
 136             {
 137                 const complex_d temp{buffer[k+step2]};
 138                 buffer[k+step2] = buffer[k] - temp;
 139                 buffer[k] += temp;
 140             }
 141
 142             const complex_d w{gArgAngle<double>[i].real(), gArgAngle<double>[i].imag()*sign};
 143             complex_d u{w};
 144             for(std::size_t j{1};j < step2;j++)
 145             {
 146                 for(std::size_t k{j};k < fftsize;k+=step)
 147                 {
 148                     const complex_d temp{buffer[k+step2] * u};
 149                     buffer[k+step2] = buffer[k] - temp;
 150                     buffer[k] += temp;
 151                 }
 152                 u *= w;
 153             }
 154         }
 155     }
 156     else
 157     {
 158         assert(log2_size < 32);
 159
 160         for(std::size_t idx{1u};idx < fftsize-1;++idx)
 161         {
 162             std::size_t revidx{idx};
 163             revidx = ((revidx&0xaaaaaaaa) >> 1) | ((revidx&0x55555555) << 1);
 164             revidx = ((revidx&0xcccccccc) >> 2) | ((revidx&0x33333333) << 2);
 165             revidx = ((revidx&0xf0f0f0f0) >> 4) | ((revidx&0x0f0f0f0f) << 4);
 166             revidx = ((revidx&0xff00ff00) >> 8) | ((revidx&0x00ff00ff) << 8);
 167             revidx = (revidx >> 16) | ((revidx&0x0000ffff) << 16);
 168             revidx >>= 32-log2_size;
 169
 170             if(idx < revidx)
 171                 std::swap(buffer[idx], buffer[revidx]);
 172         }
 173
 174         const double pi{al::numbers::pi * sign};
 175         for(std::size_t i{0};i < log2_size;++i)
 176         {
 177             const std::size_t step2{1_uz << i};
 178             const std::size_t step{2_uz << i};
 179             for(std::size_t k{0};k < fftsize;k+=step)
 180             {
 181                 const complex_d temp{buffer[k+step2]};
 182                 buffer[k+step2] = buffer[k] - temp;
 183                 buffer[k] += temp;
 184             }
 185
 186             const double arg{pi / static_cast<double>(step2)};
 187             const complex_d w{std::polar(1.0, arg)};
 188             complex_d u{w};
 189             for(std::size_t j{1};j < step2;j++)
 190             {
 191                 for(std::size_t k{j};k < fftsize;k+=step)
 192                 {
 193                     const complex_d temp{buffer[k+step2] * u};
 194                     buffer[k+step2] = buffer[k] - temp;
 195                     buffer[k] += temp;
 196                 }
 197                 u *= w;
 198             }
 199         }
 200     }
 201 }
 202
 203 void complex_hilbert(const al::span<std::complex<double>> buffer)
 204 {
 205     inverse_fft(buffer);
 206
 207     const double inverse_size = 1.0/static_cast<double>(buffer.size());
 208     auto bufiter = buffer.begin();
 209     const auto halfiter = bufiter + ptrdiff_t(buffer.size()>>1);
 210
 211     *bufiter *= inverse_size; ++bufiter;
 212     bufiter = std::transform(bufiter, halfiter, bufiter,
 213         [scale=inverse_size*2.0](std::complex<double> d){ return d * scale; });
 214     *bufiter *= inverse_size; ++bufiter;
 215
 216     std::fill(bufiter, buffer.end(), std::complex<double>{});
 217
 218     forward_fft(buffer);
 219 }