API Documentation
Loading...
Searching...
No Matches
FFT.h
Go to the documentation of this file.
1#pragma once
2#include <NDEVR/Buffer.h>
3#include <NDEVR/Angle.h>
4#include <iostream>
5#include <iomanip>
6#include <complex>
7
8namespace NDEVR
9{
10
11 template<class t_type>
12 class FFTLogic
13 {
14 protected:
15 // Fast Fourier Transform
16 static void fft(Buffer<std::complex<t_type>>& in_out, Buffer<std::complex<t_type>>& temp, uint04 offset, uint04 size)
17 {
18 // base case
19 if (size == 1)
20 return;
21 if (size % 2 != 0)// Fallback to Basic Discrete Fourier Transform
22 {
23 for (uint04 k = 0; k < size; k++)
24 temp[k + offset] = in_out[k + offset];
25 dft(temp.begin(offset), in_out.begin(offset), size, 1);
26 return;
27 }
28 lib_assert(size % 2 == 0, "Bad FFT");
29
30 uint04 even_offset = offset;
31 uint04 odd_offset = offset + size / 2;
32
33 // compute FFT of even terms
34 for (uint04 k = 0; k < size / 2; k++)
35 temp[k + even_offset] = in_out[2 * k + offset];
36
37 // compute FFT of odd terms
38 for (uint04 k = 0; k < size / 2; k++)
39 temp[k + odd_offset] = in_out[2 * k + 1 + offset];
40
41 fft(temp, in_out, even_offset, size / 2);
42 fft(temp, in_out, odd_offset, size / 2);
43
44 // combine
45 for (uint04 k = 0; k < size / 2; k++)
46 {
47 fltp08 kth = -2.0 * cast<fltp08>(k) * PI<fltp08>() / size;
48 std::complex<t_type> wk(std::cos(kth), std::sin(kth));
49 in_out[k + even_offset] = temp[k + even_offset] + (wk * temp[k + odd_offset]);
50 in_out[k + odd_offset] = temp[k + even_offset] - (wk * temp[k + odd_offset]);
51 }
52 }
53 // Basic Discrete Fourier Transform
54 static void dft(const std::complex<t_type> f[], std::complex<t_type> ftilde[], uint04 N, uint04 step)
55 {
56 std::complex<t_type> w = std::polar(1.0, -2.0 * PI<t_type>() / N);
57 std::complex<t_type> wn = 1.0;
58 for (uint04 n = 0; n < N; n++)
59 {
60 if (n > 0) wn *= w;
61 std::complex<t_type> wm = 1.0;
62 ftilde[n] = 0.0;
63 for (uint04 m = 0, mpos = 0; m < N; m++, mpos += step)
64 {
65 if (m > 0) wm *= wn;
66 ftilde[n] += f[mpos] * wm;
67 }
68 }
69 }
70
71 // compute the inverse FFT of x[], assuming its length n is a power of 2
72 static void ifft(Buffer<std::complex<t_type>>& in_out, Buffer<std::complex<t_type>>& temp, uint04 offset, uint04 size)
73 {
74 // take conjugate
75 for (uint04 i = 0; i < size; i++)
76 in_out[i + offset] = std::conj(in_out[i + offset]);
77
78 // compute forward FFT
79 fft(in_out, temp, offset, size);
80
81 // take conjugate again
82 t_type factor = 1.0 / size;
83 for (uint04 i = 0; i < size; i++)
84 in_out[i + offset] = std::conj(in_out[i + offset]) * factor;
85
86 }
87
88
89
90 //======================================================================
91
92 static void iFFT(const std::complex<t_type> ftilde[], std::complex<t_type> f[], uint04 N) // Inverse Fast Fourier Transform
93 {
94 Buffer<std::complex<t_type>> ftildeConjugate;
95 ftildeConjugate.setSize(N);
96
97 for (uint04 m = 0; m < N; m++)
98 ftildeConjugate[m] = std::conj(ftilde[m]);
99
100 fft(ftildeConjugate, f, N, 1);
101
102 t_type factor = 1.0 / N;
103 for (uint04 m = 0; m < N; m++)
104 f[m] = std::conj(f[m]) * factor;
105
106 }
107 public:
108 static void FT(Buffer<std::complex<t_type>>& in_out)
109 {
110 Buffer<std::complex<t_type>> temp;
111 temp.setSize(in_out.size());
112 fft(in_out, temp, 0, temp.size());
113 }
114 static void IFT(Buffer<std::complex<t_type>>& in_out)
115 {
116 Buffer<std::complex<t_type>> temp;
117 temp.setSize(in_out.size());
118 ifft(in_out, temp, 0, temp.size());
119 }
120 };
121
122 /*static Buffer<std::complex<t_type>> fft(const Buffer<std::complex<t_type>>& f) // Fast Fourier Transform
123 {
124 // base case
125 if (f.size() == 1)
126 {
127 return { f[0] };
128 }
129 lib_assert(f.size() % 2 == 0, "Bad FFT");
130
131 Buffer<std::complex<t_type>> buffer;
132 buffer.setSize(f.size() / 2);
133 // compute FFT of even terms
134 for (uint04 k = 0; k < f.size() / 2; k++)
135 buffer[k] = f[2 * k];
136 Buffer<std::complex<t_type>> even_fft = fft(buffer);
137
138 // compute FFT of odd terms
139 for (uint04 k = 0; k < f.size() / 2; k++)
140 buffer[k] = f[2 * k + 1];
141 Buffer<std::complex<t_type>> odd_fft = fft(buffer);
142
143 // combine
144 Buffer<std::complex<t_type>> y;
145 y.setSize(f.size());
146 for (uint04 k = 0; k < f.size() / 2; k++)
147 {
148 fltp08 kth = -2.0 * cast<fltp08>(k) * Angle::PI<fltp08>() / f.size();
149 std::complex<t_type> wk(std::cos(kth), std::sin(kth));
150 y[k ] = even_fft[k] + (wk * odd_fft[k]);
151 y[k + f.size() / 2] = even_fft[k] - (wk * odd_fft[k]);
152 }
153 return y;
154 }
155
156 // compute the inverse FFT of x[], assuming its length n is a power of 2
157 static Buffer<std::complex<t_type>> ifft(const Buffer<std::complex<t_type>>& x)
158 {
159 Buffer<std::complex<t_type>> y;
160 y.setSize(x.size());
161
162 // take conjugate
163 for (uint04 i = 0; i < x.size(); i++)
164 y[i] = std::conj(x[i]);
165
166 // compute forward FFT
167 y = fft(y);
168
169 // take conjugate again
170 t_type factor = 1.0 / x.size();
171 for (uint04 i = 0; i < x.size(); i++)
172 y[i] = std::conj(y[i]) * factor;
173 return y;
174
175 }
176
177 //======================================================================
178
179 static void DFT(const std::complex<t_type> f[], std::complex<t_type> ftilde[], uint04 N, uint04 step) // Basic Discrete Fourier Transform
180 {
181 std::complex<t_type> w = std::polar(1.0, -2.0 * Angle::PI<t_type>() / N);
182 std::complex<t_type> wn = 1.0;
183 for (uint04 n = 0; n < N; n++)
184 {
185 if (n > 0) wn *= w;
186 std::complex<t_type> wm = 1.0;
187 ftilde[n] = 0.0;
188 for (uint04 m = 0, mpos = 0; m < N; m++, mpos += step)
189 {
190 if (m > 0) wm *= wn;
191 ftilde[n] += f[mpos] * wm;
192 }
193 }
194 }
195
196 //======================================================================
197
198 static void iFFT(const std::complex<t_type> ftilde[], std::complex<t_type> f[], uint04 N) // Inverse Fast Fourier Transform
199 {
200 Buffer<std::complex<t_type>> ftildeConjugate;
201 ftildeConjugate.setSize(N);
202
203 for (uint04 m = 0; m < N; m++)
204 ftildeConjugate[m] = std::conj(ftilde[m]);
205
206 FFT(ftildeConjugate, f, N, 1);
207
208 t_type factor = 1.0 / N;
209 for (uint04 m = 0; m < N; m++)
210 f[m] = std::conj(f[m]) * factor;
211
212 }
213 public:
214 static Buffer<std::complex<t_type>> FFT(const Buffer<std::complex<t_type>>& complex, uint04 step = 1)
215 {
216 Buffer<std::complex<t_type>> output;
217 output.setSize(complex.size());
218 FFT(complex.ptr(), output.ptr(), complex.size(), step);
219 return output;
220 }
221 static Buffer<std::complex<t_type>> IFFT(const Buffer<std::complex<t_type>>& complex)
222 {
223 return ifft(complex);
224 }
225 };*/
226}
lib_assert
#define lib_assert(expression, message)
Asserts some logic in the code. Disabled in non debug mode by default. Can be re-enabled in release u...
Definition LibAssert.h:70
NDEVR::Buffer
The equivelent of std::vector but with a bit more control. The basic array unit of the library.
Definition Buffer.hpp:64
NDEVR::Buffer::size
constexpr t_index_type size() const
Definition Buffer.hpp:1461
NDEVR::Buffer::setSize
void setSize(t_index_type new_size)
Definition Buffer.hpp:1413
NDEVR::FFTLogic
Definition FFT.h:13
NDEVR::FFTLogic::iFFT
static void iFFT(const std::complex< t_type > ftilde[], std::complex< t_type > f[], uint04 N)
Definition FFT.h:92
NDEVR::FFTLogic::fft
static void fft(Buffer< std::complex< t_type > > &in_out, Buffer< std::complex< t_type > > &temp, uint04 offset, uint04 size)
Definition FFT.h:16
NDEVR::FFTLogic::IFT
static void IFT(Buffer< std::complex< t_type > > &in_out)
Definition FFT.h:114
NDEVR::FFTLogic::dft
static void dft(const std::complex< t_type > f[], std::complex< t_type > ftilde[], uint04 N, uint04 step)
Definition FFT.h:54
NDEVR::FFTLogic::ifft
static void ifft(Buffer< std::complex< t_type > > &in_out, Buffer< std::complex< t_type > > &temp, uint04 offset, uint04 size)
Definition FFT.h:72
NDEVR::FFTLogic::FT
static void FT(Buffer< std::complex< t_type > > &in_out)
Definition FFT.h:108
NDEVR
Definition ACIColor.h:37
NDEVR::uint04
uint32_t uint04
-Defines an alias representing a 4 byte, unsigned integer -Can represent exact integer values 0 throu...
Definition BaseValues.hpp:120
NDEVR::cast
constexpr t_to cast(const Angle< t_from > &value)
Definition Angle.h:514
NDEVR::fltp08
double fltp08
Defines an alias representing an 8 byte floating-point number.
Definition BaseValues.hpp:181