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