git.sesse.net Git - movit/blob - fft_pass_effect_test.cpp

   1 // Unit tests for FFTPassEffect.
   2
   3 #include <math.h>
   4 #include <stdlib.h>
   5 #include <string.h>
   6
   7 #include "effect_chain.h"
   8 #include "fft_pass_effect.h"
   9 #include "glew.h"
  10 #include "gtest/gtest.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         const int fft_size = 64;
 129         const int num_repeats = 31;  // Prime, to make things more challenging.
 130         float data[num_repeats * fft_size * 4] = { 0 };
 131         float expected_data[num_repeats * fft_size * 4], out_data[num_repeats * fft_size * 4];
 132
 133         srand(12345);
 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 TEST(FFTPassEffectTest, TwoDimensional) {  // Implicitly tests vertical.
 161         srand(1234);
 162         const int fft_size = 16;
 163         float in[fft_size * fft_size * 4], out[fft_size * fft_size * 4], expected_out[fft_size * fft_size * 4];
 164         for (int y = 0; y < fft_size; ++y) {
 165                 for (int x = 0; x < fft_size; ++x) {
 166                         in[(y * fft_size + x) * 4 + 0] =
 167                                 sin(2.0 * M_PI * (2 * x + 3 * y) / fft_size);
 168                         in[(y * fft_size + x) * 4 + 1] = 0.0;
 169                         in[(y * fft_size + x) * 4 + 2] = 0.0;
 170                         in[(y * fft_size + x) * 4 + 3] = 0.0;
 171                 }
 172         }
 173         memset(expected_out, 0, sizeof(expected_out));
 174
 175         // This result has been verified using the fft2() function in Octave,
 176         // which uses FFTW.
 177         expected_out[(3 * fft_size + 2) * 4 + 1] = -128.0;
 178         expected_out[(13 * fft_size + 14) * 4 + 1] = 128.0;
 179
 180         EffectChainTester tester(in, fft_size, fft_size, FORMAT_RGBA_PREMULTIPLIED_ALPHA, COLORSPACE_sRGB, GAMMA_LINEAR);
 181         setup_fft(tester.get_chain(), fft_size, false, false, FFTPassEffect::HORIZONTAL);
 182         setup_fft(tester.get_chain(), fft_size, false, false, FFTPassEffect::VERTICAL);
 183         tester.run(out, GL_RGBA, COLORSPACE_sRGB, GAMMA_LINEAR, OUTPUT_ALPHA_FORMAT_PREMULTIPLIED);
 184
 185         expect_equal(expected_out, out, 4 * fft_size, fft_size, 0.25, 0.0005);
 186 }
 187
 188 // The classic paper for FFT correctness testing is Funda Ergün:
 189 // “Testing Multivariate Linear Functions: Overcoming the Generator Bottleneck”
 190 // (http://www.cs.sfu.ca/~funda/PUBLICATIONS/stoc95.ps), which proves that
 191 // testing three basic properties of FFTs guarantees that the function is
 192 // correct (at least under the assumption that errors are random).
 193 //
 194 // We don't follow the paper directly, though, for a few reasons: First,
 195 // Ergün's paper really considers _self-correcting_ systems, which may
 196 // be stochastically faulty, and thus uses various relatively complicated
 197 // bounds and tests we don't really need. Second, the FFTs it considers
 198 // are all about polynomials over finite fields, which means that results
 199 // are exact and thus easy to test; we work with floats (half-floats!),
 200 // and thus need some error tolerance.
 201 //
 202 // So instead, we follow the implementation of FFTW, which is really the
 203 // gold standard when it comes to FFTs these days. They hard-code 20
 204 // testing rounds as opposed to the more complicated bounds in the paper,
 205 // and have a simpler version of the third test.
 206 //
 207 // The error bounds are set somewhat empirically, but remember that these
 208 // inputs will give frequency values as large as ~16, where 0.025 is
 209 // within the 9th bit (of 11 total mantissa bits in fp16).
 210 const int ergun_rounds = 20;
 211
 212 // Test 1: Test that FFT(a + b) = FFT(a) + FFT(b).
 213 TEST(FFTPassEffectTest, ErgunLinearityTest) {
 214         srand(1234);
 215         const int max_fft_size = 64;
 216         float a[max_fft_size * 4], b[max_fft_size * 4], sum[max_fft_size * 4];
 217         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];
 218         for (int fft_size = 2; fft_size <= max_fft_size; fft_size *= 2) {
 219                 for (int inverse = 0; inverse <= 1; ++inverse) {
 220                         for (int i = 0; i < ergun_rounds; ++i) {
 221                                 for (int j = 0; j < fft_size * 4; ++j) {
 222                                         a[j] = uniform_random();
 223                                         b[j] = uniform_random();
 224                                 }
 225                                 run_fft(a, a_out, fft_size, inverse);
 226                                 run_fft(b, b_out, fft_size, inverse);
 227
 228                                 for (int j = 0; j < fft_size * 4; ++j) {
 229                                         sum[j] = a[j] + b[j];
 230                                         expected_sum_out[j] = a_out[j] + b_out[j];
 231                                 }
 232
 233                                 run_fft(sum, sum_out, fft_size, inverse);
 234                                 expect_equal(expected_sum_out, sum_out, 4, fft_size, 0.03, 0.0005);
 235                         }
 236                 }
 237         }
 238 }
 239
 240 // Test 2: Test that FFT(delta(i)) = 1  (where delta(i) = [1 0 0 0 ...]),
 241 // or more specifically, test that FFT(a + delta(i)) - FFT(a) = 1.
 242 TEST(FFTPassEffectTest, ErgunImpulseTransform) {
 243         srand(1235);
 244         const int max_fft_size = 64;
 245         float a[max_fft_size * 4], b[max_fft_size * 4];
 246         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];
 247         for (int fft_size = 2; fft_size <= max_fft_size; fft_size *= 2) {
 248                 for (int inverse = 0; inverse <= 1; ++inverse) {
 249                         for (int i = 0; i < ergun_rounds; ++i) {
 250                                 for (int j = 0; j < fft_size * 4; ++j) {
 251                                         a[j] = uniform_random();
 252
 253                                         // Compute delta(j) - a.
 254                                         if (j < 4) {
 255                                                 b[j] = 1.0 - a[j];
 256                                         } else {
 257                                                 b[j] = -a[j];
 258                                         }
 259                                 }
 260                                 run_fft(a, a_out, fft_size, inverse);
 261                                 run_fft(b, b_out, fft_size, inverse);
 262
 263                                 for (int j = 0; j < fft_size * 4; ++j) {
 264                                         sum_out[j] = a_out[j] + b_out[j];
 265                                         expected_sum_out[j] = 1.0;
 266                                 }
 267                                 expect_equal(expected_sum_out, sum_out, 4, fft_size, 0.025, 0.0005);
 268                         }
 269                 }
 270         }
 271 }
 272
 273 // Test 3: Test the time-shift property of the FFT, in that a circular left-shift
 274 // multiplies the result by e^(j 2pi k/N) (linear phase adjustment).
 275 // As fftw_test.c says, “The paper performs more tests, but this code should be
 276 // fine too”.
 277 TEST(FFTPassEffectTest, ErgunShiftProperty) {
 278         srand(1236);
 279         const int max_fft_size = 64;
 280         float a[max_fft_size * 4], b[max_fft_size * 4];
 281         float a_out[max_fft_size * 4], b_out[max_fft_size * 4], expected_a_out[max_fft_size * 4];
 282         for (int fft_size = 2; fft_size <= max_fft_size; fft_size *= 2) {
 283                 for (int inverse = 0; inverse <= 1; ++inverse) {
 284                         for (int direction = 0; direction <= 1; ++direction) {
 285                                 for (int i = 0; i < ergun_rounds; ++i) {
 286                                         for (int j = 0; j < fft_size * 4; ++j) {
 287                                                 a[j] = uniform_random();
 288                                         }
 289
 290                                         // Circular shift left by one step.
 291                                         for (int j = 0; j < fft_size * 4; ++j) {
 292                                                 b[j] = a[(j + 4) % (fft_size * 4)];
 293                                         }
 294                                         run_fft(a, a_out, fft_size, inverse, false, FFTPassEffect::Direction(direction));
 295                                         run_fft(b, b_out, fft_size, inverse, false, FFTPassEffect::Direction(direction));
 296
 297                                         for (int j = 0; j < fft_size; ++j) {
 298                                                 double s = -sin(j * 2.0 * M_PI / fft_size);
 299                                                 double c = cos(j * 2.0 * M_PI / fft_size);
 300                                                 if (inverse) {
 301                                                         s = -s;
 302                                                 }
 303
 304                                                 expected_a_out[j * 4 + 0] = b_out[j * 4 + 0] * c - b_out[j * 4 + 1] * s;
 305                                                 expected_a_out[j * 4 + 1] = b_out[j * 4 + 0] * s + b_out[j * 4 + 1] * c;
 306
 307                                                 expected_a_out[j * 4 + 2] = b_out[j * 4 + 2] * c - b_out[j * 4 + 3] * s;
 308                                                 expected_a_out[j * 4 + 3] = b_out[j * 4 + 2] * s + b_out[j * 4 + 3] * c;
 309                                         }
 310                                         expect_equal(expected_a_out, a_out, 4, fft_size, 0.025, 0.0005);
 311                                 }
 312                         }
 313                 }
 314         }
 315 }
 316
 317 TEST(FFTPassEffectTest, BigFFTAccuracy) {
 318         srand(1234);
 319         const int max_fft_size = 2048;
 320         float in[max_fft_size * 4], out[max_fft_size * 4], out2[max_fft_size * 4];
 321         for (int fft_size = 2; fft_size <= max_fft_size; fft_size *= 2) {
 322                 for (int j = 0; j < fft_size * 4; ++j) {
 323                         in[j] = uniform_random();
 324                 }
 325                 run_fft(in, out, fft_size, false, true);  // Forward, with normalization.
 326                 run_fft(out, out2, fft_size, true);       // Reverse.
 327
 328                 // These error bounds come from
 329                 // http://en.wikipedia.org/wiki/Fast_Fourier_transform#Accuracy_and_approximations,
 330                 // with empirically estimated epsilons. Note that the calculated
 331                 // rms in expect_equal() is divided by sqrt(N), so we compensate
 332                 // similarly here.
 333                 double max_error = 0.0009 * log2(fft_size);
 334                 double rms_limit = 0.0007 * sqrt(log2(fft_size)) / sqrt(fft_size);
 335                 expect_equal(in, out2, 4, fft_size, max_error, rms_limit);
 336         }
 337 }
 338
 339 }  // namespace movit