Page MenuHomec4science
Files F66365926

fft.hh
No OneTemporary
Actions

File Metadata

Created
Mon, Jun 10, 02:24

fft.hh
View Options

#ifndef FFT_HH
#define FFT_HH
/* ------------------------------------------------------ */
#include "matrix.hh"
#include "my_types.hh"
#include <fftw3.h>
/* ------------------------------------------------------ */
struct FFT {
static Matrix<complex> transform(Matrix<complex> &m);
static Matrix<complex> itransform(Matrix<complex> &m);
static Matrix<std::complex<int>> computeFrequencies(int size);
};
/* ------------------------------------------------------ */
inline Matrix<complex> FFT::transform(Matrix<complex> &m_in) {
fftw_plan p;
int N = m_in.rows(); // same as int n1 = m_in.cols();
Matrix<complex> m_out(N);
// reinterpret cast
complex *c = m_in.data();
auto *c_in = reinterpret_cast<fftw_complex *>(c);
c = m_out.data();
auto *c_out = reinterpret_cast<fftw_complex *>(c);
p = fftw_plan_dft_2d(N, N, c_in, c_out, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p);
fftw_destroy_plan(p);
return m_out;
}
/* ------------------------------------------------------ */
inline Matrix<complex> FFT::itransform(Matrix<complex> &m_in) {
fftw_plan p;
int N = m_in.rows(); // same as int n1 = m_in.cols();
Matrix<complex> m_out(N);
// reinterpret cast
complex *c = m_in.data();
auto *c_in = reinterpret_cast<fftw_complex *>(c);
c = m_out.data();
auto *c_out = reinterpret_cast<fftw_complex *>(c);
p = fftw_plan_dft_2d(N, N, c_in, c_out, FFTW_BACKWARD, FFTW_ESTIMATE);
fftw_execute(p);
fftw_destroy_plan(p);
double rescale = 1. / (N * N);
for (auto &&entry : index(m_out)) {
auto &val = std::get<2>(entry);
val *= rescale;
}
return m_out;
}
/* ------------------------------------------------------ */
/* ------------------------------------------------------ */
inline Matrix<std::complex<int>> FFT::computeFrequencies(int size) {}
#endif // FFT_HH

Event Timeline

c4science · Help