Homec4science
Diffusion Homework_sp4e_Peruzzo_SáezUribe (master) Browse

R9490/homework3/src/fft.hh
6a14377fc7d6master

Last Change
Raw File

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) {
// changing the type from Matrix to fftw_complex but handling the memory with pointers
auto in = (fftw_complex *)m_in.data();
// create a empty output object of type matrix
Matrix<complex> m_out(m_in.size());
// changing the type from Matrix to fftw_complex
auto out = (fftw_complex *)m_out.data();
// creating the plan
fftw_plan p;
// apply the plan
p = fftw_plan_dft_2d(m_in.rows(),m_in.cols(), in, out, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p); /* repeat as needed */
fftw_destroy_plan(p);
return m_out;}
/* ------------------------------------------------------ */
inline Matrix<complex> FFT::itransform(Matrix<complex>& m_in) {
// changing the type from Matrix to fftw_complex but handling the memory with pointers
auto in = (fftw_complex *)m_in.data();
// create a empty output object of type matrix
Matrix<complex> m_out(m_in.size());
// changing the type from Matrix to fftw_complex
auto out = (fftw_complex *)m_out.data();
// creating the plan
fftw_plan p;
// apply the plan
p = fftw_plan_dft_2d(m_in.rows(),m_in.cols(), in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
fftw_execute(p); /* repeat as needed */
fftw_destroy_plan(p);
return m_out;
}
/* ------------------------------------------------------ */
/* ------------------------------------------------------ */
inline Matrix<std::complex<int>> FFT::computeFrequencies(int size) {
Matrix<std::complex<int>> freq(size);
for (auto&& entry : index(freq)) {
int i = std::get<0>(entry);
int j = std::get<1>(entry);
auto& val = std::get<2>(entry);
if (i <= size/2 - 1 && j <= size/2 - 1)
val = {i,j};
else if (i > size/2 - 1 && j <= size/2 - 1)
val = {i-size,j};
else if (i <= size/2 - 1 && j > size/2 - 1)
val = {i,j-size};
else
val = {i-size,j-size};
}
return freq;
}
#endif // FFT_HH
c4science · Help