Discrete Cosine Transform Type 4 implementation

//License: Choose GPLv2 or GPLv3

#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] * 2;
    v[N/2 - 1] = x[N - 1] * 2;

    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;
}