diff options
author | Chris Robinson <[email protected]> | 2018-11-17 04:14:57 -0800 |
---|---|---|
committer | Chris Robinson <[email protected]> | 2018-11-17 04:14:57 -0800 |
commit | ccdaca80c910047e16f710d44f640a6d6f86a195 (patch) | |
tree | 0cf95ed7a5fd6478c01cea8bd8e9c521478ccf0d /common/alcomplex.cpp | |
parent | 09943683b5872943cd1f9211ef2a77922906b906 (diff) |
Use standard complex types instead of custom
Diffstat (limited to 'common/alcomplex.cpp')
-rw-r--r-- | common/alcomplex.cpp | 76 |
1 files changed, 76 insertions, 0 deletions
diff --git a/common/alcomplex.cpp b/common/alcomplex.cpp new file mode 100644 index 00000000..de2d1be9 --- /dev/null +++ b/common/alcomplex.cpp @@ -0,0 +1,76 @@ + +#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 Danielson–Lanczos 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/(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); +} |