🛠️🐜 Antkeeper superbuild with dependencies included https://antkeeper.com
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

76 lines
1.7 KiB

  1. #include "config.h"
  2. #include "alcomplex.h"
  3. #include <cmath>
  4. namespace {
  5. constexpr double Pi{3.141592653589793238462643383279502884};
  6. } // namespace
  7. void complex_fft(std::complex<double> *FFTBuffer, int FFTSize, double Sign)
  8. {
  9. /* Bit-reversal permutation applied to a sequence of FFTSize items */
  10. for(int i{1};i < FFTSize-1;i++)
  11. {
  12. int j{0};
  13. for(int mask{1};mask < FFTSize;mask <<= 1)
  14. {
  15. if((i&mask) != 0)
  16. j++;
  17. j <<= 1;
  18. }
  19. j >>= 1;
  20. if(i < j)
  21. std::swap(FFTBuffer[i], FFTBuffer[j]);
  22. }
  23. /* Iterative form of Danielson–Lanczos lemma */
  24. int step{2};
  25. for(int i{1};i < FFTSize;i<<=1, step<<=1)
  26. {
  27. int step2{step >> 1};
  28. double arg{Pi / step2};
  29. std::complex<double> w{std::cos(arg), std::sin(arg)*Sign};
  30. std::complex<double> u{1.0, 0.0};
  31. for(int j{0};j < step2;j++)
  32. {
  33. for(int k{j};k < FFTSize;k+=step)
  34. {
  35. std::complex<double> temp{FFTBuffer[k+step2] * u};
  36. FFTBuffer[k+step2] = FFTBuffer[k] - temp;
  37. FFTBuffer[k] += temp;
  38. }
  39. u *= w;
  40. }
  41. }
  42. }
  43. void complex_hilbert(std::complex<double> *Buffer, int size)
  44. {
  45. const double inverse_size = 1.0/static_cast<double>(size);
  46. for(int i{0};i < size;i++)
  47. Buffer[i].imag(0.0);
  48. complex_fft(Buffer, size, 1.0);
  49. int todo{size>>1};
  50. int i{0};
  51. Buffer[i++] *= inverse_size;
  52. while(i < todo)
  53. Buffer[i++] *= 2.0*inverse_size;
  54. Buffer[i++] *= inverse_size;
  55. for(;i < size;i++)
  56. Buffer[i] = std::complex<double>{};
  57. complex_fft(Buffer, size, -1.0);
  58. }