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