|
@@ -0,0 +1,128 @@
|
|
|
+#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); });
|
|
|
+ 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;
|
|
|
+}
|