blob: de10ede20f838aea6eedbd1716ee5e9a2f41bf9b (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
|
#include "config.h"
#include "alcomplex.h"
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <utility>
#include "albit.h"
#include "alnumeric.h"
#include "math_defs.h"
void complex_fft(const al::span<std::complex<double>> buffer, const double sign)
{
const size_t fftsize{buffer.size()};
/* Get the number of bits used for indexing. Simplifies bit-reversal and
* the main loop count.
*/
const size_t log2_size{static_cast<size_t>(al::countr_zero(fftsize))};
/* Bit-reversal permutation applied to a sequence of fftsize items. */
for(size_t idx{1u};idx < fftsize-1;++idx)
{
size_t revidx{0u}, imask{idx};
for(size_t i{0};i < log2_size;++i)
{
revidx = (revidx<<1) | (imask&1);
imask >>= 1;
}
if(idx < revidx)
std::swap(buffer[idx], buffer[revidx]);
}
/* Iterative form of Danielson-Lanczos lemma */
size_t step2{1u};
for(size_t i{0};i < log2_size;++i)
{
const double arg{al::MathDefs<double>::Pi() / static_cast<double>(step2)};
const std::complex<double> w{std::cos(arg), std::sin(arg)*sign};
std::complex<double> u{1.0, 0.0};
const size_t step{step2 << 1};
for(size_t j{0};j < step2;j++)
{
for(size_t k{j};k < fftsize;k+=step)
{
std::complex<double> temp{buffer[k+step2] * u};
buffer[k+step2] = buffer[k] - temp;
buffer[k] += temp;
}
u *= w;
}
step2 <<= 1;
}
}
void complex_hilbert(const al::span<std::complex<double>> buffer)
{
inverse_fft(buffer);
const double inverse_size = 1.0/static_cast<double>(buffer.size());
auto bufiter = buffer.begin();
const auto halfiter = bufiter + (buffer.size()>>1);
*bufiter *= inverse_size; ++bufiter;
bufiter = std::transform(bufiter, halfiter, bufiter,
[inverse_size](const std::complex<double> &c) -> std::complex<double>
{ return c * (2.0*inverse_size); });
*bufiter *= inverse_size; ++bufiter;
std::fill(bufiter, buffer.end(), std::complex<double>{});
forward_fft(buffer);
}
|
