#include #include #include #include #include #include #include #include #include #include "wav.hpp" namespace fftw { using complex = std::complex; static_assert(sizeof(fftw::complex) == sizeof(fftw_complex), "Types should match"); struct plan { // fftw_plan fftw_plan_dft_1d(int n, fftw_complex *in, fftw_complex *out, int sign, unsigned flags); // These pointers can be equal, indicating an in-place transform. plan(int N, fftw_complex* in, fftw_complex* out){ p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE); } plan(plan &&) = delete; plan(plan const &) = delete; void operator()(){ fftw_execute(p); } ~plan(){ fftw_destroy_plan(p); } private: // is actually a fftw_plan_s* fftw_plan p; }; // the k-th output corresponds to the frequency k/n // The frequency -k/n is the same as the frequency (n-k)/n. // For real-valued data the second half is redundant: out[i] is the conjugate of out[n-i] } std::vector generate_window(unsigned int N){ std::vector window(N); auto a0 = 0.3635819; auto a1 = 0.4891775; auto a2 = 0.1365995; auto a3 = 0.0106411; for(int n = 0; n < N; ++n){ window[n] = a0 - a1*std::cos(2*M_PI*n / double(N)) + a2*std::cos(4*M_PI*n / double(N)) - a3*std::cos(6*M_PI*n / double(N)); } return window; } int main(){ const int N = 1 << 13; const int H = 1 << 13; const int W = 1 << 9; auto sqrtn = 1.0 / std::sqrt(N); std::vector input(N); auto window = generate_window(N); // in place fftw::plan p(N, reinterpret_cast(input.data()), reinterpret_cast(input.data())); png::gray_ostream image(W, H, "dmt_song.png"); wav::handle sound("samples/dmt_song.wav"); timer t("generating image"); for(int y = 0; y < H; ++y){ sound.rewind(y * 512); auto i = 0; auto end = input.size(); // copy a part to the array for(auto&& x : sound){ input[i++] = x * window[i]; if(i == end) break; } // pad with zero while(i != end){ input[i++] = 0.0; } // fft tah shit p(); // write to image auto pos = W; for(auto&& x : input){ if(!pos--) break; auto v = sqrtn * std::abs(x); image << png::gray_ostream::pixel(v); } } }