]> git.sesse.net Git - movit/blob - fft_pass_effect_test.cpp
Support other output formats than GL_RGBA16F; was easier than originally feared.
[movit] / fft_pass_effect_test.cpp
1 // Unit tests for FFTPassEffect.
2
3 #include <math.h>
4 #include <stdlib.h>
5 #include <string.h>
6 #include <epoxy/gl.h>
7 #include <gtest/gtest.h>
8
9 #include "effect_chain.h"
10 #include "fft_pass_effect.h"
11 #include "image_format.h"
12 #include "multiply_effect.h"
13 #include "test_util.h"
14
15 namespace movit {
16
17 namespace {
18
19 // Generate a random number uniformly distributed between [-1.0, 1.0].
20 float uniform_random()
21 {
22         return 2.0 * ((float)rand() / RAND_MAX - 0.5);
23 }
24
25 void setup_fft(EffectChain *chain, int fft_size, bool inverse,
26                bool add_normalizer = false,
27                FFTPassEffect::Direction direction = FFTPassEffect::HORIZONTAL)
28 {
29         assert((fft_size & (fft_size - 1)) == 0);  // Must be power of two.
30         for (int i = 1, subsize = 2; subsize <= fft_size; ++i, subsize *= 2) {
31                 Effect *fft_effect = chain->add_effect(new FFTPassEffect());
32                 bool ok = fft_effect->set_int("fft_size", fft_size);
33                 ok |= fft_effect->set_int("pass_number", i);
34                 ok |= fft_effect->set_int("inverse", inverse);
35                 ok |= fft_effect->set_int("direction", direction);
36                 assert(ok);
37         }
38
39         if (add_normalizer) {
40                 float factor[4] = { 1.0f / fft_size, 1.0f / fft_size, 1.0f / fft_size, 1.0f / fft_size };
41                 Effect *multiply_effect = chain->add_effect(new MultiplyEffect());
42                 bool ok = multiply_effect->set_vec4("factor", factor);
43                 assert(ok);
44         }
45 }
46
47 void run_fft(const float *in, float *out, int fft_size, bool inverse,
48              bool add_normalizer = false,
49              FFTPassEffect::Direction direction = FFTPassEffect::HORIZONTAL)
50 {
51         int width, height;
52         if (direction == FFTPassEffect::HORIZONTAL) {
53                 width = fft_size;
54                 height = 1;
55         } else {
56                 width = 1;
57                 height = fft_size;
58         }
59         EffectChainTester tester(in, width, height, FORMAT_RGBA_PREMULTIPLIED_ALPHA, COLORSPACE_sRGB, GAMMA_LINEAR);
60         setup_fft(tester.get_chain(), fft_size, inverse, add_normalizer, direction);
61         tester.run(out, GL_RGBA, COLORSPACE_sRGB, GAMMA_LINEAR, OUTPUT_ALPHA_FORMAT_PREMULTIPLIED);
62 }
63
64 }  // namespace
65
66 TEST(FFTPassEffectTest, ZeroStaysZero) {
67         const int fft_size = 64;
68         float data[fft_size * 4] = { 0 };
69         float out_data[fft_size * 4];
70
71         run_fft(data, out_data, fft_size, false);
72         expect_equal(data, out_data, 4, fft_size);
73
74         run_fft(data, out_data, fft_size, true);
75         expect_equal(data, out_data, 4, fft_size);
76 }
77
78 TEST(FFTPassEffectTest, Impulse) {
79         const int fft_size = 64;
80         float data[fft_size * 4] = { 0 };
81         float expected_data[fft_size * 4], out_data[fft_size * 4];
82         data[0] = 1.0;
83         data[1] = 1.2;
84         data[2] = 1.4;
85         data[3] = 3.0;
86
87         for (int i = 0; i < fft_size; ++i) {
88                 expected_data[i * 4 + 0] = data[0];
89                 expected_data[i * 4 + 1] = data[1];
90                 expected_data[i * 4 + 2] = data[2];
91                 expected_data[i * 4 + 3] = data[3];
92         }
93
94         run_fft(data, out_data, fft_size, false);
95         expect_equal(expected_data, out_data, 4, fft_size);
96
97         run_fft(data, out_data, fft_size, true);
98         expect_equal(expected_data, out_data, 4, fft_size);
99 }
100
101 TEST(FFTPassEffectTest, SingleFrequency) {
102         const int fft_size = 16;
103         float data[fft_size * 4] = { 0 };
104         float expected_data[fft_size * 4], out_data[fft_size * 4];
105         for (int i = 0; i < fft_size; ++i) {
106                 data[i * 4 + 0] = sin(2.0 * M_PI * (4.0 * i) / fft_size);
107                 data[i * 4 + 1] = 0.0;
108                 data[i * 4 + 2] = 0.0;
109                 data[i * 4 + 3] = 0.0;
110         }
111         for (int i = 0; i < fft_size; ++i) {
112                 expected_data[i * 4 + 0] = 0.0;
113                 expected_data[i * 4 + 1] = 0.0;
114                 expected_data[i * 4 + 2] = 0.0;
115                 expected_data[i * 4 + 3] = 0.0;
116         }
117         expected_data[4 * 4 + 1] = -8.0;
118         expected_data[12 * 4 + 1] = 8.0;
119
120         run_fft(data, out_data, fft_size, false, false, FFTPassEffect::HORIZONTAL);
121         expect_equal(expected_data, out_data, 4, fft_size);
122
123         run_fft(data, out_data, fft_size, false, false, FFTPassEffect::VERTICAL);
124         expect_equal(expected_data, out_data, 4, fft_size);
125 }
126
127 TEST(FFTPassEffectTest, Repeat) {
128         srand(12345);
129         for (int fft_size = 2; fft_size <= 128; fft_size *= 2) {
130                 const int num_repeats = 31;  // Prime, to make things more challenging.
131                 float data[num_repeats * fft_size * 4];
132                 float expected_data[num_repeats * fft_size * 4], out_data[num_repeats * fft_size * 4];
133
134                 for (int i = 0; i < num_repeats * fft_size * 4; ++i) {
135                         data[i] = uniform_random();
136                 }
137
138                 for (int i = 0; i < num_repeats; ++i) {
139                         run_fft(data + i * fft_size * 4, expected_data + i * fft_size * 4, fft_size, false);
140                 }
141
142                 {
143                         // Horizontal.
144                         EffectChainTester tester(data, num_repeats * fft_size, 1, FORMAT_RGBA_PREMULTIPLIED_ALPHA, COLORSPACE_sRGB, GAMMA_LINEAR);
145                         setup_fft(tester.get_chain(), fft_size, false);
146                         tester.run(out_data, GL_RGBA, COLORSPACE_sRGB, GAMMA_LINEAR, OUTPUT_ALPHA_FORMAT_PREMULTIPLIED);
147
148                         expect_equal(expected_data, out_data, 4, num_repeats * fft_size);
149                 }
150                 {
151                         // Vertical.
152                         EffectChainTester tester(data, 1, num_repeats * fft_size, FORMAT_RGBA_PREMULTIPLIED_ALPHA, COLORSPACE_sRGB, GAMMA_LINEAR);
153                         setup_fft(tester.get_chain(), fft_size, false, false, FFTPassEffect::VERTICAL);
154                         tester.run(out_data, GL_RGBA, COLORSPACE_sRGB, GAMMA_LINEAR, OUTPUT_ALPHA_FORMAT_PREMULTIPLIED);
155
156                         expect_equal(expected_data, out_data, 4, num_repeats * fft_size);
157                 }
158         }
159 }
160
161 TEST(FFTPassEffectTest, TwoDimensional) {  // Implicitly tests vertical.
162         srand(1234);
163         const int fft_size = 16;
164         float in[fft_size * fft_size * 4], out[fft_size * fft_size * 4], expected_out[fft_size * fft_size * 4];
165         for (int y = 0; y < fft_size; ++y) {
166                 for (int x = 0; x < fft_size; ++x) {
167                         in[(y * fft_size + x) * 4 + 0] =
168                                 sin(2.0 * M_PI * (2 * x + 3 * y) / fft_size);
169                         in[(y * fft_size + x) * 4 + 1] = 0.0;
170                         in[(y * fft_size + x) * 4 + 2] = 0.0;
171                         in[(y * fft_size + x) * 4 + 3] = 0.0;
172                 }
173         }
174         memset(expected_out, 0, sizeof(expected_out));
175
176         // This result has been verified using the fft2() function in Octave,
177         // which uses FFTW.
178         expected_out[(3 * fft_size + 2) * 4 + 1] = -128.0;
179         expected_out[(13 * fft_size + 14) * 4 + 1] = 128.0;
180
181         EffectChainTester tester(in, fft_size, fft_size, FORMAT_RGBA_PREMULTIPLIED_ALPHA, COLORSPACE_sRGB, GAMMA_LINEAR);
182         setup_fft(tester.get_chain(), fft_size, false, false, FFTPassEffect::HORIZONTAL);
183         setup_fft(tester.get_chain(), fft_size, false, false, FFTPassEffect::VERTICAL);
184         tester.run(out, GL_RGBA, COLORSPACE_sRGB, GAMMA_LINEAR, OUTPUT_ALPHA_FORMAT_PREMULTIPLIED);
185
186         expect_equal(expected_out, out, 4 * fft_size, fft_size, 0.25, 0.0005);
187 }
188
189 // The classic paper for FFT correctness testing is Funda Ergün:
190 // “Testing Multivariate Linear Functions: Overcoming the Generator Bottleneck”
191 // (http://www.cs.sfu.ca/~funda/PUBLICATIONS/stoc95.ps), which proves that
192 // testing three basic properties of FFTs guarantees that the function is
193 // correct (at least under the assumption that errors are random).
194 //
195 // We don't follow the paper directly, though, for a few reasons: First,
196 // Ergün's paper really considers _self-correcting_ systems, which may
197 // be stochastically faulty, and thus uses various relatively complicated
198 // bounds and tests we don't really need. Second, the FFTs it considers
199 // are all about polynomials over finite fields, which means that results
200 // are exact and thus easy to test; we work with floats (half-floats!),
201 // and thus need some error tolerance.
202 //
203 // So instead, we follow the implementation of FFTW, which is really the
204 // gold standard when it comes to FFTs these days. They hard-code 20
205 // testing rounds as opposed to the more complicated bounds in the paper,
206 // and have a simpler version of the third test.
207 //
208 // The error bounds are set somewhat empirically, but remember that these
209 // inputs will give frequency values as large as ~16, where 0.025 is
210 // within the 9th bit (of 11 total mantissa bits in fp16).
211 const int ergun_rounds = 20;
212
213 // Test 1: Test that FFT(a + b) = FFT(a) + FFT(b).
214 TEST(FFTPassEffectTest, ErgunLinearityTest) {
215         srand(1234);
216         const int max_fft_size = 64;
217         float a[max_fft_size * 4], b[max_fft_size * 4], sum[max_fft_size * 4];
218         float a_out[max_fft_size * 4], b_out[max_fft_size * 4], sum_out[max_fft_size * 4], expected_sum_out[max_fft_size * 4];
219         for (int fft_size = 2; fft_size <= max_fft_size; fft_size *= 2) {
220                 for (int inverse = 0; inverse <= 1; ++inverse) {
221                         for (int i = 0; i < ergun_rounds; ++i) {
222                                 for (int j = 0; j < fft_size * 4; ++j) {
223                                         a[j] = uniform_random();
224                                         b[j] = uniform_random();
225                                 }
226                                 run_fft(a, a_out, fft_size, inverse);
227                                 run_fft(b, b_out, fft_size, inverse);
228
229                                 for (int j = 0; j < fft_size * 4; ++j) {
230                                         sum[j] = a[j] + b[j];
231                                         expected_sum_out[j] = a_out[j] + b_out[j];
232                                 }
233
234                                 run_fft(sum, sum_out, fft_size, inverse);
235                                 expect_equal(expected_sum_out, sum_out, 4, fft_size, 0.03, 0.0005);
236                         }
237                 }
238         }
239 }
240
241 // Test 2: Test that FFT(delta(i)) = 1  (where delta(i) = [1 0 0 0 ...]),
242 // or more specifically, test that FFT(a + delta(i)) - FFT(a) = 1.
243 TEST(FFTPassEffectTest, ErgunImpulseTransform) {
244         srand(1235);
245         const int max_fft_size = 64;
246         float a[max_fft_size * 4], b[max_fft_size * 4];
247         float a_out[max_fft_size * 4], b_out[max_fft_size * 4], sum_out[max_fft_size * 4], expected_sum_out[max_fft_size * 4];
248         for (int fft_size = 2; fft_size <= max_fft_size; fft_size *= 2) {
249                 for (int inverse = 0; inverse <= 1; ++inverse) {
250                         for (int i = 0; i < ergun_rounds; ++i) {
251                                 for (int j = 0; j < fft_size * 4; ++j) {
252                                         a[j] = uniform_random();
253
254                                         // Compute delta(j) - a.
255                                         if (j < 4) {
256                                                 b[j] = 1.0 - a[j];
257                                         } else {
258                                                 b[j] = -a[j];
259                                         }
260                                 }
261                                 run_fft(a, a_out, fft_size, inverse);
262                                 run_fft(b, b_out, fft_size, inverse);
263
264                                 for (int j = 0; j < fft_size * 4; ++j) {
265                                         sum_out[j] = a_out[j] + b_out[j];
266                                         expected_sum_out[j] = 1.0;
267                                 }
268                                 expect_equal(expected_sum_out, sum_out, 4, fft_size, 0.025, 0.0005);
269                         }
270                 }
271         }
272 }
273
274 // Test 3: Test the time-shift property of the FFT, in that a circular left-shift
275 // multiplies the result by e^(j 2pi k/N) (linear phase adjustment).
276 // As fftw_test.c says, “The paper performs more tests, but this code should be
277 // fine too”.
278 TEST(FFTPassEffectTest, ErgunShiftProperty) {
279         srand(1236);
280         const int max_fft_size = 64;
281         float a[max_fft_size * 4], b[max_fft_size * 4];
282         float a_out[max_fft_size * 4], b_out[max_fft_size * 4], expected_a_out[max_fft_size * 4];
283         for (int fft_size = 2; fft_size <= max_fft_size; fft_size *= 2) {
284                 for (int inverse = 0; inverse <= 1; ++inverse) {
285                         for (int direction = 0; direction <= 1; ++direction) {
286                                 for (int i = 0; i < ergun_rounds; ++i) {
287                                         for (int j = 0; j < fft_size * 4; ++j) {
288                                                 a[j] = uniform_random();
289                                         }
290
291                                         // Circular shift left by one step.
292                                         for (int j = 0; j < fft_size * 4; ++j) {
293                                                 b[j] = a[(j + 4) % (fft_size * 4)];
294                                         }
295                                         run_fft(a, a_out, fft_size, inverse, false, FFTPassEffect::Direction(direction));
296                                         run_fft(b, b_out, fft_size, inverse, false, FFTPassEffect::Direction(direction));
297
298                                         for (int j = 0; j < fft_size; ++j) {
299                                                 double s = -sin(j * 2.0 * M_PI / fft_size);
300                                                 double c = cos(j * 2.0 * M_PI / fft_size);
301                                                 if (inverse) {
302                                                         s = -s;
303                                                 }
304
305                                                 expected_a_out[j * 4 + 0] = b_out[j * 4 + 0] * c - b_out[j * 4 + 1] * s;
306                                                 expected_a_out[j * 4 + 1] = b_out[j * 4 + 0] * s + b_out[j * 4 + 1] * c;
307
308                                                 expected_a_out[j * 4 + 2] = b_out[j * 4 + 2] * c - b_out[j * 4 + 3] * s;
309                                                 expected_a_out[j * 4 + 3] = b_out[j * 4 + 2] * s + b_out[j * 4 + 3] * c;
310                                         }
311                                         expect_equal(expected_a_out, a_out, 4, fft_size, 0.025, 0.0005);
312                                 }
313                         }
314                 }
315         }
316 }
317
318 TEST(FFTPassEffectTest, BigFFTAccuracy) {
319         srand(1234);
320         const int max_fft_size = 2048;
321         float in[max_fft_size * 4], out[max_fft_size * 4], out2[max_fft_size * 4];
322         for (int fft_size = 2; fft_size <= max_fft_size; fft_size *= 2) {
323                 for (int j = 0; j < fft_size * 4; ++j) {
324                         in[j] = uniform_random();
325                 }
326                 run_fft(in, out, fft_size, false, true);  // Forward, with normalization.
327                 run_fft(out, out2, fft_size, true);       // Reverse.
328
329                 // These error bounds come from
330                 // http://en.wikipedia.org/wiki/Fast_Fourier_transform#Accuracy_and_approximations,
331                 // with empirically estimated epsilons. Note that the calculated
332                 // rms in expect_equal() is divided by sqrt(N), so we compensate
333                 // similarly here.
334                 double max_error = 0.0009 * log2(fft_size);
335                 double rms_limit = 0.0007 * sqrt(log2(fft_size)) / sqrt(fft_size);
336                 expect_equal(in, out2, 4, fft_size, max_error, rms_limit);
337         }
338 }
339
340 }  // namespace movit