util/fft.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #include "terraingen/util/fft.h" | ||
| 2 | |||
| 3 | #include <cmath> | ||
| 4 | #include <cstddef> | ||
| 5 | #include <numbers> | ||
| 6 | #include <utility> | ||
| 7 | |||
| 8 | namespace terraingen { | ||
| 9 | |||
| 10 | namespace { | ||
| 11 | |||
| 12 | /// Reorders `values` by bit-reversed index, the input order of the iterative transform. | ||
| 13 | 57 | void BitReverse(std::vector<std::complex<double>>& values) { | |
| 14 | 57 | const std::size_t n = values.size(); | |
| 15 | 55329 | for (std::size_t i = 1, j = 0; i < n; ++i) { | |
| 16 | 55272 | std::size_t bit = n >> 1; | |
| 17 | 109996 | for (; (j & bit) != 0; bit >>= 1) { | |
| 18 | 54724 | j ^= bit; | |
| 19 | } | ||
| 20 | 55272 | j ^= bit; | |
| 21 | 55272 | if (i < j) { | |
| 22 | 26796 | std::swap(values[i], values[j]); | |
| 23 | } | ||
| 24 | } | ||
| 25 | 57 | } | |
| 26 | |||
| 27 | } // namespace | ||
| 28 | |||
| 29 | 57 | void Fft(std::vector<std::complex<double>>& values, double sign) { | |
| 30 | 57 | BitReverse(values); | |
| 31 | 57 | const std::size_t n = values.size(); | |
| 32 | 605 | for (std::size_t length = 2; length <= n; length <<= 1) { | |
| 33 | 548 | const double angle = sign * 2.0 * std::numbers::pi / static_cast<double>(length); | |
| 34 | 548 | const std::complex<double> step(std::cos(angle), std::sin(angle)); | |
| 35 | 55820 | for (std::size_t start = 0; start < n; start += length) { | |
| 36 | 55272 | std::complex<double> w(1.0, 0.0); | |
| 37 | 331816 | for (std::size_t k = 0; k < length / 2; ++k) { | |
| 38 | 276544 | const std::complex<double> even = values[start + k]; | |
| 39 | 276544 | const std::complex<double> odd = values[start + k + length / 2] * w; | |
| 40 | 276544 | values[start + k] = even + odd; | |
| 41 | 276544 | values[start + k + length / 2] = even - odd; | |
| 42 | 276544 | w *= step; | |
| 43 | } | ||
| 44 | } | ||
| 45 | } | ||
| 46 | 57 | } | |
| 47 | |||
| 48 | } // namespace terraingen | ||
| 49 |