1 /*
2 * Copyright (c) 2011 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 /*
13 * This file contains the function WebRtcSpl_ComplexFFT().
14 * The description header can be found in signal_processing_library.h
15 *
16 */
17
18 #include "common_audio/signal_processing/complex_fft_tables.h"
19 #include "common_audio/signal_processing/include/signal_processing_library.h"
20 #include "rtc_base/system/arch.h"
21
22 #define CFFTSFT 14
23 #define CFFTRND 1
24 #define CFFTRND2 16384
25
26 #define CIFFTSFT 14
27 #define CIFFTRND 1
28
29
WebRtcSpl_ComplexFFT(int16_t frfi[],int stages,int mode)30 int WebRtcSpl_ComplexFFT(int16_t frfi[], int stages, int mode)
31 {
32 int i, j, l, k, istep, n, m;
33 int16_t wr, wi;
34 int32_t tr32, ti32, qr32, qi32;
35
36 /* The 1024-value is a constant given from the size of kSinTable1024[],
37 * and should not be changed depending on the input parameter 'stages'
38 */
39 n = 1 << stages;
40 if (n > 1024)
41 return -1;
42
43 l = 1;
44 k = 10 - 1; /* Constant for given kSinTable1024[]. Do not change
45 depending on the input parameter 'stages' */
46
47 if (mode == 0)
48 {
49 // mode==0: Low-complexity and Low-accuracy mode
50 while (l < n)
51 {
52 istep = l << 1;
53
54 for (m = 0; m < l; ++m)
55 {
56 j = m << k;
57
58 /* The 256-value is a constant given as 1/4 of the size of
59 * kSinTable1024[], and should not be changed depending on the input
60 * parameter 'stages'. It will result in 0 <= j < N_SINE_WAVE/2
61 */
62 wr = kSinTable1024[j + 256];
63 wi = -kSinTable1024[j];
64
65 for (i = m; i < n; i += istep)
66 {
67 j = i + l;
68
69 tr32 = (wr * frfi[2 * j] - wi * frfi[2 * j + 1]) >> 15;
70
71 ti32 = (wr * frfi[2 * j + 1] + wi * frfi[2 * j]) >> 15;
72
73 qr32 = (int32_t)frfi[2 * i];
74 qi32 = (int32_t)frfi[2 * i + 1];
75 frfi[2 * j] = (int16_t)((qr32 - tr32) >> 1);
76 frfi[2 * j + 1] = (int16_t)((qi32 - ti32) >> 1);
77 frfi[2 * i] = (int16_t)((qr32 + tr32) >> 1);
78 frfi[2 * i + 1] = (int16_t)((qi32 + ti32) >> 1);
79 }
80 }
81
82 --k;
83 l = istep;
84
85 }
86
87 } else
88 {
89 // mode==1: High-complexity and High-accuracy mode
90 while (l < n)
91 {
92 istep = l << 1;
93
94 for (m = 0; m < l; ++m)
95 {
96 j = m << k;
97
98 /* The 256-value is a constant given as 1/4 of the size of
99 * kSinTable1024[], and should not be changed depending on the input
100 * parameter 'stages'. It will result in 0 <= j < N_SINE_WAVE/2
101 */
102 wr = kSinTable1024[j + 256];
103 wi = -kSinTable1024[j];
104
105 #ifdef WEBRTC_ARCH_ARM_V7
106 int32_t wri = 0;
107 __asm __volatile("pkhbt %0, %1, %2, lsl #16" : "=r"(wri) :
108 "r"((int32_t)wr), "r"((int32_t)wi));
109 #endif
110
111 for (i = m; i < n; i += istep)
112 {
113 j = i + l;
114
115 #ifdef WEBRTC_ARCH_ARM_V7
116 register int32_t frfi_r;
117 __asm __volatile(
118 "pkhbt %[frfi_r], %[frfi_even], %[frfi_odd],"
119 " lsl #16\n\t"
120 "smlsd %[tr32], %[wri], %[frfi_r], %[cfftrnd]\n\t"
121 "smladx %[ti32], %[wri], %[frfi_r], %[cfftrnd]\n\t"
122 :[frfi_r]"=&r"(frfi_r),
123 [tr32]"=&r"(tr32),
124 [ti32]"=r"(ti32)
125 :[frfi_even]"r"((int32_t)frfi[2*j]),
126 [frfi_odd]"r"((int32_t)frfi[2*j +1]),
127 [wri]"r"(wri),
128 [cfftrnd]"r"(CFFTRND));
129 #else
130 tr32 = wr * frfi[2 * j] - wi * frfi[2 * j + 1] + CFFTRND;
131
132 ti32 = wr * frfi[2 * j + 1] + wi * frfi[2 * j] + CFFTRND;
133 #endif
134
135 tr32 >>= 15 - CFFTSFT;
136 ti32 >>= 15 - CFFTSFT;
137
138 qr32 = ((int32_t)frfi[2 * i]) * (1 << CFFTSFT);
139 qi32 = ((int32_t)frfi[2 * i + 1]) * (1 << CFFTSFT);
140
141 frfi[2 * j] = (int16_t)(
142 (qr32 - tr32 + CFFTRND2) >> (1 + CFFTSFT));
143 frfi[2 * j + 1] = (int16_t)(
144 (qi32 - ti32 + CFFTRND2) >> (1 + CFFTSFT));
145 frfi[2 * i] = (int16_t)(
146 (qr32 + tr32 + CFFTRND2) >> (1 + CFFTSFT));
147 frfi[2 * i + 1] = (int16_t)(
148 (qi32 + ti32 + CFFTRND2) >> (1 + CFFTSFT));
149 }
150 }
151
152 --k;
153 l = istep;
154 }
155 }
156 return 0;
157 }
158
WebRtcSpl_ComplexIFFT(int16_t frfi[],int stages,int mode)159 int WebRtcSpl_ComplexIFFT(int16_t frfi[], int stages, int mode)
160 {
161 size_t i, j, l, istep, n, m;
162 int k, scale, shift;
163 int16_t wr, wi;
164 int32_t tr32, ti32, qr32, qi32;
165 int32_t tmp32, round2;
166
167 /* The 1024-value is a constant given from the size of kSinTable1024[],
168 * and should not be changed depending on the input parameter 'stages'
169 */
170 n = ((size_t)1) << stages;
171 if (n > 1024)
172 return -1;
173
174 scale = 0;
175
176 l = 1;
177 k = 10 - 1; /* Constant for given kSinTable1024[]. Do not change
178 depending on the input parameter 'stages' */
179
180 while (l < n)
181 {
182 // variable scaling, depending upon data
183 shift = 0;
184 round2 = 8192;
185
186 tmp32 = WebRtcSpl_MaxAbsValueW16(frfi, 2 * n);
187 if (tmp32 > 13573)
188 {
189 shift++;
190 scale++;
191 round2 <<= 1;
192 }
193 if (tmp32 > 27146)
194 {
195 shift++;
196 scale++;
197 round2 <<= 1;
198 }
199
200 istep = l << 1;
201
202 if (mode == 0)
203 {
204 // mode==0: Low-complexity and Low-accuracy mode
205 for (m = 0; m < l; ++m)
206 {
207 j = m << k;
208
209 /* The 256-value is a constant given as 1/4 of the size of
210 * kSinTable1024[], and should not be changed depending on the input
211 * parameter 'stages'. It will result in 0 <= j < N_SINE_WAVE/2
212 */
213 wr = kSinTable1024[j + 256];
214 wi = kSinTable1024[j];
215
216 for (i = m; i < n; i += istep)
217 {
218 j = i + l;
219
220 tr32 = (wr * frfi[2 * j] - wi * frfi[2 * j + 1]) >> 15;
221
222 ti32 = (wr * frfi[2 * j + 1] + wi * frfi[2 * j]) >> 15;
223
224 qr32 = (int32_t)frfi[2 * i];
225 qi32 = (int32_t)frfi[2 * i + 1];
226 frfi[2 * j] = (int16_t)((qr32 - tr32) >> shift);
227 frfi[2 * j + 1] = (int16_t)((qi32 - ti32) >> shift);
228 frfi[2 * i] = (int16_t)((qr32 + tr32) >> shift);
229 frfi[2 * i + 1] = (int16_t)((qi32 + ti32) >> shift);
230 }
231 }
232 } else
233 {
234 // mode==1: High-complexity and High-accuracy mode
235
236 for (m = 0; m < l; ++m)
237 {
238 j = m << k;
239
240 /* The 256-value is a constant given as 1/4 of the size of
241 * kSinTable1024[], and should not be changed depending on the input
242 * parameter 'stages'. It will result in 0 <= j < N_SINE_WAVE/2
243 */
244 wr = kSinTable1024[j + 256];
245 wi = kSinTable1024[j];
246
247 #ifdef WEBRTC_ARCH_ARM_V7
248 int32_t wri = 0;
249 __asm __volatile("pkhbt %0, %1, %2, lsl #16" : "=r"(wri) :
250 "r"((int32_t)wr), "r"((int32_t)wi));
251 #endif
252
253 for (i = m; i < n; i += istep)
254 {
255 j = i + l;
256
257 #ifdef WEBRTC_ARCH_ARM_V7
258 register int32_t frfi_r;
259 __asm __volatile(
260 "pkhbt %[frfi_r], %[frfi_even], %[frfi_odd], lsl #16\n\t"
261 "smlsd %[tr32], %[wri], %[frfi_r], %[cifftrnd]\n\t"
262 "smladx %[ti32], %[wri], %[frfi_r], %[cifftrnd]\n\t"
263 :[frfi_r]"=&r"(frfi_r),
264 [tr32]"=&r"(tr32),
265 [ti32]"=r"(ti32)
266 :[frfi_even]"r"((int32_t)frfi[2*j]),
267 [frfi_odd]"r"((int32_t)frfi[2*j +1]),
268 [wri]"r"(wri),
269 [cifftrnd]"r"(CIFFTRND)
270 );
271 #else
272
273 tr32 = wr * frfi[2 * j] - wi * frfi[2 * j + 1] + CIFFTRND;
274
275 ti32 = wr * frfi[2 * j + 1] + wi * frfi[2 * j] + CIFFTRND;
276 #endif
277 tr32 >>= 15 - CIFFTSFT;
278 ti32 >>= 15 - CIFFTSFT;
279
280 qr32 = ((int32_t)frfi[2 * i]) * (1 << CIFFTSFT);
281 qi32 = ((int32_t)frfi[2 * i + 1]) * (1 << CIFFTSFT);
282
283 frfi[2 * j] = (int16_t)(
284 (qr32 - tr32 + round2) >> (shift + CIFFTSFT));
285 frfi[2 * j + 1] = (int16_t)(
286 (qi32 - ti32 + round2) >> (shift + CIFFTSFT));
287 frfi[2 * i] = (int16_t)(
288 (qr32 + tr32 + round2) >> (shift + CIFFTSFT));
289 frfi[2 * i + 1] = (int16_t)(
290 (qi32 + ti32 + round2) >> (shift + CIFFTSFT));
291 }
292 }
293
294 }
295 --k;
296 l = istep;
297 }
298 return scale;
299 }
300