Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F86076016
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 4, 02:57
Size
4 KB
Mime Type
text/x-c
Expires
Sun, Oct 6, 02:57 (2 d)
Engine
blob
Format
Raw Data
Handle
21341608
Attached To
R7571 SP4E-TB-TL-FR
fft.hh
View Options
#ifndef FFT_HH
#define FFT_HH
/* ------------------------------------------------------ */
#include "matrix.hh"
#include "my_types.hh"
#include <fftw3.h>
/* ------------------------------------------------------ */
/**
This file defines an FFT structure that is a wrapper around the FFTW library.
In particular, it defines 3 functions to work with signals of complex numbers:
- the transform() function to compute the forward FFT of a complex signal
- the itransform() function to compute the inverse FFT of a complex signal
- the computeFrequencies() function to compute the sample frequencies of the
forward FFT of a complex signal (i.e. the coordinates of the signal in the Fourier space).
The Laplacian of these frequencies is then used to sovle the heat equation in the Fourier space.
*/
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
)
{
// Get matrix size
int
N
=
m_in
.
size
();
// Declare input and output 2D signals
fftw_complex
*
in
,
*
out
;
in
=
(
fftw_complex
*
)
fftw_malloc
(
sizeof
(
fftw_complex
)
*
N
*
N
);
out
=
(
fftw_complex
*
)
fftw_malloc
(
sizeof
(
fftw_complex
)
*
N
*
N
);
// Populate input 2D signal
for
(
int
i
=
0
;
i
<
N
;
++
i
)
{
for
(
int
j
=
0
;
j
<
N
;
++
j
)
{
in
[
i
*
N
+
j
][
0
]
=
m_in
(
i
,
j
).
real
();
in
[
i
*
N
+
j
][
1
]
=
m_in
(
i
,
j
).
imag
();
}
}
// Create, execute and destroy FFT plan
fftw_plan
p
=
fftw_plan_dft_2d
(
N
,
N
,
in
,
out
,
FFTW_FORWARD
,
FFTW_ESTIMATE
);
fftw_execute
(
p
);
// Initialize and populate transformed matrix from output 2D signal
Matrix
<
complex
>
m_out
(
N
);
for
(
int
i
=
0
;
i
<
N
;
++
i
)
{
for
(
int
j
=
0
;
j
<
N
;
++
j
)
{
m_out
(
i
,
j
)
=
complex
(
out
[
i
*
N
+
j
][
0
],
out
[
i
*
N
+
j
][
1
]);
}
}
// Destroy plan and fftw signals
fftw_destroy_plan
(
p
);
fftw_free
(
in
);
fftw_free
(
out
);
// Return FFT matrix
return
m_out
;
}
/* ------------------------------------------------------ */
inline
Matrix
<
complex
>
FFT
::
itransform
(
Matrix
<
complex
>&
m_in
)
{
// Get matrix size
int
N
=
m_in
.
size
();
// Declare input and output 2D signals
fftw_complex
*
in
,
*
out
;
in
=
(
fftw_complex
*
)
fftw_malloc
(
sizeof
(
fftw_complex
)
*
N
*
N
);
out
=
(
fftw_complex
*
)
fftw_malloc
(
sizeof
(
fftw_complex
)
*
N
*
N
);
// Populate input 2D signal
for
(
int
i
=
0
;
i
<
N
;
++
i
)
{
for
(
int
j
=
0
;
j
<
N
;
++
j
)
{
in
[
i
*
N
+
j
][
0
]
=
m_in
(
i
,
j
).
real
();
in
[
i
*
N
+
j
][
1
]
=
m_in
(
i
,
j
).
imag
();
}
}
// Create, execute and destroy inverse FFT plan
fftw_plan
p
=
fftw_plan_dft_2d
(
N
,
N
,
in
,
out
,
FFTW_BACKWARD
,
FFTW_ESTIMATE
);
fftw_execute
(
p
);
// Initialize and populate transformed matrix from output 2D signal
Matrix
<
complex
>
m_out
(
N
);
for
(
int
i
=
0
;
i
<
N
;
++
i
)
{
for
(
int
j
=
0
;
j
<
N
;
++
j
)
{
// std::cout << "Imag: " << out[i*N + j][1]/(N * N) << std::endl;
m_out
(
i
,
j
)
=
complex
(
out
[
i
*
N
+
j
][
0
],
out
[
i
*
N
+
j
][
1
]);
m_out
(
i
,
j
)
/=
(
N
*
N
);
// dividing by factor N2 needed for 2D arrays
}
}
// Destroy plan and fftw signals
fftw_destroy_plan
(
p
);
fftw_free
(
in
);
fftw_free
(
out
);
// Return FFT matrix
return
m_out
;
}
/* ------------------------------------------------------ */
/* ------------------------------------------------------ */
inline
Matrix
<
std
::
complex
<
int
>>
FFT
::
computeFrequencies
(
int
size
)
{
// Generate 1D frequency vector
std
::
vector
<
float
>
vec
(
size
);
if
(
size
%
2
==
0
)
{
// if size is even
for
(
int
i
=
0
;
i
<
size
/
2
;
++
i
)
{
vec
[
i
]
=
2
*
M_PI
*
(
float
)
i
;
}
for
(
int
i
=
size
/
2
;
i
<
size
;
++
i
)
{
vec
[
i
]
=
2
*
M_PI
*
(
-
size
/
2
+
((
float
)
i
-
size
/
2
));
}
}
else
{
// if size is odd
for
(
int
i
=
0
;
i
<
(
size
-
1
)
/
2
;
++
i
)
{
vec
[
i
]
=
2
*
M_PI
*
(
float
)
i
;
}
for
(
int
i
=
(
size
-
1
)
/
2
;
i
<
size
;
++
i
)
{
vec
[
i
]
=
2
*
M_PI
*
(
-
(
size
-
1
)
/
2
+
((
float
)
i
-
(
size
-
1
)
/
2
));
}
}
// Check the 1D coordinate vecor on Fourier space
// std::cout << "computeFrequencies, vect: " << std::endl;
// for (int i = 0; i<vec.size(); ++i){
// std::cout << "vec(" << i << "): " << vec.at(i) << std::endl;
// }
// Populate 2D matrix with frequencies vector
Matrix
<
std
::
complex
<
int
>>
m_out
(
size
);
for
(
int
i
=
0
;
i
<
size
;
++
i
)
{
for
(
int
j
=
0
;
j
<
size
;
++
j
)
{
m_out
(
i
,
j
)
=
std
::
complex
<
int
>
(
vec
[
j
],
0
);
// std::cout << "(" << i << "," << j << "): " << vec.at(j) << " ";
}
// std::cout << std::endl;
}
return
m_out
;
}
#endif
// FFT_HH
Event Timeline
Log In to Comment