You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
101 lines
2.3 KiB
101 lines
2.3 KiB
#include <iostream>
|
|
#include <vector>
|
|
#include <complex>
|
|
#include <algorithm>
|
|
#include <iterator>
|
|
#include <cmath>
|
|
|
|
#include <fftw3.h>
|
|
#include <png.hpp>
|
|
#include <utilities.hpp>
|
|
|
|
#include "wav.hpp"
|
|
|
|
namespace fftw {
|
|
using complex = std::complex<double>;
|
|
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<double> generate_window(unsigned int N){
|
|
std::vector<double> 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<fftw::complex> input(N);
|
|
auto window = generate_window(N);
|
|
|
|
// in place
|
|
fftw::plan p(N, reinterpret_cast<fftw_complex*>(input.data()), reinterpret_cast<fftw_complex*>(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);
|
|
}
|
|
}
|
|
}
|
|
|