1 /*
2 * Copyright (c) 2014 The WebRTC project authors. All Rights Reserved.
3 *
4 * Use of this source code is governed by a BSD-style license
5 * that can be found in the LICENSE file in the root of the source
6 * tree. An additional intellectual property rights grant can be found
7 * in the file PATENTS. All contributing project authors may
8 * be found in the AUTHORS file in the root of the source tree.
9 */
10
11 /*
12 * The rdft AEC algorithm, neon version of speed-critical functions.
13 *
14 * Based on the sse2 version.
15 */
16
17 #include <arm_neon.h>
18
19 #include "common_audio/third_party/ooura/fft_size_128/ooura_fft.h"
20 #include "common_audio/third_party/ooura/fft_size_128/ooura_fft_tables_common.h"
21 #include "common_audio/third_party/ooura/fft_size_128/ooura_fft_tables_neon_sse2.h"
22
23 namespace webrtc {
24
25 #if defined(WEBRTC_HAS_NEON)
cft1st_128_neon(float * a)26 void cft1st_128_neon(float* a) {
27 const float32x4_t vec_swap_sign = vld1q_f32((float32_t*)k_swap_sign);
28 int j, k2;
29
30 for (k2 = 0, j = 0; j < 128; j += 16, k2 += 4) {
31 float32x4_t a00v = vld1q_f32(&a[j + 0]);
32 float32x4_t a04v = vld1q_f32(&a[j + 4]);
33 float32x4_t a08v = vld1q_f32(&a[j + 8]);
34 float32x4_t a12v = vld1q_f32(&a[j + 12]);
35 float32x4_t a01v = vcombine_f32(vget_low_f32(a00v), vget_low_f32(a08v));
36 float32x4_t a23v = vcombine_f32(vget_high_f32(a00v), vget_high_f32(a08v));
37 float32x4_t a45v = vcombine_f32(vget_low_f32(a04v), vget_low_f32(a12v));
38 float32x4_t a67v = vcombine_f32(vget_high_f32(a04v), vget_high_f32(a12v));
39 const float32x4_t wk1rv = vld1q_f32(&rdft_wk1r[k2]);
40 const float32x4_t wk1iv = vld1q_f32(&rdft_wk1i[k2]);
41 const float32x4_t wk2rv = vld1q_f32(&rdft_wk2r[k2]);
42 const float32x4_t wk2iv = vld1q_f32(&rdft_wk2i[k2]);
43 const float32x4_t wk3rv = vld1q_f32(&rdft_wk3r[k2]);
44 const float32x4_t wk3iv = vld1q_f32(&rdft_wk3i[k2]);
45 float32x4_t x0v = vaddq_f32(a01v, a23v);
46 const float32x4_t x1v = vsubq_f32(a01v, a23v);
47 const float32x4_t x2v = vaddq_f32(a45v, a67v);
48 const float32x4_t x3v = vsubq_f32(a45v, a67v);
49 const float32x4_t x3w = vrev64q_f32(x3v);
50 float32x4_t x0w;
51 a01v = vaddq_f32(x0v, x2v);
52 x0v = vsubq_f32(x0v, x2v);
53 x0w = vrev64q_f32(x0v);
54 a45v = vmulq_f32(wk2rv, x0v);
55 a45v = vmlaq_f32(a45v, wk2iv, x0w);
56 x0v = vmlaq_f32(x1v, x3w, vec_swap_sign);
57 x0w = vrev64q_f32(x0v);
58 a23v = vmulq_f32(wk1rv, x0v);
59 a23v = vmlaq_f32(a23v, wk1iv, x0w);
60 x0v = vmlsq_f32(x1v, x3w, vec_swap_sign);
61 x0w = vrev64q_f32(x0v);
62 a67v = vmulq_f32(wk3rv, x0v);
63 a67v = vmlaq_f32(a67v, wk3iv, x0w);
64 a00v = vcombine_f32(vget_low_f32(a01v), vget_low_f32(a23v));
65 a04v = vcombine_f32(vget_low_f32(a45v), vget_low_f32(a67v));
66 a08v = vcombine_f32(vget_high_f32(a01v), vget_high_f32(a23v));
67 a12v = vcombine_f32(vget_high_f32(a45v), vget_high_f32(a67v));
68 vst1q_f32(&a[j + 0], a00v);
69 vst1q_f32(&a[j + 4], a04v);
70 vst1q_f32(&a[j + 8], a08v);
71 vst1q_f32(&a[j + 12], a12v);
72 }
73 }
74
cftmdl_128_neon(float * a)75 void cftmdl_128_neon(float* a) {
76 int j;
77 const int l = 8;
78 const float32x4_t vec_swap_sign = vld1q_f32((float32_t*)k_swap_sign);
79 float32x4_t wk1rv = vld1q_f32(cftmdl_wk1r);
80
81 for (j = 0; j < l; j += 2) {
82 const float32x2_t a_00 = vld1_f32(&a[j + 0]);
83 const float32x2_t a_08 = vld1_f32(&a[j + 8]);
84 const float32x2_t a_32 = vld1_f32(&a[j + 32]);
85 const float32x2_t a_40 = vld1_f32(&a[j + 40]);
86 const float32x4_t a_00_32 = vcombine_f32(a_00, a_32);
87 const float32x4_t a_08_40 = vcombine_f32(a_08, a_40);
88 const float32x4_t x0r0_0i0_0r1_x0i1 = vaddq_f32(a_00_32, a_08_40);
89 const float32x4_t x1r0_1i0_1r1_x1i1 = vsubq_f32(a_00_32, a_08_40);
90 const float32x2_t a_16 = vld1_f32(&a[j + 16]);
91 const float32x2_t a_24 = vld1_f32(&a[j + 24]);
92 const float32x2_t a_48 = vld1_f32(&a[j + 48]);
93 const float32x2_t a_56 = vld1_f32(&a[j + 56]);
94 const float32x4_t a_16_48 = vcombine_f32(a_16, a_48);
95 const float32x4_t a_24_56 = vcombine_f32(a_24, a_56);
96 const float32x4_t x2r0_2i0_2r1_x2i1 = vaddq_f32(a_16_48, a_24_56);
97 const float32x4_t x3r0_3i0_3r1_x3i1 = vsubq_f32(a_16_48, a_24_56);
98 const float32x4_t xx0 = vaddq_f32(x0r0_0i0_0r1_x0i1, x2r0_2i0_2r1_x2i1);
99 const float32x4_t xx1 = vsubq_f32(x0r0_0i0_0r1_x0i1, x2r0_2i0_2r1_x2i1);
100 const float32x4_t x3i0_3r0_3i1_x3r1 = vrev64q_f32(x3r0_3i0_3r1_x3i1);
101 const float32x4_t x1_x3_add =
102 vmlaq_f32(x1r0_1i0_1r1_x1i1, vec_swap_sign, x3i0_3r0_3i1_x3r1);
103 const float32x4_t x1_x3_sub =
104 vmlsq_f32(x1r0_1i0_1r1_x1i1, vec_swap_sign, x3i0_3r0_3i1_x3r1);
105 const float32x2_t yy0_a = vdup_lane_f32(vget_high_f32(x1_x3_add), 0);
106 const float32x2_t yy0_s = vdup_lane_f32(vget_high_f32(x1_x3_sub), 0);
107 const float32x4_t yy0_as = vcombine_f32(yy0_a, yy0_s);
108 const float32x2_t yy1_a = vdup_lane_f32(vget_high_f32(x1_x3_add), 1);
109 const float32x2_t yy1_s = vdup_lane_f32(vget_high_f32(x1_x3_sub), 1);
110 const float32x4_t yy1_as = vcombine_f32(yy1_a, yy1_s);
111 const float32x4_t yy0 = vmlaq_f32(yy0_as, vec_swap_sign, yy1_as);
112 const float32x4_t yy4 = vmulq_f32(wk1rv, yy0);
113 const float32x4_t xx1_rev = vrev64q_f32(xx1);
114 const float32x4_t yy4_rev = vrev64q_f32(yy4);
115
116 vst1_f32(&a[j + 0], vget_low_f32(xx0));
117 vst1_f32(&a[j + 32], vget_high_f32(xx0));
118 vst1_f32(&a[j + 16], vget_low_f32(xx1));
119 vst1_f32(&a[j + 48], vget_high_f32(xx1_rev));
120
121 a[j + 48] = -a[j + 48];
122
123 vst1_f32(&a[j + 8], vget_low_f32(x1_x3_add));
124 vst1_f32(&a[j + 24], vget_low_f32(x1_x3_sub));
125 vst1_f32(&a[j + 40], vget_low_f32(yy4));
126 vst1_f32(&a[j + 56], vget_high_f32(yy4_rev));
127 }
128
129 {
130 const int k = 64;
131 const int k1 = 2;
132 const int k2 = 2 * k1;
133 const float32x4_t wk2rv = vld1q_f32(&rdft_wk2r[k2 + 0]);
134 const float32x4_t wk2iv = vld1q_f32(&rdft_wk2i[k2 + 0]);
135 const float32x4_t wk1iv = vld1q_f32(&rdft_wk1i[k2 + 0]);
136 const float32x4_t wk3rv = vld1q_f32(&rdft_wk3r[k2 + 0]);
137 const float32x4_t wk3iv = vld1q_f32(&rdft_wk3i[k2 + 0]);
138 wk1rv = vld1q_f32(&rdft_wk1r[k2 + 0]);
139 for (j = k; j < l + k; j += 2) {
140 const float32x2_t a_00 = vld1_f32(&a[j + 0]);
141 const float32x2_t a_08 = vld1_f32(&a[j + 8]);
142 const float32x2_t a_32 = vld1_f32(&a[j + 32]);
143 const float32x2_t a_40 = vld1_f32(&a[j + 40]);
144 const float32x4_t a_00_32 = vcombine_f32(a_00, a_32);
145 const float32x4_t a_08_40 = vcombine_f32(a_08, a_40);
146 const float32x4_t x0r0_0i0_0r1_x0i1 = vaddq_f32(a_00_32, a_08_40);
147 const float32x4_t x1r0_1i0_1r1_x1i1 = vsubq_f32(a_00_32, a_08_40);
148 const float32x2_t a_16 = vld1_f32(&a[j + 16]);
149 const float32x2_t a_24 = vld1_f32(&a[j + 24]);
150 const float32x2_t a_48 = vld1_f32(&a[j + 48]);
151 const float32x2_t a_56 = vld1_f32(&a[j + 56]);
152 const float32x4_t a_16_48 = vcombine_f32(a_16, a_48);
153 const float32x4_t a_24_56 = vcombine_f32(a_24, a_56);
154 const float32x4_t x2r0_2i0_2r1_x2i1 = vaddq_f32(a_16_48, a_24_56);
155 const float32x4_t x3r0_3i0_3r1_x3i1 = vsubq_f32(a_16_48, a_24_56);
156 const float32x4_t xx = vaddq_f32(x0r0_0i0_0r1_x0i1, x2r0_2i0_2r1_x2i1);
157 const float32x4_t xx1 = vsubq_f32(x0r0_0i0_0r1_x0i1, x2r0_2i0_2r1_x2i1);
158 const float32x4_t x3i0_3r0_3i1_x3r1 = vrev64q_f32(x3r0_3i0_3r1_x3i1);
159 const float32x4_t x1_x3_add =
160 vmlaq_f32(x1r0_1i0_1r1_x1i1, vec_swap_sign, x3i0_3r0_3i1_x3r1);
161 const float32x4_t x1_x3_sub =
162 vmlsq_f32(x1r0_1i0_1r1_x1i1, vec_swap_sign, x3i0_3r0_3i1_x3r1);
163 float32x4_t xx4 = vmulq_f32(wk2rv, xx1);
164 float32x4_t xx12 = vmulq_f32(wk1rv, x1_x3_add);
165 float32x4_t xx22 = vmulq_f32(wk3rv, x1_x3_sub);
166 xx4 = vmlaq_f32(xx4, wk2iv, vrev64q_f32(xx1));
167 xx12 = vmlaq_f32(xx12, wk1iv, vrev64q_f32(x1_x3_add));
168 xx22 = vmlaq_f32(xx22, wk3iv, vrev64q_f32(x1_x3_sub));
169
170 vst1_f32(&a[j + 0], vget_low_f32(xx));
171 vst1_f32(&a[j + 32], vget_high_f32(xx));
172 vst1_f32(&a[j + 16], vget_low_f32(xx4));
173 vst1_f32(&a[j + 48], vget_high_f32(xx4));
174 vst1_f32(&a[j + 8], vget_low_f32(xx12));
175 vst1_f32(&a[j + 40], vget_high_f32(xx12));
176 vst1_f32(&a[j + 24], vget_low_f32(xx22));
177 vst1_f32(&a[j + 56], vget_high_f32(xx22));
178 }
179 }
180 }
181
reverse_order_f32x4(float32x4_t in)182 __inline static float32x4_t reverse_order_f32x4(float32x4_t in) {
183 // A B C D -> C D A B
184 const float32x4_t rev = vcombine_f32(vget_high_f32(in), vget_low_f32(in));
185 // C D A B -> D C B A
186 return vrev64q_f32(rev);
187 }
188
rftfsub_128_neon(float * a)189 void rftfsub_128_neon(float* a) {
190 const float* c = rdft_w + 32;
191 int j1, j2;
192 const float32x4_t mm_half = vdupq_n_f32(0.5f);
193
194 // Vectorized code (four at once).
195 // Note: commented number are indexes for the first iteration of the loop.
196 for (j1 = 1, j2 = 2; j2 + 7 < 64; j1 += 4, j2 += 8) {
197 // Load 'wk'.
198 const float32x4_t c_j1 = vld1q_f32(&c[j1]); // 1, 2, 3, 4,
199 const float32x4_t c_k1 = vld1q_f32(&c[29 - j1]); // 28, 29, 30, 31,
200 const float32x4_t wkrt = vsubq_f32(mm_half, c_k1); // 28, 29, 30, 31,
201 const float32x4_t wkr_ = reverse_order_f32x4(wkrt); // 31, 30, 29, 28,
202 const float32x4_t wki_ = c_j1; // 1, 2, 3, 4,
203 // Load and shuffle 'a'.
204 // 2, 4, 6, 8, 3, 5, 7, 9
205 float32x4x2_t a_j2_p = vld2q_f32(&a[0 + j2]);
206 // 120, 122, 124, 126, 121, 123, 125, 127,
207 const float32x4x2_t k2_0_4 = vld2q_f32(&a[122 - j2]);
208 // 126, 124, 122, 120
209 const float32x4_t a_k2_p0 = reverse_order_f32x4(k2_0_4.val[0]);
210 // 127, 125, 123, 121
211 const float32x4_t a_k2_p1 = reverse_order_f32x4(k2_0_4.val[1]);
212 // Calculate 'x'.
213 const float32x4_t xr_ = vsubq_f32(a_j2_p.val[0], a_k2_p0);
214 // 2-126, 4-124, 6-122, 8-120,
215 const float32x4_t xi_ = vaddq_f32(a_j2_p.val[1], a_k2_p1);
216 // 3-127, 5-125, 7-123, 9-121,
217 // Calculate product into 'y'.
218 // yr = wkr * xr - wki * xi;
219 // yi = wkr * xi + wki * xr;
220 const float32x4_t a_ = vmulq_f32(wkr_, xr_);
221 const float32x4_t b_ = vmulq_f32(wki_, xi_);
222 const float32x4_t c_ = vmulq_f32(wkr_, xi_);
223 const float32x4_t d_ = vmulq_f32(wki_, xr_);
224 const float32x4_t yr_ = vsubq_f32(a_, b_); // 2-126, 4-124, 6-122, 8-120,
225 const float32x4_t yi_ = vaddq_f32(c_, d_); // 3-127, 5-125, 7-123, 9-121,
226 // Update 'a'.
227 // a[j2 + 0] -= yr;
228 // a[j2 + 1] -= yi;
229 // a[k2 + 0] += yr;
230 // a[k2 + 1] -= yi;
231 // 126, 124, 122, 120,
232 const float32x4_t a_k2_p0n = vaddq_f32(a_k2_p0, yr_);
233 // 127, 125, 123, 121,
234 const float32x4_t a_k2_p1n = vsubq_f32(a_k2_p1, yi_);
235 // Shuffle in right order and store.
236 const float32x4_t a_k2_p0nr = vrev64q_f32(a_k2_p0n);
237 const float32x4_t a_k2_p1nr = vrev64q_f32(a_k2_p1n);
238 // 124, 125, 126, 127, 120, 121, 122, 123
239 const float32x4x2_t a_k2_n = vzipq_f32(a_k2_p0nr, a_k2_p1nr);
240 // 2, 4, 6, 8,
241 a_j2_p.val[0] = vsubq_f32(a_j2_p.val[0], yr_);
242 // 3, 5, 7, 9,
243 a_j2_p.val[1] = vsubq_f32(a_j2_p.val[1], yi_);
244 // 2, 3, 4, 5, 6, 7, 8, 9,
245 vst2q_f32(&a[0 + j2], a_j2_p);
246
247 vst1q_f32(&a[122 - j2], a_k2_n.val[1]);
248 vst1q_f32(&a[126 - j2], a_k2_n.val[0]);
249 }
250
251 // Scalar code for the remaining items.
252 for (; j2 < 64; j1 += 1, j2 += 2) {
253 const int k2 = 128 - j2;
254 const int k1 = 32 - j1;
255 const float wkr = 0.5f - c[k1];
256 const float wki = c[j1];
257 const float xr = a[j2 + 0] - a[k2 + 0];
258 const float xi = a[j2 + 1] + a[k2 + 1];
259 const float yr = wkr * xr - wki * xi;
260 const float yi = wkr * xi + wki * xr;
261 a[j2 + 0] -= yr;
262 a[j2 + 1] -= yi;
263 a[k2 + 0] += yr;
264 a[k2 + 1] -= yi;
265 }
266 }
267
rftbsub_128_neon(float * a)268 void rftbsub_128_neon(float* a) {
269 const float* c = rdft_w + 32;
270 int j1, j2;
271 const float32x4_t mm_half = vdupq_n_f32(0.5f);
272
273 a[1] = -a[1];
274 // Vectorized code (four at once).
275 // Note: commented number are indexes for the first iteration of the loop.
276 for (j1 = 1, j2 = 2; j2 + 7 < 64; j1 += 4, j2 += 8) {
277 // Load 'wk'.
278 const float32x4_t c_j1 = vld1q_f32(&c[j1]); // 1, 2, 3, 4,
279 const float32x4_t c_k1 = vld1q_f32(&c[29 - j1]); // 28, 29, 30, 31,
280 const float32x4_t wkrt = vsubq_f32(mm_half, c_k1); // 28, 29, 30, 31,
281 const float32x4_t wkr_ = reverse_order_f32x4(wkrt); // 31, 30, 29, 28,
282 const float32x4_t wki_ = c_j1; // 1, 2, 3, 4,
283 // Load and shuffle 'a'.
284 // 2, 4, 6, 8, 3, 5, 7, 9
285 float32x4x2_t a_j2_p = vld2q_f32(&a[0 + j2]);
286 // 120, 122, 124, 126, 121, 123, 125, 127,
287 const float32x4x2_t k2_0_4 = vld2q_f32(&a[122 - j2]);
288 // 126, 124, 122, 120
289 const float32x4_t a_k2_p0 = reverse_order_f32x4(k2_0_4.val[0]);
290 // 127, 125, 123, 121
291 const float32x4_t a_k2_p1 = reverse_order_f32x4(k2_0_4.val[1]);
292 // Calculate 'x'.
293 const float32x4_t xr_ = vsubq_f32(a_j2_p.val[0], a_k2_p0);
294 // 2-126, 4-124, 6-122, 8-120,
295 const float32x4_t xi_ = vaddq_f32(a_j2_p.val[1], a_k2_p1);
296 // 3-127, 5-125, 7-123, 9-121,
297 // Calculate product into 'y'.
298 // yr = wkr * xr - wki * xi;
299 // yi = wkr * xi + wki * xr;
300 const float32x4_t a_ = vmulq_f32(wkr_, xr_);
301 const float32x4_t b_ = vmulq_f32(wki_, xi_);
302 const float32x4_t c_ = vmulq_f32(wkr_, xi_);
303 const float32x4_t d_ = vmulq_f32(wki_, xr_);
304 const float32x4_t yr_ = vaddq_f32(a_, b_); // 2-126, 4-124, 6-122, 8-120,
305 const float32x4_t yi_ = vsubq_f32(c_, d_); // 3-127, 5-125, 7-123, 9-121,
306 // Update 'a'.
307 // a[j2 + 0] -= yr;
308 // a[j2 + 1] -= yi;
309 // a[k2 + 0] += yr;
310 // a[k2 + 1] -= yi;
311 // 126, 124, 122, 120,
312 const float32x4_t a_k2_p0n = vaddq_f32(a_k2_p0, yr_);
313 // 127, 125, 123, 121,
314 const float32x4_t a_k2_p1n = vsubq_f32(yi_, a_k2_p1);
315 // Shuffle in right order and store.
316 // 2, 3, 4, 5, 6, 7, 8, 9,
317 const float32x4_t a_k2_p0nr = vrev64q_f32(a_k2_p0n);
318 const float32x4_t a_k2_p1nr = vrev64q_f32(a_k2_p1n);
319 // 124, 125, 126, 127, 120, 121, 122, 123
320 const float32x4x2_t a_k2_n = vzipq_f32(a_k2_p0nr, a_k2_p1nr);
321 // 2, 4, 6, 8,
322 a_j2_p.val[0] = vsubq_f32(a_j2_p.val[0], yr_);
323 // 3, 5, 7, 9,
324 a_j2_p.val[1] = vsubq_f32(yi_, a_j2_p.val[1]);
325 // 2, 3, 4, 5, 6, 7, 8, 9,
326 vst2q_f32(&a[0 + j2], a_j2_p);
327
328 vst1q_f32(&a[122 - j2], a_k2_n.val[1]);
329 vst1q_f32(&a[126 - j2], a_k2_n.val[0]);
330 }
331
332 // Scalar code for the remaining items.
333 for (; j2 < 64; j1 += 1, j2 += 2) {
334 const int k2 = 128 - j2;
335 const int k1 = 32 - j1;
336 const float wkr = 0.5f - c[k1];
337 const float wki = c[j1];
338 const float xr = a[j2 + 0] - a[k2 + 0];
339 const float xi = a[j2 + 1] + a[k2 + 1];
340 const float yr = wkr * xr + wki * xi;
341 const float yi = wkr * xi - wki * xr;
342 a[j2 + 0] = a[j2 + 0] - yr;
343 a[j2 + 1] = yi - a[j2 + 1];
344 a[k2 + 0] = yr + a[k2 + 0];
345 a[k2 + 1] = yi - a[k2 + 1];
346 }
347 a[65] = -a[65];
348 }
349 #endif
350
351 } // namespace webrtc
352