1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
|
#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <string.h>
#include "HalideBuffer.h"
#include "fft_forward_c2c.h"
#include "fft_forward_r2c.h"
#include "fft_inverse_c2c.h"
#include "fft_inverse_c2r.h"
namespace {
const float kPi = 3.14159265358979310000f;
const int32_t kSize = 16;
} // namespace
using Halide::Runtime::Buffer;
// Note that real_buffer() is 3D (with the 3rd dimension having extent 0)
// because the fft is written generically to require 3D inputs, even when they are real.
// Hence, the resulting buffer must be accessed with buf(i, j, 0).
Buffer<float, 3> real_buffer(int32_t y_size = kSize) {
return Buffer<float, 3>::make_interleaved(kSize, y_size, 1);
}
Buffer<float, 3> complex_buffer(int32_t y_size = kSize) {
return Buffer<float, 3>::make_interleaved(kSize, y_size, 2);
}
float &re(Buffer<float, 3> &b, int x, int y) {
return b(x, y, 0);
}
float &im(Buffer<float, 3> &b, int x, int y) {
return b(x, y, 1);
}
float re(const Buffer<float, 3> &b, int x, int y) {
return b(x, y, 0);
}
float im(const Buffer<float, 3> &b, int x, int y) {
return b(x, y, 1);
}
int main(int argc, char **argv) {
std::cout << std::fixed << std::setprecision(2);
// Forward real to complex test.
{
std::cout << "Forward real to complex test.\n";
float signal_1d[kSize];
for (size_t i = 0; i < kSize; i++) {
signal_1d[i] = 0;
for (size_t k = 1; k < 5; k++) {
signal_1d[i] += cos(2 * kPi * (k * (i / (float)kSize) + (k / 16.0f)));
}
}
auto in = real_buffer();
for (int j = 0; j < kSize; j++) {
for (int i = 0; i < kSize; i++) {
in(i, j, 0) = signal_1d[i] + signal_1d[j];
}
}
auto out = complex_buffer(kSize / 2 + 1);
int halide_result;
halide_result = fft_forward_r2c(in, out);
if (halide_result != 0) {
std::cerr << "fft_forward_r2c failed returning " << halide_result << "\n";
exit(1);
}
for (size_t i = 1; i < 5; i++) {
// Check horizontal bins
float real = re(out, i, 0);
float imaginary = im(out, i, 0);
float magnitude = sqrt(real * real + imaginary * imaginary);
if (fabs(magnitude - .5f) > .001) {
std::cerr << "fft_forward_r2c bad magnitude for horizontal bin " << i << ":" << magnitude << "\n";
exit(1);
}
float phase_angle = atan2(imaginary, real);
if (fabs(phase_angle - (i / 16.0f) * 2 * kPi) > .001) {
std::cerr << "fft_forward_r2c bad phase angle for horizontal bin " << i << ": " << phase_angle << "\n";
exit(1);
}
// Check vertical bins
real = re(out, 0, i);
imaginary = im(out, 0, i);
magnitude = sqrt(real * real + imaginary * imaginary);
if (fabs(magnitude - .5f) > .001) {
std::cerr << "fft_forward_r2c bad magnitude for vertical bin " << i << ":" << magnitude << "\n";
exit(1);
}
phase_angle = atan2(imaginary, real);
if (fabs(phase_angle - (i / 16.0f) * 2 * kPi) > .001) {
std::cerr << "fft_forward_r2c bad phase angle for vertical bin " << i << ": " << phase_angle << "\n";
exit(1);
}
}
// Check all other components are close to zero.
for (size_t j = 0; j < kSize / 2 + 1; j++) {
for (size_t i = 0; i < kSize; i++) {
// The first four non-DC bins in x and y have non-zero
// values. The horizontal ones are mirrored into the
// negative frequency components as well.
if (!((j == 0 && ((i > 0 && i < 5) || (i > kSize - 5))) ||
(i == 0 && j > 0 && j < 5))) {
float real = re(out, i, j);
float imaginary = im(out, i, j);
if (fabs(real) > .001) {
std::cerr << "fft_forward_r2c real component at (" << i << ", " << j << ") is non-zero: " << real << "\n";
exit(1);
}
if (fabs(imaginary) > .001) {
std::cerr << "fft_forward_r2c imaginary component at (" << i << ", " << j << ") is non-zero: " << imaginary << "\n";
exit(1);
}
}
}
}
}
// Inverse complex to real test.
{
std::cout << "Inverse complex to real test.\n";
auto in = complex_buffer();
in.fill(0);
// There are four components that get summed to form the magnitude, which we want to be 1.
// The components are each of the positive and negative frequencies and each of the
// real and complex components. The +/- frequencies sum algebraically and the complex
// components contribute to the magnitude as the sides of triangle like any 2D vector.
float term_magnitude = 1.0f / (2.0f * sqrt(2.0f));
re(in, 1, 0) = term_magnitude;
im(in, 1, 0) = term_magnitude;
// Negative frequencies count backward from end, no DC term
re(in, kSize - 1, 0) = term_magnitude;
im(in, kSize - 1, 0) = -term_magnitude; // complex conjugate
auto out = real_buffer();
int halide_result;
halide_result = fft_inverse_c2r(in, out);
if (halide_result != 0) {
std::cerr << "fft_inverse_c2r failed returning " << halide_result << "\n";
exit(1);
}
for (size_t j = 0; j < kSize; j++) {
for (size_t i = 0; i < kSize; i++) {
float sample = out(i, j, 0);
float expected = cos(2 * kPi * (i / 16.0f + .125f));
if (fabs(sample - expected) > .001) {
std::cerr << "fft_inverse_c2r mismatch at (" << i << ", " << j << ") " << sample << " vs. " << expected << "\n";
exit(1);
}
}
}
}
// Forward complex to complex test.
{
std::cout << "Forward complex to complex test.\n";
auto in = complex_buffer();
float signal_1d_real[kSize];
float signal_1d_complex[kSize];
for (size_t i = 0; i < kSize; i++) {
signal_1d_real[i] = 0;
signal_1d_complex[i] = 0;
for (size_t k = 1; k < 5; k++) {
signal_1d_real[i] += cos(2 * kPi * (k * (i / (float)kSize) + (k / 16.0f)));
signal_1d_complex[i] += sin(2 * kPi * (k * (i / (float)kSize) + (k / 16.0f)));
}
}
for (int j = 0; j < kSize; j++) {
for (int i = 0; i < kSize; i++) {
re(in, i, j) = signal_1d_real[i] + signal_1d_real[j];
im(in, i, j) = signal_1d_complex[i] + signal_1d_complex[j];
}
}
auto out = complex_buffer();
int halide_result;
halide_result = fft_forward_c2c(in, out);
if (halide_result != 0) {
std::cerr << "fft_forward_c2c failed returning " << halide_result << "\n";
exit(1);
}
for (size_t i = 1; i < 5; i++) {
// Check horizontal bins
float real = re(out, i, 0);
float imaginary = im(out, i, 0);
float magnitude = sqrt(real * real + imaginary * imaginary);
if (fabs(magnitude - 1.0f) > .001) {
std::cerr << "fft_forward_c2c bad magnitude for horizontal bin " << i << ":" << magnitude << "\n";
exit(1);
}
float phase_angle = atan2(imaginary, real);
if (fabs(phase_angle - (i / 16.0f) * 2 * kPi) > .001) {
std::cerr << "fft_forward_c2c bad phase angle for horizontal bin " << i << ": " << phase_angle << "\n";
exit(1);
}
// Check vertical bins
real = re(out, 0, i);
imaginary = im(out, 0, i);
magnitude = sqrt(real * real + imaginary * imaginary);
if (fabs(magnitude - 1.0f) > .001) {
std::cerr << "fft_forward_c2c bad magnitude for vertical bin " << i << ":" << magnitude << "\n";
exit(1);
}
phase_angle = atan2(imaginary, real);
if (fabs(phase_angle - (i / 16.0f) * 2 * kPi) > .001) {
std::cerr << "fft_forward_c2c bad phase angle for vertical bin " << i << ": " << phase_angle << "\n";
exit(1);
}
}
// Check all other components are close to zero.
for (size_t j = 0; j < kSize; j++) {
for (size_t i = 0; i < kSize; i++) {
// The first four non-DC bins in x and y have non-zero
// values. The input is chose so the mirrored negative
// frequency components are all zero due to
// interference of the real and complex parts.
if (!((j == 0 && (i > 0 && i < 5)) ||
(i == 0 && j > 0 && j < 5))) {
float real = re(out, i, j);
float imaginary = im(out, i, j);
if (fabs(real) > .001) {
std::cerr << "fft_forward_c2c real component at (" << i << ", " << j << ") is non-zero: " << real << "\n";
exit(1);
}
if (fabs(imaginary) > .001) {
std::cerr << "fft_forward_c2c imaginary component at (" << i << ", " << j << ") is non-zero: " << imaginary << "\n";
exit(1);
}
}
}
}
}
// Inverse complex to complex test.
{
std::cout << "Inverse complex to complex test.\n";
auto in = complex_buffer();
in.fill(0);
re(in, 1, 0) = .5f;
im(in, 1, 0) = .5f;
re(in, kSize - 1, 0) = .5f;
im(in, kSize - 1, 0) = .5f; // Not conjugate. Result will not be real
auto out = complex_buffer();
int halide_result;
halide_result = fft_inverse_c2c(in, out);
if (halide_result != 0) {
std::cerr << "fft_inverse_c2c failed returning " << halide_result << "\n";
exit(1);
}
for (size_t j = 0; j < kSize; j++) {
for (size_t i = 0; i < kSize; i++) {
float real_sample = re(out, i, j);
float imaginary_sample = im(out, i, j);
float real_expected = 1 / sqrt(2) * (cos(2 * kPi * (i / 16.0f + .125)) + cos(2 * kPi * (i * (kSize - 1) / 16.0f + .125)));
float imaginary_expected = 1 / sqrt(2) * (sin(2 * kPi * (i / 16.0f + .125)) + sin(2 * kPi * (i * (kSize - 1) / 16.0f + .125)));
if (fabs(real_sample - real_expected) > .001) {
std::cerr << "fft_inverse_c2c real mismatch at (" << i << ", " << j << ") " << real_sample << " vs. " << real_expected << "\n";
exit(1);
}
if (fabs(imaginary_sample - imaginary_expected) > .001) {
std::cerr << "fft_inverse_c2c imaginary mismatch at (" << i << ", " << j << ") " << imaginary_sample << " vs. " << imaginary_expected << "\n";
exit(1);
}
}
}
}
std::cout << "Success!\n";
exit(0);
}
|
