3 * Copyright (c) 2002 Fabrice Bellard.
5 * This library is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU Lesser General Public
7 * License as published by the Free Software Foundation; either
8 * version 2 of the License, or (at your option) any later version.
10 * This library is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * Lesser General Public License for more details.
15 * You should have received a copy of the GNU Lesser General Public
16 * License along with this library; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22 * FFT/IFFT transforms.
28 * The size of the FFT is 2^nbits. If inverse is TRUE, inverse FFT is
31 int ff_fft_init(FFTContext
*s
, int nbits
, int inverse
)
34 float alpha
, c1
, s1
, s2
;
39 s
->exptab
= av_malloc((n
/ 2) * sizeof(FFTComplex
));
42 s
->revtab
= av_malloc(n
* sizeof(uint16_t));
47 s2
= inverse
? 1.0 : -1.0;
49 for(i
=0;i
<(n
/2);i
++) {
50 alpha
= 2 * M_PI
* (float)i
/ (float)n
;
56 s
->fft_calc
= ff_fft_calc_c
;
59 /* compute constant table for HAVE_SSE version */
60 #if (defined(HAVE_MMX) && (defined(HAVE_BUILTIN_VECTOR) || defined(HAVE_MM3DNOW))) || defined(HAVE_ALTIVEC)
65 has_vectors
= mm_support() & (MM_3DNOW
| MM_3DNOWEXT
| MM_SSE
| MM_SSE2
);
67 #if defined(HAVE_ALTIVEC) && !defined(ALTIVEC_USE_REFERENCE_C_CODE)
68 has_vectors
= mm_support() & MM_ALTIVEC
;
71 int np
, nblocks
, np2
, l
;
77 s
->exptab1
= av_malloc(np
* 2 * sizeof(FFTComplex
));
82 for(l
= 0; l
< np2
; l
+= 2 * nblocks
) {
84 *q
++ = s
->exptab
[l
+ nblocks
];
86 q
->re
= -s
->exptab
[l
].im
;
87 q
->im
= s
->exptab
[l
].re
;
89 q
->re
= -s
->exptab
[l
+ nblocks
].im
;
90 q
->im
= s
->exptab
[l
+ nblocks
].re
;
93 nblocks
= nblocks
>> 1;
94 } while (nblocks
!= 0);
98 if (has_vectors
& MM_3DNOWEXT
)
99 /* 3DNowEx for Athlon(XP) */
100 s
->fft_calc
= ff_fft_calc_3dn2
;
101 else if (has_vectors
& MM_3DNOW
)
102 /* 3DNow! for K6-2/3 */
103 s
->fft_calc
= ff_fft_calc_3dn
;
105 #ifdef HAVE_BUILTIN_VECTOR
106 if (has_vectors
& MM_SSE2
)
108 s
->fft_calc
= ff_fft_calc_sse
;
109 else if ((has_vectors
& MM_SSE
) &&
110 s
->fft_calc
== ff_fft_calc_c
)
112 s
->fft_calc
= ff_fft_calc_sse
;
115 s
->fft_calc
= ff_fft_calc_altivec
;
121 /* compute bit reverse table */
125 for(j
=0;j
<nbits
;j
++) {
126 m
|= ((i
>> j
) & 1) << (nbits
-j
-1);
132 av_freep(&s
->revtab
);
133 av_freep(&s
->exptab
);
134 av_freep(&s
->exptab1
);
139 #define BF(pre, pim, qre, qim, pre1, pim1, qre1, qim1) \
141 FFTSample ax, ay, bx, by;\
152 #define MUL16(a,b) ((a) * (b))
154 #define CMUL(pre, pim, are, aim, bre, bim) \
156 pre = (MUL16(are, bre) - MUL16(aim, bim));\
157 pim = (MUL16(are, bim) + MUL16(bre, aim));\
161 * Do a complex FFT with the parameters defined in ff_fft_init(). The
162 * input data must be permuted before with s->revtab table. No
163 * 1.0/sqrt(n) normalization is done.
165 void ff_fft_calc_c(FFTContext
*s
, FFTComplex
*z
)
170 register FFTComplex
*p
, *q
;
171 FFTComplex
*exptab
= s
->exptab
;
173 FFTSample tmp_re
, tmp_im
;
182 BF(p
[0].re
, p
[0].im
, p
[1].re
, p
[1].im
,
183 p
[0].re
, p
[0].im
, p
[1].re
, p
[1].im
);
194 BF(p
[0].re
, p
[0].im
, p
[2].re
, p
[2].im
,
195 p
[0].re
, p
[0].im
, p
[2].re
, p
[2].im
);
196 BF(p
[1].re
, p
[1].im
, p
[3].re
, p
[3].im
,
197 p
[1].re
, p
[1].im
, -p
[3].im
, p
[3].re
);
202 BF(p
[0].re
, p
[0].im
, p
[2].re
, p
[2].im
,
203 p
[0].re
, p
[0].im
, p
[2].re
, p
[2].im
);
204 BF(p
[1].re
, p
[1].im
, p
[3].re
, p
[3].im
,
205 p
[1].re
, p
[1].im
, p
[3].im
, -p
[3].re
);
217 for (j
= 0; j
< nblocks
; ++j
) {
218 BF(p
->re
, p
->im
, q
->re
, q
->im
,
219 p
->re
, p
->im
, q
->re
, q
->im
);
223 for(l
= nblocks
; l
< np2
; l
+= nblocks
) {
224 CMUL(tmp_re
, tmp_im
, exptab
[l
].re
, exptab
[l
].im
, q
->re
, q
->im
);
225 BF(p
->re
, p
->im
, q
->re
, q
->im
,
226 p
->re
, p
->im
, tmp_re
, tmp_im
);
234 nblocks
= nblocks
>> 1;
235 nloops
= nloops
<< 1;
236 } while (nblocks
!= 0);
240 * Do the permutation needed BEFORE calling ff_fft_calc()
242 void ff_fft_permute(FFTContext
*s
, FFTComplex
*z
)
246 const uint16_t *revtab
= s
->revtab
;
260 void ff_fft_end(FFTContext
*s
)
262 av_freep(&s
->revtab
);
263 av_freep(&s
->exptab
);
264 av_freep(&s
->exptab1
);