Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F87311669
fft.hh
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Oct 11, 22:50
Size
4 KB
Mime Type
text/x-c
Expires
Sun, Oct 13, 22:50 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
21574664
Attached To
R7581 SP4E Homework
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
];}
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
Log In to Comment