Mercurial > libavcodec.hg
annotate fft.c @ 11352:6e0af2cfdcfe libavcodec
Do MC and IDCT in coding (hilbert) order
This increases the slice size to 64 pixels, due to having to decode an
entire chroma superblock row per slice.
This can be up to 6% slower depending on clip and CPU, but is necessary
for future optimizations that gain significantly more than was lost.
author | conrad |
---|---|
date | Wed, 03 Mar 2010 23:27:40 +0000 |
parents | e45c852b6820 |
children | 4b3da727d832 |
rev | line source |
---|---|
781 | 1 /* |
2 * FFT/IFFT transforms | |
7542 | 3 * Copyright (c) 2008 Loren Merritt |
8629
04423b2f6e0b
cosmetics: Remove pointless period after copyright statement non-sentences.
diego
parents:
8590
diff
changeset
|
4 * Copyright (c) 2002 Fabrice Bellard |
7542 | 5 * Partly based on libdjbfft by D. J. Bernstein |
781 | 6 * |
3947
c8c591fe26f8
Change license headers to say 'FFmpeg' instead of 'this program/this library'
diego
parents:
3746
diff
changeset
|
7 * This file is part of FFmpeg. |
c8c591fe26f8
Change license headers to say 'FFmpeg' instead of 'this program/this library'
diego
parents:
3746
diff
changeset
|
8 * |
c8c591fe26f8
Change license headers to say 'FFmpeg' instead of 'this program/this library'
diego
parents:
3746
diff
changeset
|
9 * FFmpeg is free software; you can redistribute it and/or |
781 | 10 * modify it under the terms of the GNU Lesser General Public |
11 * License as published by the Free Software Foundation; either | |
3947
c8c591fe26f8
Change license headers to say 'FFmpeg' instead of 'this program/this library'
diego
parents:
3746
diff
changeset
|
12 * version 2.1 of the License, or (at your option) any later version. |
781 | 13 * |
3947
c8c591fe26f8
Change license headers to say 'FFmpeg' instead of 'this program/this library'
diego
parents:
3746
diff
changeset
|
14 * FFmpeg is distributed in the hope that it will be useful, |
781 | 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 | |
3947
c8c591fe26f8
Change license headers to say 'FFmpeg' instead of 'this program/this library'
diego
parents:
3746
diff
changeset
|
20 * License along with FFmpeg; if not, write to the Free Software |
3036
0b546eab515d
Update licensing information: The FSF changed postal address.
diego
parents:
2979
diff
changeset
|
21 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
781 | 22 */ |
1106 | 23 |
24 /** | |
8718
e9d9d946f213
Use full internal pathname in doxygen @file directives.
diego
parents:
8694
diff
changeset
|
25 * @file libavcodec/fft.c |
1106 | 26 * FFT/IFFT transforms. |
27 */ | |
28 | |
781 | 29 #include "dsputil.h" |
30 | |
10407
57acce8b1380
Move/add COSTABLE/SINTABLE macros to dsputil to add extern definitions
reimar
parents:
10400
diff
changeset
|
31 /* cos(2*pi*x/n) for 0<=x<=n/4, followed by its reverse */ |
57acce8b1380
Move/add COSTABLE/SINTABLE macros to dsputil to add extern definitions
reimar
parents:
10400
diff
changeset
|
32 #if !CONFIG_HARDCODED_TABLES |
57acce8b1380
Move/add COSTABLE/SINTABLE macros to dsputil to add extern definitions
reimar
parents:
10400
diff
changeset
|
33 COSTABLE(16); |
57acce8b1380
Move/add COSTABLE/SINTABLE macros to dsputil to add extern definitions
reimar
parents:
10400
diff
changeset
|
34 COSTABLE(32); |
57acce8b1380
Move/add COSTABLE/SINTABLE macros to dsputil to add extern definitions
reimar
parents:
10400
diff
changeset
|
35 COSTABLE(64); |
57acce8b1380
Move/add COSTABLE/SINTABLE macros to dsputil to add extern definitions
reimar
parents:
10400
diff
changeset
|
36 COSTABLE(128); |
57acce8b1380
Move/add COSTABLE/SINTABLE macros to dsputil to add extern definitions
reimar
parents:
10400
diff
changeset
|
37 COSTABLE(256); |
57acce8b1380
Move/add COSTABLE/SINTABLE macros to dsputil to add extern definitions
reimar
parents:
10400
diff
changeset
|
38 COSTABLE(512); |
57acce8b1380
Move/add COSTABLE/SINTABLE macros to dsputil to add extern definitions
reimar
parents:
10400
diff
changeset
|
39 COSTABLE(1024); |
57acce8b1380
Move/add COSTABLE/SINTABLE macros to dsputil to add extern definitions
reimar
parents:
10400
diff
changeset
|
40 COSTABLE(2048); |
57acce8b1380
Move/add COSTABLE/SINTABLE macros to dsputil to add extern definitions
reimar
parents:
10400
diff
changeset
|
41 COSTABLE(4096); |
57acce8b1380
Move/add COSTABLE/SINTABLE macros to dsputil to add extern definitions
reimar
parents:
10400
diff
changeset
|
42 COSTABLE(8192); |
57acce8b1380
Move/add COSTABLE/SINTABLE macros to dsputil to add extern definitions
reimar
parents:
10400
diff
changeset
|
43 COSTABLE(16384); |
57acce8b1380
Move/add COSTABLE/SINTABLE macros to dsputil to add extern definitions
reimar
parents:
10400
diff
changeset
|
44 COSTABLE(32768); |
57acce8b1380
Move/add COSTABLE/SINTABLE macros to dsputil to add extern definitions
reimar
parents:
10400
diff
changeset
|
45 COSTABLE(65536); |
10400
866dffa620d1
Use hardcoded instead of runtime-calculated ff_cos_* tables if
reimar
parents:
10391
diff
changeset
|
46 #endif |
10407
57acce8b1380
Move/add COSTABLE/SINTABLE macros to dsputil to add extern definitions
reimar
parents:
10400
diff
changeset
|
47 COSTABLE_CONST FFTSample * const ff_cos_tabs[] = { |
10492
63910f7ba293
Pad ff_cos_tabs and ff_sin_tabs so that index n points to the table for n bits.
reimar
parents:
10407
diff
changeset
|
48 NULL, NULL, NULL, NULL, |
7542 | 49 ff_cos_16, ff_cos_32, ff_cos_64, ff_cos_128, ff_cos_256, ff_cos_512, ff_cos_1024, |
50 ff_cos_2048, ff_cos_4096, ff_cos_8192, ff_cos_16384, ff_cos_32768, ff_cos_65536, | |
51 }; | |
52 | |
53 static int split_radix_permutation(int i, int n, int inverse) | |
54 { | |
55 int m; | |
56 if(n <= 2) return i&1; | |
57 m = n >> 1; | |
58 if(!(i&m)) return split_radix_permutation(i, m, inverse)*2; | |
59 m >>= 1; | |
60 if(inverse == !(i&m)) return split_radix_permutation(i, m, inverse)*4 + 1; | |
61 else return split_radix_permutation(i, m, inverse)*4 - 1; | |
62 } | |
63 | |
10496
74b0c1a0851e
Add ff_init_ff_cos_tabs function and use it in rdft.c to ensure that the
reimar
parents:
10492
diff
changeset
|
64 av_cold void ff_init_ff_cos_tabs(int index) |
74b0c1a0851e
Add ff_init_ff_cos_tabs function and use it in rdft.c to ensure that the
reimar
parents:
10492
diff
changeset
|
65 { |
74b0c1a0851e
Add ff_init_ff_cos_tabs function and use it in rdft.c to ensure that the
reimar
parents:
10492
diff
changeset
|
66 #if !CONFIG_HARDCODED_TABLES |
74b0c1a0851e
Add ff_init_ff_cos_tabs function and use it in rdft.c to ensure that the
reimar
parents:
10492
diff
changeset
|
67 int i; |
74b0c1a0851e
Add ff_init_ff_cos_tabs function and use it in rdft.c to ensure that the
reimar
parents:
10492
diff
changeset
|
68 int m = 1<<index; |
74b0c1a0851e
Add ff_init_ff_cos_tabs function and use it in rdft.c to ensure that the
reimar
parents:
10492
diff
changeset
|
69 double freq = 2*M_PI/m; |
74b0c1a0851e
Add ff_init_ff_cos_tabs function and use it in rdft.c to ensure that the
reimar
parents:
10492
diff
changeset
|
70 FFTSample *tab = ff_cos_tabs[index]; |
74b0c1a0851e
Add ff_init_ff_cos_tabs function and use it in rdft.c to ensure that the
reimar
parents:
10492
diff
changeset
|
71 for(i=0; i<=m/4; i++) |
74b0c1a0851e
Add ff_init_ff_cos_tabs function and use it in rdft.c to ensure that the
reimar
parents:
10492
diff
changeset
|
72 tab[i] = cos(i*freq); |
74b0c1a0851e
Add ff_init_ff_cos_tabs function and use it in rdft.c to ensure that the
reimar
parents:
10492
diff
changeset
|
73 for(i=1; i<m/4; i++) |
74b0c1a0851e
Add ff_init_ff_cos_tabs function and use it in rdft.c to ensure that the
reimar
parents:
10492
diff
changeset
|
74 tab[m/2-i] = tab[i]; |
74b0c1a0851e
Add ff_init_ff_cos_tabs function and use it in rdft.c to ensure that the
reimar
parents:
10492
diff
changeset
|
75 #endif |
74b0c1a0851e
Add ff_init_ff_cos_tabs function and use it in rdft.c to ensure that the
reimar
parents:
10492
diff
changeset
|
76 } |
74b0c1a0851e
Add ff_init_ff_cos_tabs function and use it in rdft.c to ensure that the
reimar
parents:
10492
diff
changeset
|
77 |
8637 | 78 av_cold int ff_fft_init(FFTContext *s, int nbits, int inverse) |
781 | 79 { |
80 int i, j, m, n; | |
81 float alpha, c1, s1, s2; | |
6503 | 82 int av_unused has_vectors; |
2967 | 83 |
7542 | 84 if (nbits < 2 || nbits > 16) |
85 goto fail; | |
781 | 86 s->nbits = nbits; |
87 n = 1 << nbits; | |
88 | |
7542 | 89 s->tmp_buf = NULL; |
8974 | 90 s->exptab = av_malloc((n / 2) * sizeof(FFTComplex)); |
781 | 91 if (!s->exptab) |
92 goto fail; | |
93 s->revtab = av_malloc(n * sizeof(uint16_t)); | |
94 if (!s->revtab) | |
95 goto fail; | |
96 s->inverse = inverse; | |
97 | |
98 s2 = inverse ? 1.0 : -1.0; | |
2967 | 99 |
7542 | 100 s->fft_permute = ff_fft_permute_c; |
8974 | 101 s->fft_calc = ff_fft_calc_c; |
11122 | 102 #if CONFIG_MDCT |
8974 | 103 s->imdct_calc = ff_imdct_calc_c; |
104 s->imdct_half = ff_imdct_half_c; | |
10161 | 105 s->mdct_calc = ff_mdct_calc_c; |
11122 | 106 #endif |
8974 | 107 s->exptab1 = NULL; |
10175
5cf49858179a
Move per-arch fft init bits into the corresponding subdirs
mru
parents:
10172
diff
changeset
|
108 s->split_radix = 1; |
781 | 109 |
10175
5cf49858179a
Move per-arch fft init bits into the corresponding subdirs
mru
parents:
10172
diff
changeset
|
110 if (ARCH_ARM) ff_fft_init_arm(s); |
5cf49858179a
Move per-arch fft init bits into the corresponding subdirs
mru
parents:
10172
diff
changeset
|
111 if (HAVE_ALTIVEC) ff_fft_init_altivec(s); |
5cf49858179a
Move per-arch fft init bits into the corresponding subdirs
mru
parents:
10172
diff
changeset
|
112 if (HAVE_MMX) ff_fft_init_mmx(s); |
781 | 113 |
10175
5cf49858179a
Move per-arch fft init bits into the corresponding subdirs
mru
parents:
10172
diff
changeset
|
114 if (s->split_radix) { |
7542 | 115 for(j=4; j<=nbits; j++) { |
10496
74b0c1a0851e
Add ff_init_ff_cos_tabs function and use it in rdft.c to ensure that the
reimar
parents:
10492
diff
changeset
|
116 ff_init_ff_cos_tabs(j); |
7542 | 117 } |
118 for(i=0; i<n; i++) | |
10172 | 119 s->revtab[-split_radix_permutation(i, n, s->inverse) & (n-1)] = i; |
7542 | 120 s->tmp_buf = av_malloc(n * sizeof(FFTComplex)); |
121 } else { | |
6504 | 122 int np, nblocks, np2, l; |
123 FFTComplex *q; | |
2967 | 124 |
7542 | 125 for(i=0; i<(n/2); i++) { |
126 alpha = 2 * M_PI * (float)i / (float)n; | |
127 c1 = cos(alpha); | |
128 s1 = sin(alpha) * s2; | |
129 s->exptab[i].re = c1; | |
130 s->exptab[i].im = s1; | |
131 } | |
132 | |
6504 | 133 np = 1 << nbits; |
134 nblocks = np >> 3; | |
135 np2 = np >> 1; | |
136 s->exptab1 = av_malloc(np * 2 * sizeof(FFTComplex)); | |
137 if (!s->exptab1) | |
138 goto fail; | |
139 q = s->exptab1; | |
140 do { | |
141 for(l = 0; l < np2; l += 2 * nblocks) { | |
142 *q++ = s->exptab[l]; | |
143 *q++ = s->exptab[l + nblocks]; | |
975
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
144 |
6504 | 145 q->re = -s->exptab[l].im; |
146 q->im = s->exptab[l].re; | |
147 q++; | |
148 q->re = -s->exptab[l + nblocks].im; | |
149 q->im = s->exptab[l + nblocks].re; | |
150 q++; | |
151 } | |
152 nblocks = nblocks >> 1; | |
153 } while (nblocks != 0); | |
154 av_freep(&s->exptab); | |
781 | 155 |
7543 | 156 /* compute bit reverse table */ |
157 for(i=0;i<n;i++) { | |
158 m=0; | |
159 for(j=0;j<nbits;j++) { | |
160 m |= ((i >> j) & 1) << (nbits-j-1); | |
161 } | |
162 s->revtab[i]=m; | |
781 | 163 } |
7542 | 164 } |
165 | |
781 | 166 return 0; |
167 fail: | |
168 av_freep(&s->revtab); | |
169 av_freep(&s->exptab); | |
170 av_freep(&s->exptab1); | |
7542 | 171 av_freep(&s->tmp_buf); |
781 | 172 return -1; |
173 } | |
174 | |
7542 | 175 void ff_fft_permute_c(FFTContext *s, FFTComplex *z) |
781 | 176 { |
177 int j, k, np; | |
178 FFTComplex tmp; | |
179 const uint16_t *revtab = s->revtab; | |
7542 | 180 np = 1 << s->nbits; |
181 | |
182 if (s->tmp_buf) { | |
183 /* TODO: handle split-radix permute in a more optimal way, probably in-place */ | |
184 for(j=0;j<np;j++) s->tmp_buf[revtab[j]] = z[j]; | |
185 memcpy(z, s->tmp_buf, np * sizeof(FFTComplex)); | |
186 return; | |
187 } | |
2967 | 188 |
781 | 189 /* reverse */ |
190 for(j=0;j<np;j++) { | |
191 k = revtab[j]; | |
192 if (k < j) { | |
193 tmp = z[k]; | |
194 z[k] = z[j]; | |
195 z[j] = tmp; | |
196 } | |
197 } | |
198 } | |
199 | |
8687 | 200 av_cold void ff_fft_end(FFTContext *s) |
781 | 201 { |
202 av_freep(&s->revtab); | |
203 av_freep(&s->exptab); | |
204 av_freep(&s->exptab1); | |
7542 | 205 av_freep(&s->tmp_buf); |
206 } | |
207 | |
208 #define sqrthalf (float)M_SQRT1_2 | |
209 | |
210 #define BF(x,y,a,b) {\ | |
211 x = a - b;\ | |
212 y = a + b;\ | |
213 } | |
214 | |
215 #define BUTTERFLIES(a0,a1,a2,a3) {\ | |
216 BF(t3, t5, t5, t1);\ | |
217 BF(a2.re, a0.re, a0.re, t5);\ | |
218 BF(a3.im, a1.im, a1.im, t3);\ | |
219 BF(t4, t6, t2, t6);\ | |
220 BF(a3.re, a1.re, a1.re, t4);\ | |
221 BF(a2.im, a0.im, a0.im, t6);\ | |
222 } | |
223 | |
224 // force loading all the inputs before storing any. | |
225 // this is slightly slower for small data, but avoids store->load aliasing | |
226 // for addresses separated by large powers of 2. | |
227 #define BUTTERFLIES_BIG(a0,a1,a2,a3) {\ | |
228 FFTSample r0=a0.re, i0=a0.im, r1=a1.re, i1=a1.im;\ | |
229 BF(t3, t5, t5, t1);\ | |
230 BF(a2.re, a0.re, r0, t5);\ | |
231 BF(a3.im, a1.im, i1, t3);\ | |
232 BF(t4, t6, t2, t6);\ | |
233 BF(a3.re, a1.re, r1, t4);\ | |
234 BF(a2.im, a0.im, i0, t6);\ | |
235 } | |
236 | |
237 #define TRANSFORM(a0,a1,a2,a3,wre,wim) {\ | |
238 t1 = a2.re * wre + a2.im * wim;\ | |
239 t2 = a2.im * wre - a2.re * wim;\ | |
240 t5 = a3.re * wre - a3.im * wim;\ | |
241 t6 = a3.im * wre + a3.re * wim;\ | |
242 BUTTERFLIES(a0,a1,a2,a3)\ | |
243 } | |
244 | |
245 #define TRANSFORM_ZERO(a0,a1,a2,a3) {\ | |
246 t1 = a2.re;\ | |
247 t2 = a2.im;\ | |
248 t5 = a3.re;\ | |
249 t6 = a3.im;\ | |
250 BUTTERFLIES(a0,a1,a2,a3)\ | |
251 } | |
252 | |
253 /* z[0...8n-1], w[1...2n-1] */ | |
254 #define PASS(name)\ | |
255 static void name(FFTComplex *z, const FFTSample *wre, unsigned int n)\ | |
256 {\ | |
257 FFTSample t1, t2, t3, t4, t5, t6;\ | |
258 int o1 = 2*n;\ | |
259 int o2 = 4*n;\ | |
260 int o3 = 6*n;\ | |
261 const FFTSample *wim = wre+o1;\ | |
262 n--;\ | |
263 \ | |
264 TRANSFORM_ZERO(z[0],z[o1],z[o2],z[o3]);\ | |
265 TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\ | |
266 do {\ | |
267 z += 2;\ | |
268 wre += 2;\ | |
269 wim -= 2;\ | |
270 TRANSFORM(z[0],z[o1],z[o2],z[o3],wre[0],wim[0]);\ | |
271 TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\ | |
272 } while(--n);\ | |
781 | 273 } |
274 | |
7542 | 275 PASS(pass) |
276 #undef BUTTERFLIES | |
277 #define BUTTERFLIES BUTTERFLIES_BIG | |
278 PASS(pass_big) | |
279 | |
280 #define DECL_FFT(n,n2,n4)\ | |
281 static void fft##n(FFTComplex *z)\ | |
282 {\ | |
283 fft##n2(z);\ | |
284 fft##n4(z+n4*2);\ | |
285 fft##n4(z+n4*3);\ | |
286 pass(z,ff_cos_##n,n4/2);\ | |
287 } | |
288 | |
289 static void fft4(FFTComplex *z) | |
290 { | |
291 FFTSample t1, t2, t3, t4, t5, t6, t7, t8; | |
292 | |
293 BF(t3, t1, z[0].re, z[1].re); | |
294 BF(t8, t6, z[3].re, z[2].re); | |
295 BF(z[2].re, z[0].re, t1, t6); | |
296 BF(t4, t2, z[0].im, z[1].im); | |
297 BF(t7, t5, z[2].im, z[3].im); | |
298 BF(z[3].im, z[1].im, t4, t8); | |
299 BF(z[3].re, z[1].re, t3, t7); | |
300 BF(z[2].im, z[0].im, t2, t5); | |
301 } | |
302 | |
303 static void fft8(FFTComplex *z) | |
304 { | |
305 FFTSample t1, t2, t3, t4, t5, t6, t7, t8; | |
306 | |
307 fft4(z); | |
308 | |
309 BF(t1, z[5].re, z[4].re, -z[5].re); | |
310 BF(t2, z[5].im, z[4].im, -z[5].im); | |
311 BF(t3, z[7].re, z[6].re, -z[7].re); | |
312 BF(t4, z[7].im, z[6].im, -z[7].im); | |
313 BF(t8, t1, t3, t1); | |
314 BF(t7, t2, t2, t4); | |
315 BF(z[4].re, z[0].re, z[0].re, t1); | |
316 BF(z[4].im, z[0].im, z[0].im, t2); | |
317 BF(z[6].re, z[2].re, z[2].re, t7); | |
318 BF(z[6].im, z[2].im, z[2].im, t8); | |
319 | |
320 TRANSFORM(z[1],z[3],z[5],z[7],sqrthalf,sqrthalf); | |
321 } | |
322 | |
8590 | 323 #if !CONFIG_SMALL |
7542 | 324 static void fft16(FFTComplex *z) |
325 { | |
326 FFTSample t1, t2, t3, t4, t5, t6; | |
327 | |
328 fft8(z); | |
329 fft4(z+8); | |
330 fft4(z+12); | |
331 | |
332 TRANSFORM_ZERO(z[0],z[4],z[8],z[12]); | |
333 TRANSFORM(z[2],z[6],z[10],z[14],sqrthalf,sqrthalf); | |
334 TRANSFORM(z[1],z[5],z[9],z[13],ff_cos_16[1],ff_cos_16[3]); | |
335 TRANSFORM(z[3],z[7],z[11],z[15],ff_cos_16[3],ff_cos_16[1]); | |
336 } | |
337 #else | |
338 DECL_FFT(16,8,4) | |
339 #endif | |
340 DECL_FFT(32,16,8) | |
341 DECL_FFT(64,32,16) | |
342 DECL_FFT(128,64,32) | |
343 DECL_FFT(256,128,64) | |
344 DECL_FFT(512,256,128) | |
8590 | 345 #if !CONFIG_SMALL |
7542 | 346 #define pass pass_big |
347 #endif | |
348 DECL_FFT(1024,512,256) | |
349 DECL_FFT(2048,1024,512) | |
350 DECL_FFT(4096,2048,1024) | |
351 DECL_FFT(8192,4096,2048) | |
352 DECL_FFT(16384,8192,4096) | |
353 DECL_FFT(32768,16384,8192) | |
354 DECL_FFT(65536,32768,16384) | |
355 | |
10391 | 356 static void (* const fft_dispatch[])(FFTComplex*) = { |
7542 | 357 fft4, fft8, fft16, fft32, fft64, fft128, fft256, fft512, fft1024, |
358 fft2048, fft4096, fft8192, fft16384, fft32768, fft65536, | |
359 }; | |
360 | |
361 void ff_fft_calc_c(FFTContext *s, FFTComplex *z) | |
362 { | |
363 fft_dispatch[s->nbits-2](z); | |
364 } | |
365 |