74 lines
1.9 KiB
C++
74 lines
1.9 KiB
C++
|
||
#include "config.h"
|
||
|
||
#include "alcomplex.h"
|
||
|
||
#include <algorithm>
|
||
#include <cmath>
|
||
#include <cstddef>
|
||
#include <utility>
|
||
|
||
#include "math_defs.h"
|
||
|
||
|
||
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{al::MathDefs<double>::Pi() / static_cast<double>(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)
|
||
{
|
||
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) - 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);
|
||
}
|