14 #include "alnumbers.h"
15 #include "alnumeric.h"
16 #include "opthelpers.h"
21 using ushort = unsigned short;
22 using ushort2 = std::pair<ushort,ushort>;
24 constexpr size_t BitReverseCounter(size_t log2_size) noexcept
26 /* Some magic math that calculates the number of swaps needed for a
27 * sequence of bit-reversed indices when index < reversed_index.
29 return (1u<<(log2_size-1)) - (1u<<((log2_size-1u)/2u));
35 static_assert(N <= sizeof(ushort)*8, "Too many bits for the bit-reversal table.");
37 ushort2 mData[BitReverseCounter(N)]{};
39 constexpr BitReverser()
41 const size_t fftsize{1u << N};
44 /* Bit-reversal permutation applied to a sequence of fftsize items. */
45 for(size_t idx{1u};idx < fftsize-1;++idx)
47 size_t revidx{0u}, imask{idx};
48 for(size_t i{0};i < N;++i)
50 revidx = (revidx<<1) | (imask&1);
56 mData[ret_i].first = static_cast<ushort>(idx);
57 mData[ret_i].second = static_cast<ushort>(revidx);
61 assert(ret_i == al::size(mData));
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.
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{{
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)
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.
102 const size_t log2_size{static_cast<size_t>(al::countr_zero(fftsize))};
104 if(log2_size >= gBitReverses.size()) UNLIKELY
106 for(size_t idx{1u};idx < fftsize-1;++idx)
108 size_t revidx{0u}, imask{idx};
109 for(size_t i{0};i < log2_size;++i)
111 revidx = (revidx<<1) | (imask&1);
116 std::swap(buffer[idx], buffer[revidx]);
119 else for(auto &rev : gBitReverses[log2_size])
120 std::swap(buffer[rev.first], buffer[rev.second]);
122 /* Iterative form of Danielson-Lanczos lemma */
123 const Real pi{al::numbers::pi_v<Real> * sign};
125 for(size_t i{0};i < log2_size;++i)
127 const Real arg{pi / static_cast<Real>(step2)};
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++)
135 for(size_t k{j};k < fftsize;k+=step)
137 std::complex<Real> temp{buffer[k+step2] * u};
138 buffer[k+step2] = buffer[k] - temp;
149 void complex_hilbert(const al::span<std::complex<double>> buffer)
151 using namespace std::placeholders;
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);
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;
164 std::fill(bufiter, buffer.end(), std::complex<double>{});
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);