Mercurial > libavcodec.hg
comparison fft.c @ 781:6f5e87957bcb libavcodec
new generic FFT/MDCT code for audio codecs
author | bellard |
---|---|
date | Mon, 28 Oct 2002 00:34:08 +0000 |
parents | |
children | 6d8b45263d6c |
comparison
equal
deleted
inserted
replaced
780:a48bb8bc63dd | 781:6f5e87957bcb |
---|---|
1 /* | |
2 * FFT/IFFT transforms | |
3 * Copyright (c) 2002 Fabrice Bellard. | |
4 * | |
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. | |
9 * | |
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. | |
14 * | |
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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA | |
18 */ | |
19 #include "dsputil.h" | |
20 | |
21 /** | |
22 * The size of the FFT is 2^nbits. If inverse is TRUE, inverse FFT is | |
23 * done | |
24 */ | |
25 int fft_init(FFTContext *s, int nbits, int inverse) | |
26 { | |
27 int i, j, m, n; | |
28 float alpha, c1, s1, s2; | |
29 | |
30 s->nbits = nbits; | |
31 n = 1 << nbits; | |
32 | |
33 s->exptab = av_malloc((n / 2) * sizeof(FFTComplex)); | |
34 if (!s->exptab) | |
35 goto fail; | |
36 s->revtab = av_malloc(n * sizeof(uint16_t)); | |
37 if (!s->revtab) | |
38 goto fail; | |
39 s->inverse = inverse; | |
40 | |
41 s2 = inverse ? 1.0 : -1.0; | |
42 | |
43 for(i=0;i<(n/2);i++) { | |
44 alpha = 2 * M_PI * (float)i / (float)n; | |
45 c1 = cos(alpha); | |
46 s1 = sin(alpha) * s2; | |
47 s->exptab[i].re = c1; | |
48 s->exptab[i].im = s1; | |
49 } | |
50 s->fft_calc = fft_calc_c; | |
51 s->exptab1 = NULL; | |
52 | |
53 /* compute constant table for HAVE_SSE version */ | |
54 #if defined(HAVE_MMX) && 0 | |
55 if (mm_flags & MM_SSE) { | |
56 int np, nblocks, np2, l; | |
57 FFTComplex *q; | |
58 | |
59 np = 1 << nbits; | |
60 nblocks = np >> 3; | |
61 np2 = np >> 1; | |
62 s->exptab1 = av_malloc(np * 2 * sizeof(FFTComplex)); | |
63 if (!s->exptab1) | |
64 goto fail; | |
65 q = s->exptab1; | |
66 do { | |
67 for(l = 0; l < np2; l += 2 * nblocks) { | |
68 *q++ = s->exptab[l]; | |
69 *q++ = s->exptab[l + nblocks]; | |
70 | |
71 q->re = -s->exptab[l].im; | |
72 q->im = s->exptab[l].re; | |
73 q++; | |
74 q->re = -s->exptab[l + nblocks].im; | |
75 q->im = s->exptab[l + nblocks].re; | |
76 q++; | |
77 } | |
78 nblocks = nblocks >> 1; | |
79 } while (nblocks != 0); | |
80 av_freep(&s->exptab); | |
81 } | |
82 #endif | |
83 | |
84 /* compute bit reverse table */ | |
85 | |
86 for(i=0;i<n;i++) { | |
87 m=0; | |
88 for(j=0;j<nbits;j++) { | |
89 m |= ((i >> j) & 1) << (nbits-j-1); | |
90 } | |
91 s->revtab[i]=m; | |
92 } | |
93 return 0; | |
94 fail: | |
95 av_freep(&s->revtab); | |
96 av_freep(&s->exptab); | |
97 av_freep(&s->exptab1); | |
98 return -1; | |
99 } | |
100 | |
101 /* butter fly op */ | |
102 #define BF(pre, pim, qre, qim, pre1, pim1, qre1, qim1) \ | |
103 {\ | |
104 FFTSample ax, ay, bx, by;\ | |
105 bx=pre1;\ | |
106 by=pim1;\ | |
107 ax=qre1;\ | |
108 ay=qim1;\ | |
109 pre = (bx + ax);\ | |
110 pim = (by + ay);\ | |
111 qre = (bx - ax);\ | |
112 qim = (by - ay);\ | |
113 } | |
114 | |
115 #define MUL16(a,b) ((a) * (b)) | |
116 | |
117 #define CMUL(pre, pim, are, aim, bre, bim) \ | |
118 {\ | |
119 pre = (MUL16(are, bre) - MUL16(aim, bim));\ | |
120 pim = (MUL16(are, bim) + MUL16(bre, aim));\ | |
121 } | |
122 | |
123 /** | |
124 * Do a complex FFT with the parameters defined in fft_init(). The | |
125 * input data must be permuted before with s->revtab table. No | |
126 * 1.0/sqrt(n) normalization is done. | |
127 */ | |
128 void fft_calc_c(FFTContext *s, FFTComplex *z) | |
129 { | |
130 int ln = s->nbits; | |
131 int j, np, np2; | |
132 int nblocks, nloops; | |
133 register FFTComplex *p, *q; | |
134 FFTComplex *exptab = s->exptab; | |
135 int l; | |
136 FFTSample tmp_re, tmp_im; | |
137 | |
138 np = 1 << ln; | |
139 | |
140 /* pass 0 */ | |
141 | |
142 p=&z[0]; | |
143 j=(np >> 1); | |
144 do { | |
145 BF(p[0].re, p[0].im, p[1].re, p[1].im, | |
146 p[0].re, p[0].im, p[1].re, p[1].im); | |
147 p+=2; | |
148 } while (--j != 0); | |
149 | |
150 /* pass 1 */ | |
151 | |
152 | |
153 p=&z[0]; | |
154 j=np >> 2; | |
155 if (s->inverse) { | |
156 do { | |
157 BF(p[0].re, p[0].im, p[2].re, p[2].im, | |
158 p[0].re, p[0].im, p[2].re, p[2].im); | |
159 BF(p[1].re, p[1].im, p[3].re, p[3].im, | |
160 p[1].re, p[1].im, -p[3].im, p[3].re); | |
161 p+=4; | |
162 } while (--j != 0); | |
163 } else { | |
164 do { | |
165 BF(p[0].re, p[0].im, p[2].re, p[2].im, | |
166 p[0].re, p[0].im, p[2].re, p[2].im); | |
167 BF(p[1].re, p[1].im, p[3].re, p[3].im, | |
168 p[1].re, p[1].im, p[3].im, -p[3].re); | |
169 p+=4; | |
170 } while (--j != 0); | |
171 } | |
172 /* pass 2 .. ln-1 */ | |
173 | |
174 nblocks = np >> 3; | |
175 nloops = 1 << 2; | |
176 np2 = np >> 1; | |
177 do { | |
178 p = z; | |
179 q = z + nloops; | |
180 for (j = 0; j < nblocks; ++j) { | |
181 BF(p->re, p->im, q->re, q->im, | |
182 p->re, p->im, q->re, q->im); | |
183 | |
184 p++; | |
185 q++; | |
186 for(l = nblocks; l < np2; l += nblocks) { | |
187 CMUL(tmp_re, tmp_im, exptab[l].re, exptab[l].im, q->re, q->im); | |
188 BF(p->re, p->im, q->re, q->im, | |
189 p->re, p->im, tmp_re, tmp_im); | |
190 p++; | |
191 q++; | |
192 } | |
193 | |
194 p += nloops; | |
195 q += nloops; | |
196 } | |
197 nblocks = nblocks >> 1; | |
198 nloops = nloops << 1; | |
199 } while (nblocks != 0); | |
200 } | |
201 | |
202 /** | |
203 * Do the permutation needed BEFORE calling fft_calc() | |
204 */ | |
205 void fft_permute(FFTContext *s, FFTComplex *z) | |
206 { | |
207 int j, k, np; | |
208 FFTComplex tmp; | |
209 const uint16_t *revtab = s->revtab; | |
210 | |
211 /* reverse */ | |
212 np = 1 << s->nbits; | |
213 for(j=0;j<np;j++) { | |
214 k = revtab[j]; | |
215 if (k < j) { | |
216 tmp = z[k]; | |
217 z[k] = z[j]; | |
218 z[j] = tmp; | |
219 } | |
220 } | |
221 } | |
222 | |
223 void fft_end(FFTContext *s) | |
224 { | |
225 av_freep(&s->revtab); | |
226 av_freep(&s->exptab); | |
227 av_freep(&s->exptab1); | |
228 } | |
229 |