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>
fft(const InputType * data,std::complex<float> * result,int n)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
transpose(std::vector<std::complex<float>> * data,int n)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>
fft2d(const InputType * input,int n)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);
FFTTestArg__anon75c8d0c20111::FFTTestArg80 FFTTestArg(int n_in, tform_fun_t fft_in) : n(n_in), fft(fft_in) {}
81 };
82
operator <<(std::ostream & os,const FFTTestArg & test_arg)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:
SetUp()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 }
TearDown()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
TEST_P(FFT2DTest,Correct)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
TEST_P(FFT2DTest,Benchmark)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;
IFFTTestArg__anon75c8d0c20111::IFFTTestArg164 IFFTTestArg(int n_in, tform_fun_t ifft_in) : n(n_in), ifft(ifft_in) {}
165 };
166
operator <<(std::ostream & os,const IFFTTestArg & test_arg)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:
SetUp()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 }
TearDown()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
TEST_P(IFFT2DTest,Correctness)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
TEST_P(IFFT2DTest,Benchmark)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
257