Page MenuHomec4science
Files F60145253

fft.hh
No OneTemporary
Actions

File Metadata

Created
Sat, Apr 27, 20:26

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) {
// Implementation for the Forward FFT
UInt N = m_in.cols(); // The sampled size
UInt Sz = N*N; // Size of the Matrix NxN
fftw_plan P; // Initialize a plan
// 2D Complex Matricies alocation
complex *m_out_f = new complex[Sz], *m_in_f= new complex[Sz];
// Casting Values
for (auto&& entry : index(m_in)){
int i = std::get<0>(entry);
int j = std::get<1>(entry);
auto& val = std::get<2>(entry);
m_in_f[i*N+j] = val;}
/// I've obtained this (following lines) solution from this link (Not sure but should work):
/// http://www.fftw.org/fftw3_doc/Complex-numbers.html
/// https://www.reddit.com/r/cpp/comments/jfk9q/is_there_a_quick_way_to_convert_fftw_complex_to_a/
P = fftw_plan_dft_2d(
N,
N,
reinterpret_cast<fftw_complex*>(m_in_f),
reinterpret_cast<fftw_complex*>(m_out_f),
FFTW_FORWARD,
FFTW_ESTIMATE // Your flags here.
);
fftw_execute(P);
// Casting the final output matrix -> m_out_f_f
Matrix<complex> m_out_f_f(N);
for (auto&& entry : index(m_out_f_f)) {
int i = std::get<0>(entry);
int j = std::get<1>(entry);
auto& val = std::get<2>(entry);
val = m_out_f[i*N+j];}
// Dealocation (FFTW Destructor) -> Not sure about this!
fftw_destroy_plan(P);
fftw_free(m_in_f);
fftw_free(m_out_f);
// Alternative method that should de-allocate the vectors
//delete [] m_in_f;
//delete [] m_out_f;
return m_out_f_f;
}
/* ------------------------------------------------------ */
inline Matrix<complex> FFT::itransform(Matrix<complex>& m_in) {
// Goes the same as the previouse one
// The results are not normalized in this function (Divide by N^2 the output)
UInt N = m_in.cols();
UInt Sz = N*N;
fftw_plan P;
complex *m_out_f = new complex[Sz], *m_in_f= new complex[Sz];
for (auto&& entry : index(m_in)){
int i = std::get<0>(entry);
int j = std::get<1>(entry);
auto& val = std::get<2>(entry);
m_in_f[i*N+j] = val;}
P = fftw_plan_dft_2d(
N,
N,
reinterpret_cast<fftw_complex*>(m_in_f),
reinterpret_cast<fftw_complex*>(m_out_f),
FFTW_BACKWARD,
FFTW_ESTIMATE
);
fftw_execute(P);
Matrix<complex> m_out_f_f(N);
for (auto&& entry : index(m_out_f_f)) {
int i = std::get<0>(entry);
int j = std::get<1>(entry);
auto& val = std::get<2>(entry);
val = m_out_f[i*N+j];} // Normalization needed
//Normalizing the results
m_out_f_f/=Sz;
fftw_destroy_plan(P);
fftw_free(m_in_f);
fftw_free(m_out_f);
// Alternative method that should de-allocate the vectors
//delete [] m_in_f;
//delete [] m_out_f;
return m_out_f_f;
}
/* ------------------------------------------------------ */
/* ------------------------------------------------------ */
inline Matrix<std::complex<int>> FFT::computeFrequencies(int size) {
/// The i stands for the real part, whereas the j stands for the imaginary part
/// Explanation of procedure on the link:
/// http://www.fftw.org/fftw3.pdf?fbclid=IwAR3bgzRkoSBACIJnMNfk79RX29rsVkIhWgGYg_iJcPvTM_8wLLN045ntywk
// specially the figure 2.1 -- The last member is treated spatially!
// x2 to run over the whole size
// to account for the mirror effect of DFT symmetric matrix! (To be verified)
int indicies;
if (size%2 == 0)
indicies = size/2;
else
indicies = size/2+1;
// New output matrix
Matrix<std::complex<int>> final_frequencies(size);
for (auto&& entry : index(final_frequencies)) {
int i = std::get<0>(entry);
int j = std::get<1>(entry);
auto& val = std::get<2>(entry);
if(i < indicies)
val.real(i);
else
val.real(i - size);
if(j < indicies)
val.imag(j);
else
val.imag(j - size);
}
return final_frequencies;
}
#endif // FFT_HH

Event Timeline

c4science · Help