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