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