third_party/aom/test/fft_test.cc

   1 /*
   2  * Copyright (c) 2018, Alliance for Open Media. All rights reserved
   3  *
   4  * This source code is subject to the terms of the BSD 2 Clause License and
   5  * the Alliance for Open Media Patent License 1.0. If the BSD 2 Clause License
   6  * was not distributed with this source code in the LICENSE file, you can
   7  * obtain it at www.aomedia.org/license/software. If the Alliance for Open
   8  * Media Patent License 1.0 was not distributed with this source code in the
   9  * PATENTS file, you can obtain it at www.aomedia.org/license/patent.
  10  */
  11
  12 #include <math.h>
  13
  14 #include <algorithm>
  15 #include <complex>
  16 #include <vector>
  17
  18 #include "aom_dsp/fft_common.h"
  19 #include "aom_mem/aom_mem.h"
  20 #include "av1/common/common.h"
  21 #include "config/aom_dsp_rtcd.h"
  22 #include "test/acm_random.h"
  23 #include "third_party/googletest/src/googletest/include/gtest/gtest.h"
  24
  25 namespace {
  26
  27 typedef void (*tform_fun_t)(const float *input, float *temp, float *output);
  28
  29 // Simple 1D FFT implementation
  30 template <typename InputType>
  31 void fft(const InputType *data, std::complex<float> *result, int n) {
  32   if (n == 1) {
  33     result[0] = data[0];
  34     return;
  35   }
  36   std::vector<InputType> temp(n);
  37   for (int k = 0; k < n / 2; ++k) {
  38     temp[k] = data[2 * k];
  39     temp[n / 2 + k] = data[2 * k + 1];
  40   }
  41   fft(&temp[0], result, n / 2);
  42   fft(&temp[n / 2], result + n / 2, n / 2);
  43   for (int k = 0; k < n / 2; ++k) {
  44     std::complex<float> w = std::complex<float>((float)cos(2. * PI * k / n),
  45                                                 (float)-sin(2. * PI * k / n));
  46     std::complex<float> a = result[k];
  47     std::complex<float> b = result[n / 2 + k];
  48     result[k] = a + w * b;
  49     result[n / 2 + k] = a - w * b;
  50   }
  51 }
  52
  53 void transpose(std::vector<std::complex<float> > *data, int n) {
  54   for (int y = 0; y < n; ++y) {
  55     for (int x = y + 1; x < n; ++x) {
  56       std::swap((*data)[y * n + x], (*data)[x * n + y]);
  57     }
  58   }
  59 }
  60
  61 // Simple 2D FFT implementation
  62 template <class InputType>
  63 std::vector<std::complex<float> > fft2d(const InputType *input, int n) {
  64   std::vector<std::complex<float> > rowfft(n * n);
  65   std::vector<std::complex<float> > result(n * n);
  66   for (int y = 0; y < n; ++y) {
  67     fft(input + y * n, &rowfft[y * n], n);
  68   }
  69   transpose(&rowfft, n);
  70   for (int y = 0; y < n; ++y) {
  71     fft(&rowfft[y * n], &result[y * n], n);
  72   }
  73   transpose(&result, n);
  74   return result;
  75 }
  76
  77 struct FFTTestArg {
  78   int n;
  79   void (*fft)(const float *input, float *temp, float *output);
  80   FFTTestArg(int n_in, tform_fun_t fft_in) : n(n_in), fft(fft_in) {}
  81 };
  82
  83 std::ostream &operator<<(std::ostream &os, const FFTTestArg &test_arg) {
  84   return os << "fft_arg { n:" << test_arg.n << " fft:" << test_arg.fft << " }";
  85 }
  86
  87 class FFT2DTest : public ::testing::TestWithParam<FFTTestArg> {
  88  protected:
  89   void SetUp() {
  90     int n = GetParam().n;
  91     input_ = (float *)aom_memalign(32, sizeof(*input_) * n * n);
  92     temp_ = (float *)aom_memalign(32, sizeof(*temp_) * n * n);
  93     output_ = (float *)aom_memalign(32, sizeof(*output_) * n * n * 2);
  94     memset(input_, 0, sizeof(*input_) * n * n);
  95     memset(temp_, 0, sizeof(*temp_) * n * n);
  96     memset(output_, 0, sizeof(*output_) * n * n * 2);
  97   }
  98   void TearDown() {
  99     aom_free(input_);
 100     aom_free(temp_);
 101     aom_free(output_);
 102   }
 103   float *input_;
 104   float *temp_;
 105   float *output_;
 106 };
 107
 108 TEST_P(FFT2DTest, Correct) {
 109   int n = GetParam().n;
 110   for (int i = 0; i < n * n; ++i) {
 111     input_[i] = 1;
 112     std::vector<std::complex<float> > expected = fft2d<float>(&input_[0], n);
 113     GetParam().fft(&input_[0], &temp_[0], &output_[0]);
 114     for (int y = 0; y < n; ++y) {
 115       for (int x = 0; x < (n / 2) + 1; ++x) {
 116         EXPECT_NEAR(expected[y * n + x].real(), output_[2 * (y * n + x)], 1e-5);
 117         EXPECT_NEAR(expected[y * n + x].imag(), output_[2 * (y * n + x) + 1],
 118                     1e-5);
 119       }
 120     }
 121     input_[i] = 0;
 122   }
 123 }
 124
 125 TEST_P(FFT2DTest, Benchmark) {
 126   int n = GetParam().n;
 127   float sum = 0;
 128   for (int i = 0; i < 1000 * (64 - n); ++i) {
 129     input_[i % (n * n)] = 1;
 130     GetParam().fft(&input_[0], &temp_[0], &output_[0]);
 131     sum += output_[0];
 132     input_[i % (n * n)] = 0;
 133   }
 134 }
 135
 136 INSTANTIATE_TEST_CASE_P(C, FFT2DTest,
 137                         ::testing::Values(FFTTestArg(2, aom_fft2x2_float_c),
 138                                           FFTTestArg(4, aom_fft4x4_float_c),
 139                                           FFTTestArg(8, aom_fft8x8_float_c),
 140                                           FFTTestArg(16, aom_fft16x16_float_c),
 141                                           FFTTestArg(32,
 142                                                      aom_fft32x32_float_c)));
 143 #if ARCH_X86 || ARCH_X86_64
 144 #if HAVE_SSE2
 145 INSTANTIATE_TEST_CASE_P(
 146     SSE2, FFT2DTest,
 147     ::testing::Values(FFTTestArg(4, aom_fft4x4_float_sse2),
 148                       FFTTestArg(8, aom_fft8x8_float_sse2),
 149                       FFTTestArg(16, aom_fft16x16_float_sse2),
 150                       FFTTestArg(32, aom_fft32x32_float_sse2)));
 151 #endif  // HAVE_SSE2
 152 #if HAVE_AVX2
 153 INSTANTIATE_TEST_CASE_P(
 154     AVX2, FFT2DTest,
 155     ::testing::Values(FFTTestArg(8, aom_fft8x8_float_avx2),
 156                       FFTTestArg(16, aom_fft16x16_float_avx2),
 157                       FFTTestArg(32, aom_fft32x32_float_avx2)));
 158 #endif  // HAVE_AVX2
 159 #endif  // ARCH_X86 || ARCH_X86_64
 160
 161 struct IFFTTestArg {
 162   int n;
 163   tform_fun_t ifft;
 164   IFFTTestArg(int n_in, tform_fun_t ifft_in) : n(n_in), ifft(ifft_in) {}
 165 };
 166
 167 std::ostream &operator<<(std::ostream &os, const IFFTTestArg &test_arg) {
 168   return os << "ifft_arg { n:" << test_arg.n << " fft:" << test_arg.ifft
 169             << " }";
 170 }
 171
 172 class IFFT2DTest : public ::testing::TestWithParam<IFFTTestArg> {
 173  protected:
 174   void SetUp() {
 175     int n = GetParam().n;
 176     input_ = (float *)aom_memalign(32, sizeof(*input_) * n * n * 2);
 177     temp_ = (float *)aom_memalign(32, sizeof(*temp_) * n * n * 2);
 178     output_ = (float *)aom_memalign(32, sizeof(*output_) * n * n);
 179     memset(input_, 0, sizeof(*input_) * n * n * 2);
 180     memset(temp_, 0, sizeof(*temp_) * n * n * 2);
 181     memset(output_, 0, sizeof(*output_) * n * n);
 182   }
 183   void TearDown() {
 184     aom_free(input_);
 185     aom_free(temp_);
 186     aom_free(output_);
 187   }
 188   float *input_;
 189   float *temp_;
 190   float *output_;
 191 };
 192
 193 TEST_P(IFFT2DTest, Correctness) {
 194   int n = GetParam().n;
 195   ASSERT_GE(n, 2);
 196   std::vector<float> expected(n * n);
 197   std::vector<float> actual(n * n);
 198   // Do forward transform then invert to make sure we get back expected
 199   for (int y = 0; y < n; ++y) {
 200     for (int x = 0; x < n; ++x) {
 201       expected[y * n + x] = 1;
 202       std::vector<std::complex<float> > input_c = fft2d(&expected[0], n);
 203       for (int i = 0; i < n * n; ++i) {
 204         input_[2 * i + 0] = input_c[i].real();
 205         input_[2 * i + 1] = input_c[i].imag();
 206       }
 207       GetParam().ifft(&input_[0], &temp_[0], &output_[0]);
 208
 209       for (int yy = 0; yy < n; ++yy) {
 210         for (int xx = 0; xx < n; ++xx) {
 211           EXPECT_NEAR(expected[yy * n + xx], output_[yy * n + xx] / (n * n),
 212                       1e-5);
 213         }
 214       }
 215       expected[y * n + x] = 0;
 216     }
 217   }
 218 };
 219
 220 TEST_P(IFFT2DTest, Benchmark) {
 221   int n = GetParam().n;
 222   float sum = 0;
 223   for (int i = 0; i < 1000 * (64 - n); ++i) {
 224     input_[i % (n * n)] = 1;
 225     GetParam().ifft(&input_[0], &temp_[0], &output_[0]);
 226     sum += output_[0];
 227     input_[i % (n * n)] = 0;
 228   }
 229 }
 230 INSTANTIATE_TEST_CASE_P(
 231     C, IFFT2DTest,
 232     ::testing::Values(IFFTTestArg(2, aom_ifft2x2_float_c),
 233                       IFFTTestArg(4, aom_ifft4x4_float_c),
 234                       IFFTTestArg(8, aom_ifft8x8_float_c),
 235                       IFFTTestArg(16, aom_ifft16x16_float_c),
 236                       IFFTTestArg(32, aom_ifft32x32_float_c)));
 237 #if ARCH_X86 || ARCH_X86_64
 238 #if HAVE_SSE2
 239 INSTANTIATE_TEST_CASE_P(
 240     SSE2, IFFT2DTest,
 241     ::testing::Values(IFFTTestArg(4, aom_ifft4x4_float_sse2),
 242                       IFFTTestArg(8, aom_ifft8x8_float_sse2),
 243                       IFFTTestArg(16, aom_ifft16x16_float_sse2),
 244                       IFFTTestArg(32, aom_ifft32x32_float_sse2)));
 245 #endif  // HAVE_SSE2
 246
 247 #if HAVE_AVX2
 248 INSTANTIATE_TEST_CASE_P(
 249     AVX2, IFFT2DTest,
 250     ::testing::Values(IFFTTestArg(8, aom_ifft8x8_float_avx2),
 251                       IFFTTestArg(16, aom_ifft16x16_float_avx2),
 252                       IFFTTestArg(32, aom_ifft32x32_float_avx2)));
 253 #endif  // HAVE_AVX2
 254 #endif  // ARCH_X86 || ARCH_X86_64
 255
 256 }  // namespace