// Copyright 2022 Google LLC // // This source code is licensed under the BSD-style license found in the // LICENSE file in the root directory of this source tree. #pragma once #include #include #include #include #include #include #include #include #include #include #include #include // twiddle table for bfly4 for fft size 256 (complex numbers) // Even numbers are numpy.floor(0.5 + 32767 * numpy.cos(-2*pi*numpy.linspace(0, 255, num=256) / 256)).astype(numpy.int16).tolist() // Odd numbers are numpy.floor(0.5 + 32767 * numpy.sin(-2*pi*numpy.linspace(0, 255, num=256) / 256)).astype(numpy.int16).tolist() static const int16_t xnn_reference_table_fft256_twiddle[512] = { 32767, 0, 32757, -804, 32728, -1608, 32678, -2410, 32609, -3212, 32521, -4011, 32412, -4808, 32285, -5602, 32137, -6393, 31971, -7179, 31785, -7962, 31580, -8739, 31356, -9512, 31113,-10278, 30852,-11039, 30571,-11793, 30273,-12539, 29956,-13279, 29621,-14010, 29268,-14732, 28898,-15446, 28510,-16151, 28105,-16846, 27683,-17530, 27245,-18204, 26790,-18868, 26319,-19519, 25832,-20159, 25329,-20787, 24811,-21403, 24279,-22005, 23731,-22594, 23170,-23170, 22594,-23731, 22005,-24279, 21403,-24811, 20787,-25329, 20159,-25832, 19519,-26319, 18868,-26790, 18204,-27245, 17530,-27683, 16846,-28105, 16151,-28510, 15446,-28898, 14732,-29268, 14010,-29621, 13279,-29956, 12539,-30273, 11793,-30571, 11039,-30852, 10278,-31113, 9512,-31356, 8739,-31580, 7962,-31785, 7179,-31971, 6393,-32137, 5602,-32285, 4808,-32412, 4011,-32521, 3212,-32609, 2410,-32678, 1608,-32728, 804,-32757, 0,-32767, -804,-32757, -1608,-32728, -2410,-32678, -3212,-32609, -4011,-32521, -4808,-32412, -5602,-32285, -6393,-32137, -7179,-31971, -7962,-31785, -8739,-31580, -9512,-31356, -10278,-31113, -11039,-30852, -11793,-30571, -12539,-30273, -13279,-29956, -14010,-29621, -14732,-29268, -15446,-28898, -16151,-28510, -16846,-28105, -17530,-27683, -18204,-27245, -18868,-26790, -19519,-26319, -20159,-25832, -20787,-25329, -21403,-24811, -22005,-24279, -22594,-23731, -23170,-23170, -23731,-22594, -24279,-22005, -24811,-21403, -25329,-20787, -25832,-20159, -26319,-19519, -26790,-18868, -27245,-18204, -27683,-17530, -28105,-16846, -28510,-16151, -28898,-15446, -29268,-14732, -29621,-14010, -29956,-13279, -30273,-12539, -30571,-11793, -30852,-11039, -31113,-10278, -31356, -9512, -31580, -8739, -31785, -7962, -31971, -7179, -32137, -6393, -32285, -5602, -32412, -4808, -32521, -4011, -32609, -3212, -32678, -2410, -32728, -1608, -32757, -804, -32767, 0, -32757, 804, -32728, 1608, -32678, 2410, -32609, 3212, -32521, 4011, -32412, 4808, -32285, 5602, -32137, 6393, -31971, 7179, -31785, 7962, -31580, 8739, -31356, 9512, -31113, 10278, -30852, 11039, -30571, 11793, -30273, 12539, -29956, 13279, -29621, 14010, -29268, 14732, -28898, 15446, -28510, 16151, -28105, 16846, -27683, 17530, -27245, 18204, -26790, 18868, -26319, 19519, -25832, 20159, -25329, 20787, -24811, 21403, -24279, 22005, -23731, 22594, -23170, 23170, -22594, 23731, -22005, 24279, -21403, 24811, -20787, 25329, -20159, 25832, -19519, 26319, -18868, 26790, -18204, 27245, -17530, 27683, -16846, 28105, -16151, 28510, -15446, 28898, -14732, 29268, -14010, 29621, -13279, 29956, -12539, 30273, -11793, 30571, -11039, 30852, -10278, 31113, -9512, 31356, -8739, 31580, -7962, 31785, -7179, 31971, -6393, 32137, -5602, 32285, -4808, 32412, -4011, 32521, -3212, 32609, -2410, 32678, -1608, 32728, -804, 32757, 0, 32767, 804, 32757, 1608, 32728, 2410, 32678, 3212, 32609, 4011, 32521, 4808, 32412, 5602, 32285, 6393, 32137, 7179, 31971, 7962, 31785, 8739, 31580, 9512, 31356, 10278, 31113, 11039, 30852, 11793, 30571, 12539, 30273, 13279, 29956, 14010, 29621, 14732, 29268, 15446, 28898, 16151, 28510, 16846, 28105, 17530, 27683, 18204, 27245, 18868, 26790, 19519, 26319, 20159, 25832, 20787, 25329, 21403, 24811, 22005, 24279, 22594, 23731, 23170, 23170, 23731, 22594, 24279, 22005, 24811, 21403, 25329, 20787, 25832, 20159, 26319, 19519, 26790, 18868, 27245, 18204, 27683, 17530, 28105, 16846, 28510, 16151, 28898, 15446, 29268, 14732, 29621, 14010, 29956, 13279, 30273, 12539, 30571, 11793, 30852, 11039, 31113, 10278, 31356, 9512, 31580, 8739, 31785, 7962, 31971, 7179, 32137, 6393, 32285, 5602, 32412, 4808, 32521, 4011, 32609, 3212, 32678, 2410, 32728, 1608, 32757, 804 }; void xnn_cs16_bfly4_reference( size_t samples, int16_t* data, const size_t stride, const int16_t* twiddle) { const int16_t* tw1 = twiddle; const int16_t* tw2 = tw1; const int16_t* tw3 = tw1; int16_t* out0 = data; int16_t* out1 = data + samples * 2; int16_t* out2 = data + samples * 4; int16_t* out3 = data + samples * 6; assert(samples != 0); assert(stride != 0); assert(twiddle != NULL); assert(data != NULL); do { int32_t vout0_r = (int32_t) out0[0]; int32_t vout0_i = (int32_t) out0[1]; int32_t vout1_r = (int32_t) out1[0]; int32_t vout1_i = (int32_t) out1[1]; int32_t vout2_r = (int32_t) out2[0]; int32_t vout2_i = (int32_t) out2[1]; int32_t vout3_r = (int32_t) out3[0]; int32_t vout3_i = (int32_t) out3[1]; const int32_t tw1_r = (const int32_t) tw1[0]; const int32_t tw1_i = (const int32_t) tw1[1]; const int32_t tw2_r = (const int32_t) tw2[0]; const int32_t tw2_i = (const int32_t) tw2[1]; const int32_t tw3_r = (const int32_t) tw3[0]; const int32_t tw3_i = (const int32_t) tw3[1]; // Note 32767 / 4 = 8191. Should be 8192. vout0_r = (vout0_r * 8191 + 16384) >> 15; vout0_i = (vout0_i * 8191 + 16384) >> 15; vout1_r = (vout1_r * 8191 + 16384) >> 15; vout1_i = (vout1_i * 8191 + 16384) >> 15; vout2_r = (vout2_r * 8191 + 16384) >> 15; vout2_i = (vout2_i * 8191 + 16384) >> 15; vout3_r = (vout3_r * 8191 + 16384) >> 15; vout3_i = (vout3_i * 8191 + 16384) >> 15; const int32_t vtmp0_r = math_asr_s32(vout1_r * tw1_r - vout1_i * tw1_i + 16384, 15); const int32_t vtmp0_i = math_asr_s32(vout1_r * tw1_i + vout1_i * tw1_r + 16384, 15); const int32_t vtmp1_r = math_asr_s32(vout2_r * tw2_r - vout2_i * tw2_i + 16384, 15); const int32_t vtmp1_i = math_asr_s32(vout2_r * tw2_i + vout2_i * tw2_r + 16384, 15); const int32_t vtmp2_r = math_asr_s32(vout3_r * tw3_r - vout3_i * tw3_i + 16384, 15); const int32_t vtmp2_i = math_asr_s32(vout3_r * tw3_i + vout3_i * tw3_r + 16384, 15); const int32_t vtmp5_r = vout0_r - vtmp1_r; const int32_t vtmp5_i = vout0_i - vtmp1_i; vout0_r += vtmp1_r; vout0_i += vtmp1_i; const int32_t vtmp3_r = vtmp0_r + vtmp2_r; const int32_t vtmp3_i = vtmp0_i + vtmp2_i; const int32_t vtmp4_r = vtmp0_r - vtmp2_r; const int32_t vtmp4_i = vtmp0_i - vtmp2_i; vout2_r = vout0_r - vtmp3_r; vout2_i = vout0_i - vtmp3_i; tw1 += stride * 2; tw2 += stride * 4; tw3 += stride * 6; vout0_r += vtmp3_r; vout0_i += vtmp3_i; vout1_r = vtmp5_r + vtmp4_i; vout1_i = vtmp5_i - vtmp4_r; vout3_r = vtmp5_r - vtmp4_i; vout3_i = vtmp5_i + vtmp4_r; out0[0] = (int16_t) vout0_r; out0[1] = (int16_t) vout0_i; out1[0] = (int16_t) vout1_r; out1[1] = (int16_t) vout1_i; out2[0] = (int16_t) vout2_r; out2[1] = (int16_t) vout2_i; out3[0] = (int16_t) vout3_r; out3[1] = (int16_t) vout3_i; out0 += 2; out1 += 2; out2 += 2; out3 += 2; } while(--samples != 0); } class BFly4MicrokernelTester { public: inline BFly4MicrokernelTester& samples(size_t samples) { assert(samples != 0); this->samples_ = samples; return *this; } inline size_t samples() const { return this->samples_; } inline BFly4MicrokernelTester& stride(uint32_t stride) { this->stride_ = stride; return *this; } inline uint32_t stride() const { return this->stride_; } inline BFly4MicrokernelTester& iterations(size_t iterations) { this->iterations_ = iterations; return *this; } inline size_t iterations() const { return this->iterations_; } void Test(xnn_cs16_bfly4_ukernel_function bfly4) const { std::random_device random_device; auto rng = std::mt19937(random_device()); auto i16rng = std::bind(std::uniform_int_distribution(), std::ref(rng)); const size_t fft_size = (samples() == 1 ? 1 : (samples() * stride())) * 4; // 4 for bfly4. // 256 complex numbers = fft_size * 2 = 512 std::vector y(fft_size * 2); std::vector y_ref(fft_size * 2); for (size_t iteration = 0; iteration < iterations(); iteration++) { std::generate(y.begin(), y.end(), std::ref(i16rng)); y_ref = y; // Compute reference results. xnn_cs16_bfly4_reference(samples(), y_ref.data(), stride(), xnn_reference_table_fft256_twiddle); // Call optimized micro-kernel. bfly4(samples(), y.data(), stride(), xnn_reference_table_fft256_twiddle); // Verify results. for (size_t n = 0; n < fft_size * 2; n++) { ASSERT_EQ(y[n], y_ref[n]) << "at sample " << n << " / " << fft_size << "\nsamples " << samples() << "\nstride " << stride(); } } } private: size_t samples_{1}; uint32_t stride_{1}; size_t iterations_{15}; };