1 // Unit tests for FFTPassEffect.
7 #include "effect_chain.h"
8 #include "fft_pass_effect.h"
10 #include "gtest/gtest.h"
11 #include "image_format.h"
12 #include "multiply_effect.h"
13 #include "test_util.h"
17 // Generate a random number uniformly distributed between [-1.0, 1.0].
18 float uniform_random()
20 return 2.0 * ((float)rand() / RAND_MAX - 0.5);
23 void setup_fft(EffectChain *chain, int fft_size, bool inverse,
24 bool add_normalizer = false,
25 FFTPassEffect::Direction direction = FFTPassEffect::HORIZONTAL)
27 assert((fft_size & (fft_size - 1)) == 0); // Must be power of two.
28 for (int i = 1, subsize = 2; subsize <= fft_size; ++i, subsize *= 2) {
29 Effect *fft_effect = chain->add_effect(new FFTPassEffect());
30 bool ok = fft_effect->set_int("fft_size", fft_size);
31 ok |= fft_effect->set_int("pass_number", i);
32 ok |= fft_effect->set_int("inverse", inverse);
33 ok |= fft_effect->set_int("direction", direction);
38 float factor[4] = { 1.0f / fft_size, 1.0f / fft_size, 1.0f / fft_size, 1.0f / fft_size };
39 Effect *multiply_effect = chain->add_effect(new MultiplyEffect());
40 bool ok = multiply_effect->set_vec4("factor", factor);
45 void run_fft(const float *in, float *out, int fft_size, bool inverse,
46 bool add_normalizer = false,
47 FFTPassEffect::Direction direction = FFTPassEffect::HORIZONTAL)
50 if (direction == FFTPassEffect::HORIZONTAL) {
57 EffectChainTester tester(in, width, height, FORMAT_RGBA_PREMULTIPLIED_ALPHA, COLORSPACE_sRGB, GAMMA_LINEAR);
58 setup_fft(tester.get_chain(), fft_size, inverse, add_normalizer, direction);
59 tester.run(out, GL_RGBA, COLORSPACE_sRGB, GAMMA_LINEAR, OUTPUT_ALPHA_FORMAT_PREMULTIPLIED);
64 TEST(FFTPassEffectTest, ZeroStaysZero) {
65 const int fft_size = 64;
66 float data[fft_size * 4] = { 0 };
67 float out_data[fft_size * 4];
69 run_fft(data, out_data, fft_size, false);
70 expect_equal(data, out_data, 4, fft_size);
72 run_fft(data, out_data, fft_size, true);
73 expect_equal(data, out_data, 4, fft_size);
76 TEST(FFTPassEffectTest, Impulse) {
77 const int fft_size = 64;
78 float data[fft_size * 4] = { 0 };
79 float expected_data[fft_size * 4], out_data[fft_size * 4];
85 for (int i = 0; i < fft_size; ++i) {
86 expected_data[i * 4 + 0] = data[0];
87 expected_data[i * 4 + 1] = data[1];
88 expected_data[i * 4 + 2] = data[2];
89 expected_data[i * 4 + 3] = data[3];
92 run_fft(data, out_data, fft_size, false);
93 expect_equal(expected_data, out_data, 4, fft_size);
95 run_fft(data, out_data, fft_size, true);
96 expect_equal(expected_data, out_data, 4, fft_size);
99 TEST(FFTPassEffectTest, SingleFrequency) {
100 const int fft_size = 16;
101 float data[fft_size * 4] = { 0 };
102 float expected_data[fft_size * 4], out_data[fft_size * 4];
103 for (int i = 0; i < fft_size; ++i) {
104 data[i * 4 + 0] = sin(2.0 * M_PI * (4.0 * i) / fft_size);
105 data[i * 4 + 1] = 0.0;
106 data[i * 4 + 2] = 0.0;
107 data[i * 4 + 3] = 0.0;
109 for (int i = 0; i < fft_size; ++i) {
110 expected_data[i * 4 + 0] = 0.0;
111 expected_data[i * 4 + 1] = 0.0;
112 expected_data[i * 4 + 2] = 0.0;
113 expected_data[i * 4 + 3] = 0.0;
115 expected_data[4 * 4 + 1] = -8.0;
116 expected_data[12 * 4 + 1] = 8.0;
118 run_fft(data, out_data, fft_size, false, false, FFTPassEffect::HORIZONTAL);
119 expect_equal(expected_data, out_data, 4, fft_size);
121 run_fft(data, out_data, fft_size, false, false, FFTPassEffect::VERTICAL);
122 expect_equal(expected_data, out_data, 4, fft_size);
125 TEST(FFTPassEffectTest, Repeat) {
126 const int fft_size = 64;
127 const int num_repeats = 31; // Prime, to make things more challenging.
128 float data[num_repeats * fft_size * 4] = { 0 };
129 float expected_data[num_repeats * fft_size * 4], out_data[num_repeats * fft_size * 4];
132 for (int i = 0; i < num_repeats * fft_size * 4; ++i) {
133 data[i] = uniform_random();
136 for (int i = 0; i < num_repeats; ++i) {
137 run_fft(data + i * fft_size * 4, expected_data + i * fft_size * 4, fft_size, false);
142 EffectChainTester tester(data, num_repeats * fft_size, 1, FORMAT_RGBA_PREMULTIPLIED_ALPHA, COLORSPACE_sRGB, GAMMA_LINEAR);
143 setup_fft(tester.get_chain(), fft_size, false);
144 tester.run(out_data, GL_RGBA, COLORSPACE_sRGB, GAMMA_LINEAR, OUTPUT_ALPHA_FORMAT_PREMULTIPLIED);
146 expect_equal(expected_data, out_data, 4, num_repeats * fft_size);
150 EffectChainTester tester(data, 1, num_repeats * fft_size, FORMAT_RGBA_PREMULTIPLIED_ALPHA, COLORSPACE_sRGB, GAMMA_LINEAR);
151 setup_fft(tester.get_chain(), fft_size, false, false, FFTPassEffect::VERTICAL);
152 tester.run(out_data, GL_RGBA, COLORSPACE_sRGB, GAMMA_LINEAR, OUTPUT_ALPHA_FORMAT_PREMULTIPLIED);
154 expect_equal(expected_data, out_data, 4, num_repeats * fft_size);
158 TEST(FFTPassEffectTest, TwoDimensional) { // Implicitly tests vertical.
160 const int fft_size = 16;
161 float in[fft_size * fft_size * 4], out[fft_size * fft_size * 4], expected_out[fft_size * fft_size * 4];
162 for (int y = 0; y < fft_size; ++y) {
163 for (int x = 0; x < fft_size; ++x) {
164 in[(y * fft_size + x) * 4 + 0] =
165 sin(2.0 * M_PI * (2 * x + 3 * y) / fft_size);
166 in[(y * fft_size + x) * 4 + 1] = 0.0;
167 in[(y * fft_size + x) * 4 + 2] = 0.0;
168 in[(y * fft_size + x) * 4 + 3] = 0.0;
171 memset(expected_out, 0, sizeof(expected_out));
173 // This result has been verified using the fft2() function in Octave,
175 expected_out[(3 * fft_size + 2) * 4 + 1] = -128.0;
176 expected_out[(13 * fft_size + 14) * 4 + 1] = 128.0;
178 EffectChainTester tester(in, fft_size, fft_size, FORMAT_RGBA_PREMULTIPLIED_ALPHA, COLORSPACE_sRGB, GAMMA_LINEAR);
179 setup_fft(tester.get_chain(), fft_size, false, false, FFTPassEffect::HORIZONTAL);
180 setup_fft(tester.get_chain(), fft_size, false, false, FFTPassEffect::VERTICAL);
181 tester.run(out, GL_RGBA, COLORSPACE_sRGB, GAMMA_LINEAR, OUTPUT_ALPHA_FORMAT_PREMULTIPLIED);
183 expect_equal(expected_out, out, 4 * fft_size, fft_size, 0.25, 0.0005);
186 // The classic paper for FFT correctness testing is Funda Ergün:
187 // “Testing Multivariate Linear Functions: Overcoming the Generator Bottleneck”
188 // (http://www.cs.sfu.ca/~funda/PUBLICATIONS/stoc95.ps), which proves that
189 // testing three basic properties of FFTs guarantees that the function is
190 // correct (at least under the assumption that errors are random).
192 // We don't follow the paper directly, though, for a few reasons: First,
193 // Ergün's paper really considers _self-correcting_ systems, which may
194 // be stochastically faulty, and thus uses various relatively complicated
195 // bounds and tests we don't really need. Second, the FFTs it considers
196 // are all about polynomials over finite fields, which means that results
197 // are exact and thus easy to test; we work with floats (half-floats!),
198 // and thus need some error tolerance.
200 // So instead, we follow the implementation of FFTW, which is really the
201 // gold standard when it comes to FFTs these days. They hard-code 20
202 // testing rounds as opposed to the more complicated bounds in the paper,
203 // and have a simpler version of the third test.
205 // The error bounds are set somewhat empirically, but remember that these
206 // inputs will give frequency values as large as ~16, where 0.025 is
207 // within the 9th bit (of 11 total mantissa bits in fp16).
208 const int ergun_rounds = 20;
210 // Test 1: Test that FFT(a + b) = FFT(a) + FFT(b).
211 TEST(FFTPassEffectTest, ErgunLinearityTest) {
213 const int max_fft_size = 64;
214 float a[max_fft_size * 4], b[max_fft_size * 4], sum[max_fft_size * 4];
215 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];
216 for (int fft_size = 2; fft_size <= max_fft_size; fft_size *= 2) {
217 for (int inverse = 0; inverse <= 1; ++inverse) {
218 for (int i = 0; i < ergun_rounds; ++i) {
219 for (int j = 0; j < fft_size * 4; ++j) {
220 a[j] = uniform_random();
221 b[j] = uniform_random();
223 run_fft(a, a_out, fft_size, inverse);
224 run_fft(b, b_out, fft_size, inverse);
226 for (int j = 0; j < fft_size * 4; ++j) {
227 sum[j] = a[j] + b[j];
228 expected_sum_out[j] = a_out[j] + b_out[j];
231 run_fft(sum, sum_out, fft_size, inverse);
232 expect_equal(expected_sum_out, sum_out, 4, fft_size, 0.03, 0.0005);
238 // Test 2: Test that FFT(delta(i)) = 1 (where delta(i) = [1 0 0 0 ...]),
239 // or more specifically, test that FFT(a + delta(i)) - FFT(a) = 1.
240 TEST(FFTPassEffectTest, ErgunImpulseTransform) {
242 const int max_fft_size = 64;
243 float a[max_fft_size * 4], b[max_fft_size * 4];
244 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];
245 for (int fft_size = 2; fft_size <= max_fft_size; fft_size *= 2) {
246 for (int inverse = 0; inverse <= 1; ++inverse) {
247 for (int i = 0; i < ergun_rounds; ++i) {
248 for (int j = 0; j < fft_size * 4; ++j) {
249 a[j] = uniform_random();
251 // Compute delta(j) - a.
258 run_fft(a, a_out, fft_size, inverse);
259 run_fft(b, b_out, fft_size, inverse);
261 for (int j = 0; j < fft_size * 4; ++j) {
262 sum_out[j] = a_out[j] + b_out[j];
263 expected_sum_out[j] = 1.0;
265 expect_equal(expected_sum_out, sum_out, 4, fft_size, 0.025, 0.0005);
271 // Test 3: Test the time-shift property of the FFT, in that a circular left-shift
272 // multiplies the result by e^(j 2pi k/N) (linear phase adjustment).
273 // As fftw_test.c says, “The paper performs more tests, but this code should be
275 TEST(FFTPassEffectTest, ErgunShiftProperty) {
277 const int max_fft_size = 64;
278 float a[max_fft_size * 4], b[max_fft_size * 4];
279 float a_out[max_fft_size * 4], b_out[max_fft_size * 4], expected_a_out[max_fft_size * 4];
280 for (int fft_size = 2; fft_size <= max_fft_size; fft_size *= 2) {
281 for (int inverse = 0; inverse <= 1; ++inverse) {
282 for (int direction = 0; direction <= 1; ++direction) {
283 for (int i = 0; i < ergun_rounds; ++i) {
284 for (int j = 0; j < fft_size * 4; ++j) {
285 a[j] = uniform_random();
288 // Circular shift left by one step.
289 for (int j = 0; j < fft_size * 4; ++j) {
290 b[j] = a[(j + 4) % (fft_size * 4)];
292 run_fft(a, a_out, fft_size, inverse, false, FFTPassEffect::Direction(direction));
293 run_fft(b, b_out, fft_size, inverse, false, FFTPassEffect::Direction(direction));
295 for (int j = 0; j < fft_size; ++j) {
296 double s = -sin(j * 2.0 * M_PI / fft_size);
297 double c = cos(j * 2.0 * M_PI / fft_size);
302 expected_a_out[j * 4 + 0] = b_out[j * 4 + 0] * c - b_out[j * 4 + 1] * s;
303 expected_a_out[j * 4 + 1] = b_out[j * 4 + 0] * s + b_out[j * 4 + 1] * c;
305 expected_a_out[j * 4 + 2] = b_out[j * 4 + 2] * c - b_out[j * 4 + 3] * s;
306 expected_a_out[j * 4 + 3] = b_out[j * 4 + 2] * s + b_out[j * 4 + 3] * c;
308 expect_equal(expected_a_out, a_out, 4, fft_size, 0.025, 0.0005);
315 TEST(FFTPassEffectTest, BigFFTAccuracy) {
317 const int max_fft_size = 2048;
318 float in[max_fft_size * 4], out[max_fft_size * 4], out2[max_fft_size * 4];
319 for (int fft_size = 2; fft_size <= max_fft_size; fft_size *= 2) {
320 for (int j = 0; j < fft_size * 4; ++j) {
321 in[j] = uniform_random();
323 run_fft(in, out, fft_size, false, true); // Forward, with normalization.
324 run_fft(out, out2, fft_size, true); // Reverse.
326 // These error bounds come from
327 // http://en.wikipedia.org/wiki/Fast_Fourier_transform#Accuracy_and_approximations,
328 // with empirically estimated epsilons. Note that the calculated
329 // rms in expect_equal() is divided by sqrt(N), so we compensate
331 double max_error = 0.0009 * log2(fft_size);
332 double rms_limit = 0.0007 * sqrt(log2(fft_size)) / sqrt(fft_size);
333 expect_equal(in, out2, 4, fft_size, max_error, rms_limit);