diff options
Diffstat (limited to 'apps/codecs/libwmapro/fft.c')
-rw-r--r-- | apps/codecs/libwmapro/fft.c | 385 |
1 files changed, 0 insertions, 385 deletions
diff --git a/apps/codecs/libwmapro/fft.c b/apps/codecs/libwmapro/fft.c deleted file mode 100644 index 80dd35b2b6..0000000000 --- a/apps/codecs/libwmapro/fft.c +++ /dev/null | |||
@@ -1,385 +0,0 @@ | |||
1 | /* | ||
2 | * FFT/IFFT transforms | ||
3 | * Copyright (c) 2008 Loren Merritt | ||
4 | * Copyright (c) 2002 Fabrice Bellard | ||
5 | * Partly based on libdjbfft by D. J. Bernstein | ||
6 | * | ||
7 | * This file is part of FFmpeg. | ||
8 | * | ||
9 | * FFmpeg is free software; you can redistribute it and/or | ||
10 | * modify it under the terms of the GNU Lesser General Public | ||
11 | * License as published by the Free Software Foundation; either | ||
12 | * version 2.1 of the License, or (at your option) any later version. | ||
13 | * | ||
14 | * FFmpeg is distributed in the hope that it will be useful, | ||
15 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
16 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | ||
17 | * Lesser General Public License for more details. | ||
18 | * | ||
19 | * You should have received a copy of the GNU Lesser General Public | ||
20 | * License along with FFmpeg; if not, write to the Free Software | ||
21 | * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA | ||
22 | */ | ||
23 | |||
24 | /** | ||
25 | * @file libavcodec/fft.c | ||
26 | * FFT/IFFT transforms. | ||
27 | */ | ||
28 | |||
29 | #include "dsputil.h" | ||
30 | #include "fft.h" | ||
31 | #ifndef M_PI | ||
32 | #define M_PI 3.14159265358979323846 /* pi */ | ||
33 | #endif | ||
34 | #ifndef M_SQRT1_2 | ||
35 | #define M_SQRT1_2 0.70710678118654752440 /* 1/sqrt(2) */ | ||
36 | #endif | ||
37 | #ifndef M_SQRT2 | ||
38 | #define M_SQRT2 1.41421356237309504880 /* sqrt(2) */ | ||
39 | #endif | ||
40 | |||
41 | /* cos(2*pi*x/n) for 0<=x<=n/4, followed by its reverse */ | ||
42 | DECLARE_ALIGNED_16(FFTSample, ff_cos_16[8]); | ||
43 | DECLARE_ALIGNED_16(FFTSample, ff_cos_32[16]); | ||
44 | DECLARE_ALIGNED_16(FFTSample, ff_cos_64[32]); | ||
45 | DECLARE_ALIGNED_16(FFTSample, ff_cos_128[64]); | ||
46 | DECLARE_ALIGNED_16(FFTSample, ff_cos_256[128]); | ||
47 | DECLARE_ALIGNED_16(FFTSample, ff_cos_512[256]); | ||
48 | DECLARE_ALIGNED_16(FFTSample, ff_cos_1024[512]); | ||
49 | DECLARE_ALIGNED_16(FFTSample, ff_cos_2048[1024]); | ||
50 | DECLARE_ALIGNED_16(FFTSample, ff_cos_4096[2048]); | ||
51 | DECLARE_ALIGNED_16(FFTSample, ff_cos_8192[4096]); | ||
52 | DECLARE_ALIGNED_16(FFTSample, ff_cos_16384[8192]); | ||
53 | DECLARE_ALIGNED_16(FFTSample, ff_cos_32768[16384]); | ||
54 | DECLARE_ALIGNED_16(FFTSample, ff_cos_65536[32768]); | ||
55 | |||
56 | FFTSample *ff_cos_tabs[] = { | ||
57 | ff_cos_16, ff_cos_32, ff_cos_64, ff_cos_128, ff_cos_256, ff_cos_512, ff_cos_1024, | ||
58 | ff_cos_2048, ff_cos_4096, ff_cos_8192, ff_cos_16384, ff_cos_32768, ff_cos_65536, | ||
59 | }; | ||
60 | |||
61 | static int split_radix_permutation(int i, int n, int inverse) | ||
62 | { | ||
63 | int m; | ||
64 | if(n <= 2) return i&1; | ||
65 | m = n >> 1; | ||
66 | if(!(i&m)) return split_radix_permutation(i, m, inverse)*2; | ||
67 | m >>= 1; | ||
68 | if(inverse == !(i&m)) return split_radix_permutation(i, m, inverse)*4 + 1; | ||
69 | else return split_radix_permutation(i, m, inverse)*4 - 1; | ||
70 | } | ||
71 | |||
72 | av_cold int ff_fft_init(FFTContext *s, int nbits, int inverse) | ||
73 | { | ||
74 | int i, j, m, n; | ||
75 | float alpha, c1, s1, s2; | ||
76 | int split_radix = 1; | ||
77 | int av_unused has_vectors; | ||
78 | |||
79 | if (nbits < 2 || nbits > 16) | ||
80 | goto fail; | ||
81 | s->nbits = nbits; | ||
82 | n = 1 << nbits; | ||
83 | |||
84 | s->tmp_buf = NULL; | ||
85 | s->exptab = av_malloc((n / 2) * sizeof(FFTComplex)); | ||
86 | if (!s->exptab) | ||
87 | goto fail; | ||
88 | s->revtab = av_malloc(n * sizeof(uint16_t)); | ||
89 | if (!s->revtab) | ||
90 | goto fail; | ||
91 | s->inverse = inverse; | ||
92 | |||
93 | s2 = inverse ? 1.0 : -1.0; | ||
94 | |||
95 | s->fft_permute = ff_fft_permute_c; | ||
96 | s->fft_calc = ff_fft_calc_c; | ||
97 | s->imdct_calc = ff_imdct_calc_c; | ||
98 | s->imdct_half = ff_imdct_half_c; | ||
99 | s->exptab1 = NULL; | ||
100 | |||
101 | #if HAVE_MMX && HAVE_YASM | ||
102 | has_vectors = mm_support(); | ||
103 | if (has_vectors & FF_MM_SSE && HAVE_SSE) { | ||
104 | /* SSE for P3/P4/K8 */ | ||
105 | s->imdct_calc = ff_imdct_calc_sse; | ||
106 | s->imdct_half = ff_imdct_half_sse; | ||
107 | s->fft_permute = ff_fft_permute_sse; | ||
108 | s->fft_calc = ff_fft_calc_sse; | ||
109 | } else if (has_vectors & FF_MM_3DNOWEXT && HAVE_AMD3DNOWEXT) { | ||
110 | /* 3DNowEx for K7 */ | ||
111 | s->imdct_calc = ff_imdct_calc_3dn2; | ||
112 | s->imdct_half = ff_imdct_half_3dn2; | ||
113 | s->fft_calc = ff_fft_calc_3dn2; | ||
114 | } else if (has_vectors & FF_MM_3DNOW && HAVE_AMD3DNOW) { | ||
115 | /* 3DNow! for K6-2/3 */ | ||
116 | s->imdct_calc = ff_imdct_calc_3dn; | ||
117 | s->imdct_half = ff_imdct_half_3dn; | ||
118 | s->fft_calc = ff_fft_calc_3dn; | ||
119 | } | ||
120 | #elif HAVE_ALTIVEC | ||
121 | has_vectors = mm_support(); | ||
122 | if (has_vectors & FF_MM_ALTIVEC) { | ||
123 | s->fft_calc = ff_fft_calc_altivec; | ||
124 | split_radix = 0; | ||
125 | } | ||
126 | #endif | ||
127 | |||
128 | if (split_radix) { | ||
129 | for(j=4; j<=nbits; j++) { | ||
130 | int m = 1<<j; | ||
131 | double freq = 2*M_PI/m; | ||
132 | FFTSample *tab = ff_cos_tabs[j-4]; | ||
133 | for(i=0; i<=m/4; i++) | ||
134 | tab[i] = cos(i*freq); | ||
135 | for(i=1; i<m/4; i++) | ||
136 | tab[m/2-i] = tab[i]; | ||
137 | } | ||
138 | for(i=0; i<n; i++) | ||
139 | s->revtab[-split_radix_permutation(i, n, s->inverse) & (n-1)] = i; | ||
140 | s->tmp_buf = av_malloc(n * sizeof(FFTComplex)); | ||
141 | } else { | ||
142 | int np, nblocks, np2, l; | ||
143 | FFTComplex *q; | ||
144 | |||
145 | for(i=0; i<(n/2); i++) { | ||
146 | alpha = 2 * M_PI * (float)i / (float)n; | ||
147 | c1 = cos(alpha); | ||
148 | s1 = sin(alpha) * s2; | ||
149 | s->exptab[i].re = c1; | ||
150 | s->exptab[i].im = s1; | ||
151 | } | ||
152 | |||
153 | np = 1 << nbits; | ||
154 | nblocks = np >> 3; | ||
155 | np2 = np >> 1; | ||
156 | s->exptab1 = av_malloc(np * 2 * sizeof(FFTComplex)); | ||
157 | if (!s->exptab1) | ||
158 | goto fail; | ||
159 | q = s->exptab1; | ||
160 | do { | ||
161 | for(l = 0; l < np2; l += 2 * nblocks) { | ||
162 | *q++ = s->exptab[l]; | ||
163 | *q++ = s->exptab[l + nblocks]; | ||
164 | |||
165 | q->re = -s->exptab[l].im; | ||
166 | q->im = s->exptab[l].re; | ||
167 | q++; | ||
168 | q->re = -s->exptab[l + nblocks].im; | ||
169 | q->im = s->exptab[l + nblocks].re; | ||
170 | q++; | ||
171 | } | ||
172 | nblocks = nblocks >> 1; | ||
173 | } while (nblocks != 0); | ||
174 | av_freep(&s->exptab); | ||
175 | |||
176 | /* compute bit reverse table */ | ||
177 | for(i=0;i<n;i++) { | ||
178 | m=0; | ||
179 | for(j=0;j<nbits;j++) { | ||
180 | m |= ((i >> j) & 1) << (nbits-j-1); | ||
181 | } | ||
182 | s->revtab[i]=m; | ||
183 | } | ||
184 | } | ||
185 | |||
186 | return 0; | ||
187 | fail: | ||
188 | av_freep(&s->revtab); | ||
189 | av_freep(&s->exptab); | ||
190 | av_freep(&s->exptab1); | ||
191 | av_freep(&s->tmp_buf); | ||
192 | return -1; | ||
193 | } | ||
194 | |||
195 | void ff_fft_permute_c(FFTContext *s, FFTComplex *z) | ||
196 | { | ||
197 | int j, k, np; | ||
198 | FFTComplex tmp; | ||
199 | const uint16_t *revtab = s->revtab; | ||
200 | np = 1 << s->nbits; | ||
201 | |||
202 | if (s->tmp_buf) { | ||
203 | /* TODO: handle split-radix permute in a more optimal way, probably in-place */ | ||
204 | for(j=0;j<np;j++) s->tmp_buf[revtab[j]] = z[j]; | ||
205 | memcpy(z, s->tmp_buf, np * sizeof(FFTComplex)); | ||
206 | return; | ||
207 | } | ||
208 | |||
209 | /* reverse */ | ||
210 | for(j=0;j<np;j++) { | ||
211 | k = revtab[j]; | ||
212 | if (k < j) { | ||
213 | tmp = z[k]; | ||
214 | z[k] = z[j]; | ||
215 | z[j] = tmp; | ||
216 | } | ||
217 | } | ||
218 | } | ||
219 | |||
220 | av_cold void ff_fft_end(FFTContext *s) | ||
221 | { | ||
222 | av_freep(&s->revtab); | ||
223 | av_freep(&s->exptab); | ||
224 | av_freep(&s->exptab1); | ||
225 | av_freep(&s->tmp_buf); | ||
226 | } | ||
227 | |||
228 | #define sqrthalf (float)M_SQRT1_2 | ||
229 | |||
230 | #define BF(x,y,a,b) {\ | ||
231 | x = a - b;\ | ||
232 | y = a + b;\ | ||
233 | } | ||
234 | |||
235 | #define BUTTERFLIES(a0,a1,a2,a3) {\ | ||
236 | BF(t3, t5, t5, t1);\ | ||
237 | BF(a2.re, a0.re, a0.re, t5);\ | ||
238 | BF(a3.im, a1.im, a1.im, t3);\ | ||
239 | BF(t4, t6, t2, t6);\ | ||
240 | BF(a3.re, a1.re, a1.re, t4);\ | ||
241 | BF(a2.im, a0.im, a0.im, t6);\ | ||
242 | } | ||
243 | |||
244 | // force loading all the inputs before storing any. | ||
245 | // this is slightly slower for small data, but avoids store->load aliasing | ||
246 | // for addresses separated by large powers of 2. | ||
247 | #define BUTTERFLIES_BIG(a0,a1,a2,a3) {\ | ||
248 | FFTSample r0=a0.re, i0=a0.im, r1=a1.re, i1=a1.im;\ | ||
249 | BF(t3, t5, t5, t1);\ | ||
250 | BF(a2.re, a0.re, r0, t5);\ | ||
251 | BF(a3.im, a1.im, i1, t3);\ | ||
252 | BF(t4, t6, t2, t6);\ | ||
253 | BF(a3.re, a1.re, r1, t4);\ | ||
254 | BF(a2.im, a0.im, i0, t6);\ | ||
255 | } | ||
256 | |||
257 | #define TRANSFORM(a0,a1,a2,a3,wre,wim) {\ | ||
258 | t1 = a2.re * wre + a2.im * wim;\ | ||
259 | t2 = a2.im * wre - a2.re * wim;\ | ||
260 | t5 = a3.re * wre - a3.im * wim;\ | ||
261 | t6 = a3.im * wre + a3.re * wim;\ | ||
262 | BUTTERFLIES(a0,a1,a2,a3)\ | ||
263 | } | ||
264 | |||
265 | #define TRANSFORM_ZERO(a0,a1,a2,a3) {\ | ||
266 | t1 = a2.re;\ | ||
267 | t2 = a2.im;\ | ||
268 | t5 = a3.re;\ | ||
269 | t6 = a3.im;\ | ||
270 | BUTTERFLIES(a0,a1,a2,a3)\ | ||
271 | } | ||
272 | |||
273 | /* z[0...8n-1], w[1...2n-1] */ | ||
274 | #define PASS(name)\ | ||
275 | static void name(FFTComplex *z, const FFTSample *wre, unsigned int n)\ | ||
276 | {\ | ||
277 | FFTSample t1, t2, t3, t4, t5, t6;\ | ||
278 | int o1 = 2*n;\ | ||
279 | int o2 = 4*n;\ | ||
280 | int o3 = 6*n;\ | ||
281 | const FFTSample *wim = wre+o1;\ | ||
282 | n--;\ | ||
283 | \ | ||
284 | TRANSFORM_ZERO(z[0],z[o1],z[o2],z[o3]);\ | ||
285 | TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\ | ||
286 | do {\ | ||
287 | z += 2;\ | ||
288 | wre += 2;\ | ||
289 | wim -= 2;\ | ||
290 | TRANSFORM(z[0],z[o1],z[o2],z[o3],wre[0],wim[0]);\ | ||
291 | TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\ | ||
292 | } while(--n);\ | ||
293 | } | ||
294 | |||
295 | PASS(pass) | ||
296 | #undef BUTTERFLIES | ||
297 | #define BUTTERFLIES BUTTERFLIES_BIG | ||
298 | PASS(pass_big) | ||
299 | |||
300 | #define DECL_FFT(n,n2,n4)\ | ||
301 | static void fft##n(FFTComplex *z)\ | ||
302 | {\ | ||
303 | fft##n2(z);\ | ||
304 | fft##n4(z+n4*2);\ | ||
305 | fft##n4(z+n4*3);\ | ||
306 | pass(z,ff_cos_##n,n4/2);\ | ||
307 | } | ||
308 | |||
309 | static void fft4(FFTComplex *z) | ||
310 | { | ||
311 | FFTSample t1, t2, t3, t4, t5, t6, t7, t8; | ||
312 | |||
313 | BF(t3, t1, z[0].re, z[1].re); | ||
314 | BF(t8, t6, z[3].re, z[2].re); | ||
315 | BF(z[2].re, z[0].re, t1, t6); | ||
316 | BF(t4, t2, z[0].im, z[1].im); | ||
317 | BF(t7, t5, z[2].im, z[3].im); | ||
318 | BF(z[3].im, z[1].im, t4, t8); | ||
319 | BF(z[3].re, z[1].re, t3, t7); | ||
320 | BF(z[2].im, z[0].im, t2, t5); | ||
321 | } | ||
322 | |||
323 | static void fft8(FFTComplex *z) | ||
324 | { | ||
325 | FFTSample t1, t2, t3, t4, t5, t6, t7, t8; | ||
326 | |||
327 | fft4(z); | ||
328 | |||
329 | BF(t1, z[5].re, z[4].re, -z[5].re); | ||
330 | BF(t2, z[5].im, z[4].im, -z[5].im); | ||
331 | BF(t3, z[7].re, z[6].re, -z[7].re); | ||
332 | BF(t4, z[7].im, z[6].im, -z[7].im); | ||
333 | BF(t8, t1, t3, t1); | ||
334 | BF(t7, t2, t2, t4); | ||
335 | BF(z[4].re, z[0].re, z[0].re, t1); | ||
336 | BF(z[4].im, z[0].im, z[0].im, t2); | ||
337 | BF(z[6].re, z[2].re, z[2].re, t7); | ||
338 | BF(z[6].im, z[2].im, z[2].im, t8); | ||
339 | |||
340 | TRANSFORM(z[1],z[3],z[5],z[7],sqrthalf,sqrthalf); | ||
341 | } | ||
342 | |||
343 | #if !CONFIG_SMALL | ||
344 | static void fft16(FFTComplex *z) | ||
345 | { | ||
346 | FFTSample t1, t2, t3, t4, t5, t6; | ||
347 | |||
348 | fft8(z); | ||
349 | fft4(z+8); | ||
350 | fft4(z+12); | ||
351 | |||
352 | TRANSFORM_ZERO(z[0],z[4],z[8],z[12]); | ||
353 | TRANSFORM(z[2],z[6],z[10],z[14],sqrthalf,sqrthalf); | ||
354 | TRANSFORM(z[1],z[5],z[9],z[13],ff_cos_16[1],ff_cos_16[3]); | ||
355 | TRANSFORM(z[3],z[7],z[11],z[15],ff_cos_16[3],ff_cos_16[1]); | ||
356 | } | ||
357 | #else | ||
358 | DECL_FFT(16,8,4) | ||
359 | #endif | ||
360 | DECL_FFT(32,16,8) | ||
361 | DECL_FFT(64,32,16) | ||
362 | DECL_FFT(128,64,32) | ||
363 | DECL_FFT(256,128,64) | ||
364 | DECL_FFT(512,256,128) | ||
365 | #if !CONFIG_SMALL | ||
366 | #define pass pass_big | ||
367 | #endif | ||
368 | DECL_FFT(1024,512,256) | ||
369 | DECL_FFT(2048,1024,512) | ||
370 | DECL_FFT(4096,2048,1024) | ||
371 | DECL_FFT(8192,4096,2048) | ||
372 | DECL_FFT(16384,8192,4096) | ||
373 | DECL_FFT(32768,16384,8192) | ||
374 | DECL_FFT(65536,32768,16384) | ||
375 | |||
376 | static void (*fft_dispatch[])(FFTComplex*) = { | ||
377 | fft4, fft8, fft16, fft32, fft64, fft128, fft256, fft512, fft1024, | ||
378 | fft2048, fft4096, fft8192, fft16384, fft32768, fft65536, | ||
379 | }; | ||
380 | |||
381 | void ff_fft_calc_c(FFTContext *s, FFTComplex *z) | ||
382 | { | ||
383 | fft_dispatch[s->nbits-2](z); | ||
384 | } | ||
385 | |||