/src/opus/celt/kiss_fft.c
Line | Count | Source |
1 | | /*Copyright (c) 2003-2004, Mark Borgerding |
2 | | Lots of modifications by Jean-Marc Valin |
3 | | Copyright (c) 2005-2007, Xiph.Org Foundation |
4 | | Copyright (c) 2008, Xiph.Org Foundation, CSIRO |
5 | | |
6 | | All rights reserved. |
7 | | |
8 | | Redistribution and use in source and binary forms, with or without |
9 | | modification, are permitted provided that the following conditions are met: |
10 | | |
11 | | * Redistributions of source code must retain the above copyright notice, |
12 | | this list of conditions and the following disclaimer. |
13 | | * Redistributions in binary form must reproduce the above copyright notice, |
14 | | this list of conditions and the following disclaimer in the |
15 | | documentation and/or other materials provided with the distribution. |
16 | | |
17 | | THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" |
18 | | AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
19 | | IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
20 | | ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE |
21 | | LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR |
22 | | CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF |
23 | | SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS |
24 | | INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN |
25 | | CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) |
26 | | ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE |
27 | | POSSIBILITY OF SUCH DAMAGE.*/ |
28 | | |
29 | | /* This code is originally from Mark Borgerding's KISS-FFT but has been |
30 | | heavily modified to better suit Opus */ |
31 | | |
32 | | #ifndef SKIP_CONFIG_H |
33 | | # ifdef HAVE_CONFIG_H |
34 | | # include "config.h" |
35 | | # endif |
36 | | #endif |
37 | | |
38 | | #include "_kiss_fft_guts.h" |
39 | | #include "arch.h" |
40 | | #include "os_support.h" |
41 | | #include "mathops.h" |
42 | | #include "stack_alloc.h" |
43 | | |
44 | | #if defined(CUSTOM_MODES) || defined(ENABLE_OPUS_CUSTOM_API) |
45 | | #define OPUS_CUSTOM 1 |
46 | | #endif |
47 | | |
48 | | |
49 | | #if defined(ENABLE_PFA) |
50 | | extern void opus_fft_pfa_c(const kiss_fft_state *st, const kiss_fft_cpx *fin, kiss_fft_cpx *fout ARG_FIXED(int downshift)); |
51 | | extern void opus_ifft_pfa_c(const kiss_fft_state *st, const kiss_fft_cpx *fin, kiss_fft_cpx *fout ARG_FIXED(int fft_shift)); |
52 | | #endif |
53 | | |
54 | | #ifndef M_PI |
55 | | #define M_PI 3.141592653 |
56 | | #endif |
57 | | |
58 | | /* The guts header contains all the multiplication and addition macros that are defined for |
59 | | complex numbers. It also declares the kf_ internal functions. |
60 | | */ |
61 | | |
62 | | #if !defined(ENABLE_PFA) || defined(OPUS_CUSTOM) || defined(ENABLE_DEEP_PLC) |
63 | | static void kf_bfly2( |
64 | | kiss_fft_cpx * Fout, |
65 | | int m, |
66 | | int N |
67 | | ) |
68 | 77.9k | { |
69 | 77.9k | kiss_fft_cpx * Fout2; |
70 | 77.9k | int i; |
71 | 77.9k | (void)m; |
72 | | #ifdef CUSTOM_MODES |
73 | | if (m==1) |
74 | | { |
75 | | celt_assert(m==1); |
76 | | for (i=0;i<N;i++) |
77 | | { |
78 | | kiss_fft_cpx t; |
79 | | Fout2 = Fout + 1; |
80 | | t = *Fout2; |
81 | | C_SUB( *Fout2 , *Fout , t ); |
82 | | C_ADDTO( *Fout , t ); |
83 | | Fout += 2; |
84 | | } |
85 | | } else |
86 | | #endif |
87 | 77.9k | { |
88 | 77.9k | celt_coef tw; |
89 | 77.9k | tw = QCONST32(0.7071067812f, COEF_SHIFT-1); |
90 | | /* We know that m==4 here because the radix-2 is just after a radix-4 */ |
91 | 77.9k | celt_assert(m==4); |
92 | 3.40M | for (i=0;i<N;i++) |
93 | 3.32M | { |
94 | 3.32M | kiss_fft_cpx t; |
95 | 3.32M | Fout2 = Fout + 4; |
96 | 3.32M | t = Fout2[0]; |
97 | 3.32M | C_SUB( Fout2[0] , Fout[0] , t ); |
98 | 3.32M | C_ADDTO( Fout[0] , t ); |
99 | | |
100 | 3.32M | t.r = S_MUL(ADD32_ovflw(Fout2[1].r, Fout2[1].i), tw); |
101 | 3.32M | t.i = S_MUL(SUB32_ovflw(Fout2[1].i, Fout2[1].r), tw); |
102 | 3.32M | C_SUB( Fout2[1] , Fout[1] , t ); |
103 | 3.32M | C_ADDTO( Fout[1] , t ); |
104 | | |
105 | 3.32M | t.r = Fout2[2].i; |
106 | 3.32M | t.i = NEG32_ovflw(Fout2[2].r); |
107 | 3.32M | C_SUB( Fout2[2] , Fout[2] , t ); |
108 | 3.32M | C_ADDTO( Fout[2] , t ); |
109 | | |
110 | 3.32M | t.r = S_MUL(SUB32_ovflw(Fout2[3].i, Fout2[3].r), tw); |
111 | 3.32M | t.i = S_MUL(NEG32_ovflw(ADD32_ovflw(Fout2[3].i, Fout2[3].r)), tw); |
112 | 3.32M | C_SUB( Fout2[3] , Fout[3] , t ); |
113 | 3.32M | C_ADDTO( Fout[3] , t ); |
114 | 3.32M | Fout += 8; |
115 | 3.32M | } |
116 | 77.9k | } |
117 | 77.9k | } |
118 | | |
119 | | static void kf_bfly4( |
120 | | kiss_fft_cpx * Fout, |
121 | | const size_t fstride, |
122 | | const kiss_fft_state *st, |
123 | | int m, |
124 | | int N, |
125 | | int mm |
126 | | ) |
127 | 542k | { |
128 | 542k | int i; |
129 | | |
130 | 542k | if (m==1) |
131 | 325k | { |
132 | | /* Degenerate case where all the twiddles are 1. */ |
133 | 18.2M | for (i=0;i<N;i++) |
134 | 17.9M | { |
135 | 17.9M | kiss_fft_cpx scratch0, scratch1; |
136 | | |
137 | 17.9M | C_SUB( scratch0 , *Fout, Fout[2] ); |
138 | 17.9M | C_ADDTO(*Fout, Fout[2]); |
139 | 17.9M | C_ADD( scratch1 , Fout[1] , Fout[3] ); |
140 | 17.9M | C_SUB( Fout[2], *Fout, scratch1 ); |
141 | 17.9M | C_ADDTO( *Fout , scratch1 ); |
142 | 17.9M | C_SUB( scratch1 , Fout[1] , Fout[3] ); |
143 | | |
144 | 17.9M | Fout[1].r = ADD32_ovflw(scratch0.r, scratch1.i); |
145 | 17.9M | Fout[1].i = SUB32_ovflw(scratch0.i, scratch1.r); |
146 | 17.9M | Fout[3].r = SUB32_ovflw(scratch0.r, scratch1.i); |
147 | 17.9M | Fout[3].i = ADD32_ovflw(scratch0.i, scratch1.r); |
148 | 17.9M | Fout+=4; |
149 | 17.9M | } |
150 | 325k | } else { |
151 | 217k | int j; |
152 | 217k | kiss_fft_cpx scratch[6]; |
153 | 217k | const kiss_twiddle_cpx *tw1,*tw2,*tw3; |
154 | 217k | const int m2=2*m; |
155 | 217k | const int m3=3*m; |
156 | 217k | kiss_fft_cpx * Fout_beg = Fout; |
157 | 3.47M | for (i=0;i<N;i++) |
158 | 3.25M | { |
159 | 3.25M | Fout = Fout_beg + i*mm; |
160 | 3.25M | tw3 = tw2 = tw1 = st->twiddles; |
161 | | /* m is guaranteed to be a multiple of 4. */ |
162 | 19.1M | for (j=0;j<m;j++) |
163 | 15.9M | { |
164 | 15.9M | C_MUL(scratch[0],Fout[m] , *tw1 ); |
165 | 15.9M | C_MUL(scratch[1],Fout[m2] , *tw2 ); |
166 | 15.9M | C_MUL(scratch[2],Fout[m3] , *tw3 ); |
167 | | |
168 | 15.9M | C_SUB( scratch[5] , *Fout, scratch[1] ); |
169 | 15.9M | C_ADDTO(*Fout, scratch[1]); |
170 | 15.9M | C_ADD( scratch[3] , scratch[0] , scratch[2] ); |
171 | 15.9M | C_SUB( scratch[4] , scratch[0] , scratch[2] ); |
172 | 15.9M | C_SUB( Fout[m2], *Fout, scratch[3] ); |
173 | 15.9M | tw1 += fstride; |
174 | 15.9M | tw2 += fstride*2; |
175 | 15.9M | tw3 += fstride*3; |
176 | 15.9M | C_ADDTO( *Fout , scratch[3] ); |
177 | | |
178 | 15.9M | Fout[m].r = ADD32_ovflw(scratch[5].r, scratch[4].i); |
179 | 15.9M | Fout[m].i = SUB32_ovflw(scratch[5].i, scratch[4].r); |
180 | 15.9M | Fout[m3].r = SUB32_ovflw(scratch[5].r, scratch[4].i); |
181 | 15.9M | Fout[m3].i = ADD32_ovflw(scratch[5].i, scratch[4].r); |
182 | 15.9M | ++Fout; |
183 | 15.9M | } |
184 | 3.25M | } |
185 | 217k | } |
186 | 542k | } |
187 | | |
188 | | |
189 | | #ifndef RADIX_TWO_ONLY |
190 | | |
191 | | static void kf_bfly3( |
192 | | kiss_fft_cpx * Fout, |
193 | | const size_t fstride, |
194 | | const kiss_fft_state *st, |
195 | | int m, |
196 | | int N, |
197 | | int mm |
198 | | ) |
199 | 325k | { |
200 | 325k | int i; |
201 | 325k | size_t k; |
202 | 325k | const size_t m2 = 2*m; |
203 | 325k | const kiss_twiddle_cpx *tw1,*tw2; |
204 | 325k | kiss_fft_cpx scratch[5]; |
205 | 325k | kiss_twiddle_cpx epi3; |
206 | | |
207 | 325k | kiss_fft_cpx * Fout_beg = Fout; |
208 | | #ifdef FIXED_POINT |
209 | | /*epi3.r = -16384;*/ /* Unused */ |
210 | | epi3.i = -QCONST32(0.86602540f, COEF_SHIFT-1); |
211 | | #else |
212 | 325k | epi3 = st->twiddles[fstride*m]; |
213 | 325k | #endif |
214 | 1.95M | for (i=0;i<N;i++) |
215 | 1.62M | { |
216 | 1.62M | Fout = Fout_beg + i*mm; |
217 | 1.62M | tw1=tw2=st->twiddles; |
218 | | /* For non-custom modes, m is guaranteed to be a multiple of 4. */ |
219 | 1.62M | k=m; |
220 | 23.9M | do { |
221 | | |
222 | 23.9M | C_MUL(scratch[1],Fout[m] , *tw1); |
223 | 23.9M | C_MUL(scratch[2],Fout[m2] , *tw2); |
224 | | |
225 | 23.9M | C_ADD(scratch[3],scratch[1],scratch[2]); |
226 | 23.9M | C_SUB(scratch[0],scratch[1],scratch[2]); |
227 | 23.9M | tw1 += fstride; |
228 | 23.9M | tw2 += fstride*2; |
229 | | |
230 | 23.9M | Fout[m].r = SUB32_ovflw(Fout->r, HALF_OF(scratch[3].r)); |
231 | 23.9M | Fout[m].i = SUB32_ovflw(Fout->i, HALF_OF(scratch[3].i)); |
232 | | |
233 | 23.9M | C_MULBYSCALAR( scratch[0] , epi3.i ); |
234 | | |
235 | 23.9M | C_ADDTO(*Fout,scratch[3]); |
236 | | |
237 | 23.9M | Fout[m2].r = ADD32_ovflw(Fout[m].r, scratch[0].i); |
238 | 23.9M | Fout[m2].i = SUB32_ovflw(Fout[m].i, scratch[0].r); |
239 | | |
240 | 23.9M | Fout[m].r = SUB32_ovflw(Fout[m].r, scratch[0].i); |
241 | 23.9M | Fout[m].i = ADD32_ovflw(Fout[m].i, scratch[0].r); |
242 | | |
243 | 23.9M | ++Fout; |
244 | 23.9M | } while(--k); |
245 | 1.62M | } |
246 | 325k | } |
247 | | |
248 | | |
249 | | #ifndef OVERRIDE_kf_bfly5 |
250 | | static void kf_bfly5( |
251 | | kiss_fft_cpx * Fout, |
252 | | const size_t fstride, |
253 | | const kiss_fft_state *st, |
254 | | int m, |
255 | | int N, |
256 | | int mm |
257 | | ) |
258 | 325k | { |
259 | 325k | kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4; |
260 | 325k | int i, u; |
261 | 325k | kiss_fft_cpx scratch[13]; |
262 | 325k | const kiss_twiddle_cpx *tw; |
263 | 325k | kiss_twiddle_cpx ya,yb; |
264 | 325k | kiss_fft_cpx * Fout_beg = Fout; |
265 | | |
266 | | #ifdef FIXED_POINT |
267 | | ya.r = QCONST32(0.30901699f, COEF_SHIFT-1); |
268 | | ya.i = -QCONST32(0.95105652f, COEF_SHIFT-1); |
269 | | yb.r = -QCONST32(0.80901699f, COEF_SHIFT-1); |
270 | | yb.i = -QCONST32(0.58778525f, COEF_SHIFT-1); |
271 | | #else |
272 | 325k | ya = st->twiddles[fstride*m]; |
273 | 325k | yb = st->twiddles[fstride*2*m]; |
274 | 325k | #endif |
275 | 325k | tw=st->twiddles; |
276 | | |
277 | 650k | for (i=0;i<N;i++) |
278 | 325k | { |
279 | 325k | Fout = Fout_beg + i*mm; |
280 | 325k | Fout0=Fout; |
281 | 325k | Fout1=Fout0+m; |
282 | 325k | Fout2=Fout0+2*m; |
283 | 325k | Fout3=Fout0+3*m; |
284 | 325k | Fout4=Fout0+4*m; |
285 | | |
286 | | /* For non-custom modes, m is guaranteed to be a multiple of 4. */ |
287 | 14.7M | for ( u=0; u<m; ++u ) { |
288 | 14.3M | scratch[0] = *Fout0; |
289 | | |
290 | 14.3M | C_MUL(scratch[1] ,*Fout1, tw[u*fstride]); |
291 | 14.3M | C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]); |
292 | 14.3M | C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]); |
293 | 14.3M | C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]); |
294 | | |
295 | 14.3M | C_ADD( scratch[7],scratch[1],scratch[4]); |
296 | 14.3M | C_SUB( scratch[10],scratch[1],scratch[4]); |
297 | 14.3M | C_ADD( scratch[8],scratch[2],scratch[3]); |
298 | 14.3M | C_SUB( scratch[9],scratch[2],scratch[3]); |
299 | | |
300 | 14.3M | Fout0->r = ADD32_ovflw(Fout0->r, ADD32_ovflw(scratch[7].r, scratch[8].r)); |
301 | 14.3M | Fout0->i = ADD32_ovflw(Fout0->i, ADD32_ovflw(scratch[7].i, scratch[8].i)); |
302 | | |
303 | 14.3M | scratch[5].r = ADD32_ovflw(scratch[0].r, ADD32_ovflw(S_MUL(scratch[7].r,ya.r), S_MUL(scratch[8].r,yb.r))); |
304 | 14.3M | scratch[5].i = ADD32_ovflw(scratch[0].i, ADD32_ovflw(S_MUL(scratch[7].i,ya.r), S_MUL(scratch[8].i,yb.r))); |
305 | | |
306 | 14.3M | scratch[6].r = ADD32_ovflw(S_MUL(scratch[10].i,ya.i), S_MUL(scratch[9].i,yb.i)); |
307 | 14.3M | scratch[6].i = NEG32_ovflw(ADD32_ovflw(S_MUL(scratch[10].r,ya.i), S_MUL(scratch[9].r,yb.i))); |
308 | | |
309 | 14.3M | C_SUB(*Fout1,scratch[5],scratch[6]); |
310 | 14.3M | C_ADD(*Fout4,scratch[5],scratch[6]); |
311 | | |
312 | 14.3M | scratch[11].r = ADD32_ovflw(scratch[0].r, ADD32_ovflw(S_MUL(scratch[7].r,yb.r), S_MUL(scratch[8].r,ya.r))); |
313 | 14.3M | scratch[11].i = ADD32_ovflw(scratch[0].i, ADD32_ovflw(S_MUL(scratch[7].i,yb.r), S_MUL(scratch[8].i,ya.r))); |
314 | 14.3M | scratch[12].r = SUB32_ovflw(S_MUL(scratch[9].i,ya.i), S_MUL(scratch[10].i,yb.i)); |
315 | 14.3M | scratch[12].i = SUB32_ovflw(S_MUL(scratch[10].r,yb.i), S_MUL(scratch[9].r,ya.i)); |
316 | | |
317 | 14.3M | C_ADD(*Fout2,scratch[11],scratch[12]); |
318 | 14.3M | C_SUB(*Fout3,scratch[11],scratch[12]); |
319 | | |
320 | 14.3M | ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4; |
321 | 14.3M | } |
322 | 325k | } |
323 | 325k | } |
324 | | #endif /* OVERRIDE_kf_bfly5 */ |
325 | | |
326 | | |
327 | | #endif |
328 | | |
329 | | |
330 | | #ifdef CUSTOM_MODES |
331 | | |
332 | | static |
333 | | void compute_bitrev_table( |
334 | | int Fout, |
335 | | opus_int16 *f, |
336 | | const size_t fstride, |
337 | | int in_stride, |
338 | | opus_int16 * factors, |
339 | | const kiss_fft_state *st |
340 | | ) |
341 | | { |
342 | | const int p=*factors++; /* the radix */ |
343 | | const int m=*factors++; /* stage's fft length/p */ |
344 | | |
345 | | /*printf ("fft %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N);*/ |
346 | | if (m==1) |
347 | | { |
348 | | int j; |
349 | | for (j=0;j<p;j++) |
350 | | { |
351 | | *f = Fout+j; |
352 | | f += fstride*in_stride; |
353 | | } |
354 | | } else { |
355 | | int j; |
356 | | for (j=0;j<p;j++) |
357 | | { |
358 | | compute_bitrev_table( Fout , f, fstride*p, in_stride, factors,st); |
359 | | f += fstride*in_stride; |
360 | | Fout += m; |
361 | | } |
362 | | } |
363 | | } |
364 | | |
365 | | /* facbuf is populated by p1,m1,p2,m2, ... |
366 | | where |
367 | | p[i] * m[i] = m[i-1] |
368 | | m0 = n */ |
369 | | static |
370 | | int kf_factor(int n,opus_int16 * facbuf) |
371 | | { |
372 | | int p=4; |
373 | | int i; |
374 | | int stages=0; |
375 | | int nbak = n; |
376 | | |
377 | | /*factor out powers of 4, powers of 2, then any remaining primes */ |
378 | | do { |
379 | | while (n % p) { |
380 | | switch (p) { |
381 | | case 4: p = 2; break; |
382 | | case 2: p = 3; break; |
383 | | default: p += 2; break; |
384 | | } |
385 | | if (p>32000 || (opus_int32)p*(opus_int32)p > n) |
386 | | p = n; /* no more factors, skip to end */ |
387 | | } |
388 | | n /= p; |
389 | | #ifdef RADIX_TWO_ONLY |
390 | | if (p!=2 && p != 4) |
391 | | #else |
392 | | if (p>5) |
393 | | #endif |
394 | | { |
395 | | return 0; |
396 | | } |
397 | | facbuf[2*stages] = p; |
398 | | if (p==2 && stages > 1) |
399 | | { |
400 | | facbuf[2*stages] = 4; |
401 | | facbuf[2] = 2; |
402 | | } |
403 | | stages++; |
404 | | } while (n > 1); |
405 | | n = nbak; |
406 | | /* Reverse the order to get the radix 4 at the end, so we can use the |
407 | | fast degenerate case. It turns out that reversing the order also |
408 | | improves the noise behaviour. */ |
409 | | for (i=0;i<stages/2;i++) |
410 | | { |
411 | | int tmp; |
412 | | tmp = facbuf[2*i]; |
413 | | facbuf[2*i] = facbuf[2*(stages-i-1)]; |
414 | | facbuf[2*(stages-i-1)] = tmp; |
415 | | } |
416 | | for (i=0;i<stages;i++) |
417 | | { |
418 | | n /= facbuf[2*i]; |
419 | | facbuf[2*i+1] = n; |
420 | | } |
421 | | return 1; |
422 | | } |
423 | | |
424 | | static void compute_twiddles(kiss_twiddle_cpx *twiddles, int nfft) |
425 | | { |
426 | | int i; |
427 | | #ifdef FIXED_POINT |
428 | | for (i=0;i<nfft;++i) { |
429 | | opus_val32 phase = -i; |
430 | | #ifdef ENABLE_QEXT |
431 | | twiddles[i].r = (int)MIN32(2147483647, floor(.5+2147483648*cos((2*M_PI/nfft)*phase))); |
432 | | twiddles[i].i = (int)MIN32(2147483647, floor(.5+2147483648*sin((2*M_PI/nfft)*phase))); |
433 | | #else |
434 | | kf_cexp2(twiddles+i, DIV32(SHL32(phase,17),nfft)); |
435 | | #endif |
436 | | } |
437 | | #else |
438 | | for (i=0;i<nfft;++i) { |
439 | | const double pi=3.14159265358979323846264338327; |
440 | | double phase = ( -2*pi /nfft ) * i; |
441 | | kf_cexp(twiddles+i, phase ); |
442 | | } |
443 | | #endif |
444 | | } |
445 | | |
446 | | int opus_fft_alloc_arch_c(kiss_fft_state *st) { |
447 | | (void)st; |
448 | | return 0; |
449 | | } |
450 | | |
451 | | /* |
452 | | * |
453 | | * Allocates all necessary storage space for the fft and ifft. |
454 | | * The return value is a contiguous block of memory. As such, |
455 | | * It can be freed with free(). |
456 | | * */ |
457 | | kiss_fft_state *opus_fft_alloc_twiddles(int nfft,void * mem,size_t * lenmem, |
458 | | const kiss_fft_state *base, int arch) |
459 | | { |
460 | | kiss_fft_state *st=NULL; |
461 | | size_t memneeded = sizeof(struct kiss_fft_state); /* twiddle factors*/ |
462 | | |
463 | | if ( lenmem==NULL ) { |
464 | | st = ( kiss_fft_state*)KISS_FFT_MALLOC( memneeded ); |
465 | | }else{ |
466 | | if (mem != NULL && *lenmem >= memneeded) |
467 | | st = (kiss_fft_state*)mem; |
468 | | *lenmem = memneeded; |
469 | | } |
470 | | if (st) { |
471 | | opus_int16 *bitrev; |
472 | | kiss_twiddle_cpx *twiddles; |
473 | | |
474 | | st->nfft=nfft; |
475 | | #ifdef FIXED_POINT |
476 | | st->scale_shift = celt_ilog2(st->nfft); |
477 | | # ifdef ENABLE_QEXT |
478 | | if (st->nfft == 1<<st->scale_shift) |
479 | | st->scale = QCONST32(1.0f, 30); |
480 | | else |
481 | | st->scale = (((opus_int64)1073741824<<st->scale_shift)+st->nfft/2)/st->nfft; |
482 | | # else |
483 | | if (st->nfft == 1<<st->scale_shift) |
484 | | st->scale = Q15ONE; |
485 | | else |
486 | | st->scale = (1073741824+st->nfft/2)/st->nfft>>(15-st->scale_shift); |
487 | | # endif |
488 | | #else |
489 | | st->scale = 1.f/nfft; |
490 | | #endif |
491 | | #if defined(ENABLE_PFA) |
492 | | if (nfft == 60 || nfft == 120 || nfft == 240 || nfft == 480 |
493 | | #if defined(ENABLE_QEXT) |
494 | | || nfft == 960 |
495 | | #endif |
496 | | ) |
497 | | { |
498 | | st->twiddles = NULL; |
499 | | st->shift = -1; |
500 | | st->bitrev = NULL; |
501 | | } else |
502 | | #endif |
503 | | { |
504 | | if (base != NULL && base->twiddles != NULL) |
505 | | { |
506 | | st->twiddles = base->twiddles; |
507 | | st->shift = 0; |
508 | | while (st->shift < 32 && nfft<<st->shift != base->nfft) |
509 | | st->shift++; |
510 | | if (st->shift>=32) |
511 | | goto fail; |
512 | | } else { |
513 | | st->twiddles = twiddles = (kiss_twiddle_cpx*)KISS_FFT_MALLOC(sizeof(kiss_twiddle_cpx)*nfft); |
514 | | compute_twiddles(twiddles, nfft); |
515 | | st->shift = -1; |
516 | | } |
517 | | if (!kf_factor(nfft,st->factors)) |
518 | | { |
519 | | goto fail; |
520 | | } |
521 | | |
522 | | /* bitrev */ |
523 | | st->bitrev = bitrev = (opus_int16*)KISS_FFT_MALLOC(sizeof(opus_int16)*nfft); |
524 | | if (st->bitrev==NULL) |
525 | | goto fail; |
526 | | compute_bitrev_table(0, bitrev, 1,1, st->factors,st); |
527 | | } |
528 | | |
529 | | /* Initialize architecture specific fft parameters */ |
530 | | if (opus_fft_alloc_arch(st, arch)) |
531 | | goto fail; |
532 | | } |
533 | | return st; |
534 | | fail: |
535 | | opus_fft_free(st, arch); |
536 | | return NULL; |
537 | | } |
538 | | |
539 | | kiss_fft_state *opus_fft_alloc(int nfft,void * mem,size_t * lenmem, int arch) |
540 | | { |
541 | | return opus_fft_alloc_twiddles(nfft, mem, lenmem, NULL, arch); |
542 | | } |
543 | | |
544 | | void opus_fft_free_arch_c(kiss_fft_state *st) { |
545 | | (void)st; |
546 | | } |
547 | | |
548 | | void opus_fft_free(const kiss_fft_state *cfg, int arch) |
549 | | { |
550 | | if (cfg) |
551 | | { |
552 | | opus_fft_free_arch((kiss_fft_state *)cfg, arch); |
553 | | opus_free((opus_int16*)cfg->bitrev); |
554 | | if (cfg->shift < 0) |
555 | | opus_free((kiss_twiddle_cpx*)cfg->twiddles); |
556 | | opus_free((kiss_fft_state*)cfg); |
557 | | } |
558 | | } |
559 | | |
560 | | #endif /* CUSTOM_MODES */ |
561 | | |
562 | | #if !defined(ENABLE_PFA) || defined(OPUS_CUSTOM) || defined(ENABLE_DEEP_PLC) |
563 | | #ifdef FIXED_POINT |
564 | | #ifndef OVERRIDE_fft_downshift |
565 | | static void fft_downshift(kiss_fft_cpx *x, int N, int *total, int step) { |
566 | | int shift; |
567 | | shift = IMIN(step, *total); |
568 | | *total -= shift; |
569 | | if (shift == 1) { |
570 | | int i; |
571 | | for (i=0;i<N;i++) { |
572 | | x[i].r = SHR32(x[i].r, 1); |
573 | | x[i].i = SHR32(x[i].i, 1); |
574 | | } |
575 | | } else if (shift>0) { |
576 | | int i; |
577 | | for (i=0;i<N;i++) { |
578 | | x[i].r = PSHR32(x[i].r, shift); |
579 | | x[i].i = PSHR32(x[i].i, shift); |
580 | | } |
581 | | } |
582 | | } |
583 | | #endif /* OVERRIDE_fft_downshift */ |
584 | | #else |
585 | | #define fft_downshift(x, N, total, step) |
586 | | #endif |
587 | | |
588 | | void opus_fft_impl(const kiss_fft_state *st,kiss_fft_cpx *fout ARG_FIXED(int downshift)) |
589 | 325k | { |
590 | 325k | int m2, m; |
591 | 325k | int p; |
592 | 325k | int L; |
593 | 325k | int fstride[MAXFACTORS]; |
594 | 325k | int i; |
595 | 325k | int shift; |
596 | | |
597 | | /* st->shift can be -1 */ |
598 | 325k | shift = st->shift>0 ? st->shift : 0; |
599 | | |
600 | 325k | fstride[0] = 1; |
601 | 325k | L=0; |
602 | 1.27M | do { |
603 | 1.27M | p = st->factors[2*L]; |
604 | 1.27M | m = st->factors[2*L+1]; |
605 | 1.27M | fstride[L+1] = fstride[L]*p; |
606 | 1.27M | L++; |
607 | 1.27M | } while(m!=1); |
608 | 325k | m = st->factors[2*L-1]; |
609 | 1.59M | for (i=L-1;i>=0;i--) |
610 | 1.27M | { |
611 | 1.27M | if (i!=0) |
612 | 945k | m2 = st->factors[2*i-1]; |
613 | 325k | else |
614 | 325k | m2 = 1; |
615 | 1.27M | switch (st->factors[2*i]) |
616 | 1.27M | { |
617 | 77.9k | case 2: |
618 | 77.9k | fft_downshift(fout, st->nfft, &downshift, 1); |
619 | 77.9k | kf_bfly2(fout, m, fstride[i]); |
620 | 77.9k | break; |
621 | 542k | case 4: |
622 | 542k | fft_downshift(fout, st->nfft, &downshift, 2); |
623 | 542k | kf_bfly4(fout,fstride[i]<<shift,st,m, fstride[i], m2); |
624 | 542k | break; |
625 | 0 | #ifndef RADIX_TWO_ONLY |
626 | 325k | case 3: |
627 | 325k | fft_downshift(fout, st->nfft, &downshift, 2); |
628 | 325k | kf_bfly3(fout,fstride[i]<<shift,st,m, fstride[i], m2); |
629 | 325k | break; |
630 | 325k | case 5: |
631 | 325k | fft_downshift(fout, st->nfft, &downshift, 3); |
632 | 325k | kf_bfly5(fout,fstride[i]<<shift,st,m, fstride[i], m2); |
633 | 325k | break; |
634 | 1.27M | #endif |
635 | 1.27M | } |
636 | 1.27M | m = m2; |
637 | 1.27M | } |
638 | 325k | fft_downshift(fout, st->nfft, &downshift, downshift); |
639 | 325k | } |
640 | | #endif /* !ENABLE_PFA || OPUS_CUSTOM || ENABLE_DEEP_PLC */ |
641 | | #endif |
642 | | |
643 | | void opus_fft_c(const kiss_fft_state *st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout) |
644 | 0 | { |
645 | | #if defined(ENABLE_PFA) |
646 | | if (st->nfft == 60 || st->nfft == 120 || st->nfft == 240 || st->nfft == 480 |
647 | | #if defined(ENABLE_QEXT) |
648 | | || st->nfft == 960 |
649 | | #endif |
650 | | ) |
651 | | { |
652 | | int i; |
653 | | celt_coef scale = st->scale; |
654 | | VARDECL(kiss_fft_cpx, tmp); |
655 | | SAVE_STACK; |
656 | | ALLOC(tmp, st->nfft, kiss_fft_cpx); |
657 | | for (i = 0; i < st->nfft; i++) { |
658 | | #ifndef FIXED_POINT |
659 | | tmp[i].r = fin[i].r * scale; |
660 | | tmp[i].i = fin[i].i * scale; |
661 | | #else |
662 | | tmp[i].r = S_MUL2(fin[i].r, scale); |
663 | | tmp[i].i = S_MUL2(fin[i].i, scale); |
664 | | #endif |
665 | | } |
666 | | opus_fft_pfa_c(st, tmp, fout ARG_FIXED(st->scale_shift - 1)); |
667 | | RESTORE_STACK; |
668 | | return; |
669 | | } |
670 | | #if !defined(OPUS_CUSTOM) && !defined(ENABLE_DEEP_PLC) |
671 | | celt_assert2(0, "PFA FFT called with unsupported size in non-custom mode"); |
672 | | #endif |
673 | | #endif |
674 | |
|
675 | 0 | #if !defined(ENABLE_PFA) || defined(OPUS_CUSTOM) || defined(ENABLE_DEEP_PLC) |
676 | 0 | int i; |
677 | 0 | celt_coef scale; |
678 | | #ifdef FIXED_POINT |
679 | | /* Allows us to scale with MULT16_32_Q16(), which is faster than |
680 | | MULT16_32_Q15() on ARM. */ |
681 | | int scale_shift = st->scale_shift-1; |
682 | | #endif |
683 | 0 | scale = st->scale; |
684 | |
|
685 | 0 | celt_assert2 (fin != fout, "In-place FFT not supported"); |
686 | | /* Bit-reverse the input */ |
687 | 0 | for (i=0;i<st->nfft;i++) |
688 | 0 | { |
689 | 0 | kiss_fft_cpx x = fin[i]; |
690 | 0 | fout[st->bitrev[i]].r = S_MUL2(x.r, scale); |
691 | 0 | fout[st->bitrev[i]].i = S_MUL2(x.i, scale); |
692 | 0 | } |
693 | 0 | opus_fft_impl(st, fout ARG_FIXED(scale_shift)); |
694 | 0 | #endif |
695 | 0 | } |
696 | | |
697 | | |
698 | | void opus_ifft_c(const kiss_fft_state *st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout) |
699 | 0 | { |
700 | | #if defined(ENABLE_PFA) |
701 | | if (st->nfft == 60 || st->nfft == 120 || st->nfft == 240 || st->nfft == 480 |
702 | | #if defined(ENABLE_QEXT) |
703 | | || st->nfft == 960 |
704 | | #endif |
705 | | ) |
706 | | { |
707 | | opus_ifft_pfa_c(st, fin, fout ARG_FIXED(0)); |
708 | | return; |
709 | | } |
710 | | #if !defined(OPUS_CUSTOM) && !defined(ENABLE_DEEP_PLC) |
711 | | celt_assert2(0, "PFA IFFT called with unsupported size in non-custom mode"); |
712 | | #endif |
713 | | #endif |
714 | |
|
715 | 0 | #if !defined(ENABLE_PFA) || defined(OPUS_CUSTOM) || defined(ENABLE_DEEP_PLC) |
716 | 0 | int i; |
717 | 0 | celt_assert2 (fin != fout, "In-place FFT not supported"); |
718 | | /* Bit-reverse the input */ |
719 | 0 | for (i=0;i<st->nfft;i++) |
720 | 0 | fout[st->bitrev[i]] = fin[i]; |
721 | 0 | for (i=0;i<st->nfft;i++) |
722 | 0 | fout[i].i = -fout[i].i; |
723 | 0 | opus_fft_impl(st, fout ARG_FIXED(0)); |
724 | 0 | for (i=0;i<st->nfft;i++) |
725 | 0 | fout[i].i = -fout[i].i; |
726 | 0 | #endif |
727 | 0 | } |