Created
February 7, 2023 19:33
-
-
Save omnp/eb8d5551dd209cda26edc3e7a0361043 to your computer and use it in GitHub Desktop.
something like DFT (discrete fourier transform) in C++
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #pragma once | |
| #ifdef _WIN64 | |
| #define _USE_MATH_DEFINES | |
| #endif | |
| #include <cstdint> | |
| #include <complex> | |
| #include <cmath> | |
| #include <vector> | |
| #include <map> | |
| template <typename T> class dft | |
| { | |
| using complex = std::complex<T>; | |
| using vector = std::vector<complex>; | |
| public: | |
| dft(const uint32_t size) { | |
| uint32_t key = 1u; | |
| for (uint32_t k = 0; k < size+1; k++) { | |
| factors[key] = std::exp<T>(complex(0.0, -2.0) * complex(M_PI) / complex(key)); | |
| key <<= 1u; | |
| } | |
| } | |
| vector transform(vector const xs) { | |
| const uint32_t n = xs.size(); | |
| if (n <= 1u) | |
| return xs; | |
| const uint32_t n2 = n >> 1; | |
| auto xs1 = vector(xs.begin(), xs.begin() + n2); | |
| auto xs2 = vector(xs.begin() + n2, xs.end()); | |
| vector A = xs1; | |
| vector B = xs1; | |
| vector C(n); | |
| const complex factor = factors.at(n); | |
| #pragma omp for simd | |
| for (uint32_t i = 0; i < n2; i++) | |
| A[i] += xs2[i]; | |
| #pragma omp for simd | |
| for (uint32_t i = 0; i < n2; i++) | |
| B[i] -= xs2[i]; | |
| complex multiplier{1.0,0.0}; | |
| for (uint32_t i = 0; i < n2; i++) { | |
| B[i] *= multiplier; | |
| multiplier *= factor; | |
| } | |
| const vector & U = transform(A); | |
| const vector & V = transform(B); | |
| #pragma omp for simd | |
| for (uint32_t i = 0; i < n2; i++) | |
| C[2*i] = U[i]; | |
| #pragma omp for simd | |
| for (uint32_t i = 0; i < n2; i++) | |
| C[2*i+1] = V[i]; | |
| return C; | |
| } | |
| vector inverse(vector const & xs) { | |
| vector ys = xs; | |
| const complex n{ys.size()}; | |
| #pragma omp for simd | |
| for (uint32_t i = 0; i < xs.size(); i++) | |
| ys[i] = std::conj<T>(ys[i]) / n; | |
| return transform(ys); | |
| } | |
| private: | |
| std::map<uint32_t, complex> factors; | |
| }; |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment