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
|
#include "config.h"
#include "alcomplex.h"
#include <cmath>
#include <algorithm>
namespace {
constexpr double Pi{3.141592653589793238462643383279502884};
} // namespace
void complex_fft(const al::span<std::complex<double>> buffer, const double sign)
{
const size_t fftsize{buffer.size()};
/* Bit-reversal permutation applied to a sequence of FFTSize items */
for(size_t i{1u};i < fftsize-1;i++)
{
size_t j{0u};
for(size_t mask{1u};mask < fftsize;mask <<= 1)
{
if((i&mask) != 0)
j++;
j <<= 1;
}
j >>= 1;
if(i < j)
std::swap(buffer[i], buffer[j]);
}
/* Iterative form of DanielsonLanczos lemma */
size_t step{2u};
for(size_t i{1u};i < fftsize;i<<=1, step<<=1)
{
const size_t step2{step >> 1};
double arg{Pi / step2};
std::complex<double> w{std::cos(arg), std::sin(arg)*sign};
std::complex<double> u{1.0, 0.0};
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;
}
}
}
void complex_hilbert(const al::span<std::complex<double>> buffer)
{
std::for_each(buffer.begin(), buffer.end(), [](std::complex<double> &c) { c.imag(0.0); });
complex_fft(buffer, 1.0);
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>{});
complex_fft(buffer, -1.0);
}
|
