libsimple_compress/unit_tests/fft.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);

template <typename Waveform>
void DFT(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);

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

	for(std::size_t i = 0; i != wave.size(); ++i)
	{
		for(std::size_t j = 0; j != wave.size(); ++j)
		{
			control[i].real( control[i].real() + wave[j] * std::cos(tau*j*i/wave.size()) / control.size());
			control[i].imag( control[i].imag() + wave[j] * std::sin(tau*j*i/wave.size()) / control.size());
		}
	}

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

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

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

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

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


#if defined SIMPLE_SUPPORT_DEBUG_HPP
	if(size <= 32)
	{
		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(std::abs(x), " ");
		simple::support::print("]\n");

		simple::support::print("CNR: [ ");
		for(auto&& x : control)
			simple::support::print(std::abs(x), " ");
		simple::support::print("]\n");
	}
#endif
}

int main()
{

	DFT(2, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x*2); });
	DFT(2, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x); });
	DFT(2, [&](auto x) { return std::sin(tau*x); });
	DFT(2, [&](auto x) { return std::cos(tau*x); });
	DFT(4, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x*2); });
	DFT(4, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x); });
	DFT(4, [&](auto x) { return std::sin(tau*x); });
	DFT(4, [&](auto x) { return std::cos(tau*x); });
	DFT(8, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x*2); });
	DFT(8, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x); });
	DFT(8, [&](auto x) { return std::sin(tau*x); });
	DFT(8, [&](auto x) { return std::cos(tau*x); });
	DFT(16, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x*2); });
	DFT(16, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x); });
	DFT(16, [&](auto x) { return std::sin(tau*x); });
	DFT(16, [&](auto x) { return std::cos(tau*x); });
	DFT(1024, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x*2); });
	DFT(1024, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x); });
	DFT(1024, [&](auto x) { return std::sin(tau*x); });
	DFT(1024, [&](auto x) { return std::cos(tau*x); });
	DFT(2048, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x*2); });
	DFT(2048, [&](auto x) { return std::sin(tau*x) + std::cos(tau*x); });
	DFT(2048, [&](auto x) { return std::sin(tau*x); });
	DFT(2048, [&](auto x) { return std::cos(tau*x); });

	DFT(2048, [&](auto x) { return std::sin(tau*13*x) + std::cos(tau*345*x); });
	DFT(2048, [&](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 << "DFT 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{};

	DFT(2048, [&](auto) { return number(random); });

	return 0;
}