|
|
-
- #include "config.h"
-
- #include "alcomplex.h"
-
- #include <cmath>
-
- namespace {
-
- constexpr double Pi{3.141592653589793238462643383279502884};
-
- } // namespace
-
- void complex_fft(std::complex<double> *FFTBuffer, int FFTSize, double Sign)
- {
- /* Bit-reversal permutation applied to a sequence of FFTSize items */
- for(int i{1};i < FFTSize-1;i++)
- {
- int j{0};
- for(int mask{1};mask < FFTSize;mask <<= 1)
- {
- if((i&mask) != 0)
- j++;
- j <<= 1;
- }
- j >>= 1;
-
- if(i < j)
- std::swap(FFTBuffer[i], FFTBuffer[j]);
- }
-
- /* Iterative form of DanielsonLanczos lemma */
- int step{2};
- for(int i{1};i < FFTSize;i<<=1, step<<=1)
- {
- int 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(int j{0};j < step2;j++)
- {
- for(int k{j};k < FFTSize;k+=step)
- {
- std::complex<double> temp{FFTBuffer[k+step2] * u};
- FFTBuffer[k+step2] = FFTBuffer[k] - temp;
- FFTBuffer[k] += temp;
- }
-
- u *= w;
- }
- }
- }
-
- void complex_hilbert(std::complex<double> *Buffer, int size)
- {
- const double inverse_size = 1.0/static_cast<double>(size);
-
- for(int i{0};i < size;i++)
- Buffer[i].imag(0.0);
-
- complex_fft(Buffer, size, 1.0);
-
- int todo{size>>1};
- int i{0};
-
- Buffer[i++] *= inverse_size;
- while(i < todo)
- Buffer[i++] *= 2.0*inverse_size;
- Buffer[i++] *= inverse_size;
-
- for(;i < size;i++)
- Buffer[i] = std::complex<double>{};
-
- complex_fft(Buffer, size, -1.0);
- }
|