15 #include "alnumbers.h"
16 #include "alnumeric.h"
17 #include "opthelpers.h"
22 using ushort
= unsigned short;
23 using ushort2
= std::pair
<ushort
,ushort
>;
24 using complex_d
= std::complex<double>;
26 constexpr std::size_t BitReverseCounter(std::size_t log2_size
) noexcept
28 /* Some magic math that calculates the number of swaps needed for a
29 * sequence of bit-reversed indices when index < reversed_index.
31 return (1_zu
<<(log2_size
-1)) - (1_zu
<<((log2_size
-1_zu
)/2_zu
));
35 template<std::size_t N
>
37 static_assert(N
<= sizeof(ushort
)*8, "Too many bits for the bit-reversal table.");
39 std::array
<ushort2
,BitReverseCounter(N
)> mData
{};
41 constexpr BitReverser()
43 const std::size_t fftsize
{1u << N
};
46 /* Bit-reversal permutation applied to a sequence of fftsize items. */
47 for(std::size_t idx
{1u};idx
< fftsize
-1;++idx
)
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);
59 mData
[ret_i
].first
= static_cast<ushort
>(idx
);
60 mData
[ret_i
].second
= static_cast<ushort
>(revidx
);
64 assert(ret_i
== std::size(mData
));
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.
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
{{
96 /* Lookup table for std::polar(1, pi / (1<<index)); */
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)}
114 void complex_fft(const al::span
<std::complex<double>> buffer
, const double sign
)
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.
120 const std::size_t log2_size
{static_cast<std::size_t>(al::countr_zero(fftsize
))};
122 if(log2_size
< gBitReverses
.size()) LIKELY
124 for(auto &rev
: gBitReverses
[log2_size
])
125 std::swap(buffer
[rev
.first
], buffer
[rev
.second
]);
127 /* Iterative form of Danielson-Lanczos lemma */
128 for(std::size_t i
{0};i
< log2_size
;++i
)
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.
135 for(std::size_t k
{0};k
< fftsize
;k
+=step
)
137 const complex_d temp
{buffer
[k
+step2
]};
138 buffer
[k
+step2
] = buffer
[k
] - temp
;
142 const complex_d w
{gArgAngle
<double>[i
].real(), gArgAngle
<double>[i
].imag()*sign
};
144 for(std::size_t j
{1};j
< step2
;j
++)
146 for(std::size_t k
{j
};k
< fftsize
;k
+=step
)
148 const complex_d temp
{buffer
[k
+step2
] * u
};
149 buffer
[k
+step2
] = buffer
[k
] - temp
;
158 assert(log2_size
< 32);
160 for(std::size_t idx
{1u};idx
< fftsize
-1;++idx
)
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
;
171 std::swap(buffer
[idx
], buffer
[revidx
]);
174 const double pi
{al::numbers::pi
* sign
};
175 for(std::size_t i
{0};i
< log2_size
;++i
)
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
)
181 const complex_d temp
{buffer
[k
+step2
]};
182 buffer
[k
+step2
] = buffer
[k
] - temp
;
186 const double arg
{pi
/ static_cast<double>(step2
)};
187 const complex_d w
{std::polar(1.0, arg
)};
189 for(std::size_t j
{1};j
< step2
;j
++)
191 for(std::size_t k
{j
};k
< fftsize
;k
+=step
)
193 const complex_d temp
{buffer
[k
+step2
] * u
};
194 buffer
[k
+step2
] = buffer
[k
] - temp
;
203 void complex_hilbert(const al::span
<std::complex<double>> buffer
)
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);
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
;
216 std::fill(bufiter
, buffer
.end(), std::complex<double>{});