libsimple_compress/unit_tests/dct.cpp

#include "simple/compress/fft.hpp"
#include "simple/support/random.hpp"

#include <cassert>
#include <vector>
#include <cmath>
#include <algorithm>
#include <numeric>
#include <complex>
#include <random>
#include <iostream>
#include <iomanip>

using namespace simple::compress;

const float tau = 2*std::acos(-1);

// DCT more widely used than DFT, as it has several practical advantages
// the end result is the same size as input, instead of double the size
// cosine being even function behaves better at endpoints of the window
template <typename Waveform>
void DCT_1(std::size_t size, Waveform waveform)
{

	std::vector<float> wave;
	wave.resize(size);

	std::iota(wave.begin(), wave.end(), 0);
	std::transform(wave.begin(), wave.end(), wave.begin(), [&](auto x) { return x/size; });
	std::transform(wave.begin(), wave.end(), wave.begin(), waveform);

	const auto dft_size = 2*(size-1);

	wave.resize(dft_size);
	// symmetric extension of the waveform results in canceling of sine/imaginary components
	std::copy(wave.rbegin()+size-1, wave.rend()-1, wave.begin() + size);

	std::vector<std::complex<float>> frequencies;
	frequencies.resize(dft_size);

	fft(std::multiplies<>{}, std::polar(1.f,tau/wave.size()), wave, frequencies);
	for(auto&& x : frequencies)
		x = x / float(dft_size);

	for(std::size_t i = 0; i != size; ++i)
		assert(std::abs(frequencies[i].imag()) < 0.001);

	for(auto&& x : frequencies)
		x.imag(0);

	// no idea how or why symmetric extension of the frequencies works to do the inverse... magic...
	std::copy(frequencies.rbegin()+size-1, frequencies.rend()-1, frequencies.begin() + size);

	std::vector<std::complex<float>> aaanback;
	aaanback.resize(dft_size);
	fft(std::multiplies<>{}, std::polar(1.f,tau/wave.size()), frequencies, aaanback);

	for(std::size_t i = 0; i != size; ++i)
		assert(std::abs(aaanback[i].real() - wave[i]) < 0.001);

	wave.resize(size);
	frequencies.resize(size);
	aaanback.resize(size);

#if defined SIMPLE_SUPPORT_DEBUG_HPP
	if(size <= 33)
	{
		std::cerr << std::fixed;
		std::cerr << std::setprecision(3);
		simple::support::print('\n');

		simple::support::print("IN : [ ");
		for(auto&& x : wave)
			simple::support::print(x, " ");
		simple::support::print("]\n");

		simple::support::print("DEC: [ ");
		for(auto&& x : aaanback)
			simple::support::print(x.real(), " ");
		simple::support::print("]\n");

		simple::support::print("ENC: [ ");
		for(auto&& x : frequencies)
			simple::support::print(x.real(), " ");
		simple::support::print("]\n");
	}
#endif
}

int main()
{

	DCT_1(5, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x*2); });
	DCT_1(5, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x); });
	DCT_1(5, [&](auto x) { return std::sin(tau*x); });
	DCT_1(5, [&](auto x) { return std::cos(tau*x); });
	DCT_1(9, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x*2); });
	DCT_1(9, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x); });
	DCT_1(9, [&](auto x) { return std::sin(tau*x); });
	DCT_1(9, [&](auto x) { return std::cos(tau*x); });
	DCT_1(17, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x*2); });
	DCT_1(17, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x); });
	DCT_1(17, [&](auto x) { return std::sin(tau*x); });
	DCT_1(17, [&](auto x) { return std::cos(tau*x); });
	DCT_1(33, [&](auto x) { return std::sin(tau*x); });
	DCT_1(33, [&](auto x) { return std::cos(tau*x); });
	DCT_1(33, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x*2); });
	DCT_1(1025, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x*2); });
	DCT_1(1025, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x); });
	DCT_1(1025, [&](auto x) { return std::sin(tau*x); });
	DCT_1(1025, [&](auto x) { return std::cos(tau*x); });
	DCT_1(2049, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x*2); });
	DCT_1(2049, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x); });
	DCT_1(2049, [&](auto x) { return std::sin(tau*x)/2; });
	DCT_1(2049, [&](auto x) { return std::sin(tau*x); });
	DCT_1(2049, [&](auto x) { return std::cos(tau*x); });

	DCT_1(2049, [&](auto x) { return std::sin(tau*13*x) + std::cos(tau*345*x); });
	DCT_1(2049, [&](auto x) { return std::sin(tau*134*x) + std::sin(tau*345*x) + std::cos(tau*789*x); });

	auto seed = std::random_device{}();
	std::cout << "DCT_1 random test seed: " << std::hex << std::showbase << seed << std::dec << std::endl;
	simple::support::random::engine::tiny<unsigned long long> random{seed};
	std::uniform_real_distribution<double> number{};

	DCT_1(2049, [&](auto) { return number(random); });

	return 0;
}