Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <vector>
- #include <cassert>
- #include <complex>
- #include <cmath>
- #include "fftw3.h"
- #ifndef M_PI
- #define M_PI 3.14159265358979323846
- #endif
- std::complex<double> compute_twiddle(size_t i, size_t len, bool inverse) {
- double constant;
- if(inverse) {
- constant = 2 * M_PI / len;
- } else {
- constant = -2 * M_PI / len;
- }
- double real = std::cos(i * constant);
- double imag = std::sin(i * constant);
- return std::complex<double>(real, imag);
- }
- std::vector<double> compute_dct4(std::vector<double> x) {
- assert(x.size()%2 == 0);
- size_t N = x.size();
- //pre-process the input by splitting into into two arrays, "w" and "v"
- std::vector<double> w(N/2), v(N/2);
- w[0] = x[0];
- v[N/2 - 1] = x[N - 1];
- for(size_t k = 1; k < N/2; k++) {
- w[k] = x[2 * k - 1] + x[2 * k];
- v[k - 1] = x[2 * k - 1] - x[2 * k];
- //NOTE: in the paper, the array "v" is 1-indexed rather than 0-indexed, so we have to compensate for this by subtracting 1 from every index
- //IE, the loop starts with k = 1, but since v is 1-indexed, this corresponds to v[0] in a 0-indexed array
- }
- //using FFTW, run a DCT-3 on array "w", and a DST-3 on array "v"
- std::vector<double> W(N/2), V(N/2);
- auto dct3_plan = fftw_plan_r2r_1d(N/2, w.data(), W.data(), FFTW_REDFT01, FFTW_ESTIMATE);
- auto dst3_plan = fftw_plan_r2r_1d(N/2, v.data(), V.data(), FFTW_RODFT01, FFTW_ESTIMATE);
- fftw_execute(dct3_plan);
- fftw_execute(dst3_plan);
- //post-process the data by combining it back into a single array
- std::vector<std::complex<double>> Z(N*2);
- for(size_t k = 0; k < N/2; k++) {
- Z[2 * k + 1] = std::complex<double>(W[k], V[k]);
- Z[2 * (x.size() - 1 - k) + 1] = std::complex<double>(W[k], -V[k]);
- }
- std::vector<double> final_result(N);
- for(size_t k = 0; k < N; k++) {
- std::complex<double> twiddle = compute_twiddle(2 * k + 1, 8*N, false);
- final_result[k] = (Z[2 * k + 1] * twiddle).real();
- }
- return final_result;
- }
- void print_vec(std::vector<double> vec) {
- if(vec.size() > 0) {
- std::cout << vec[0];
- }
- for(size_t i = 1; i < vec.size(); i++) {
- std::cout << ", " << vec[i];
- }
- std::cout << std::endl;
- }
- int main()
- {
- std::vector<double> input{-1.847365, 4.444461, -18.9616, 20.626879};
- std::vector<double> actual = compute_dct4(input);
- std::vector<double> expected(input.size());
- auto expected_plan = fftw_plan_r2r_1d(input.size(), input.data(), expected.data(), FFTW_REDFT11, FFTW_ESTIMATE);
- fftw_execute(expected_plan);
- std::cout << "Input: ";
- print_vec(input);
- std::cout << "Expected: ";
- print_vec(expected);
- std::cout << "Actual: ";
- print_vec(actual);
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement