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