Mercurial > libavcodec.hg
annotate fft.c @ 4166:eced83504436 libavcodec
mp3 header (de)compression bitstream filter
this will make mp3 frames 4 bytes smaller, it will not give you binary identical mp3 files, but it will give you mp3 files which decode to binary identical output
this will only work in containers providing at least packet size, sample_rate and number of channels
bugreports about mp3 files for which this fails are welcome
and this is experimental (dont expect compatibility and dont even expect to be able to decompress what you compressed, hell dont even expect this to work without editing the source a little)
author | michael |
---|---|
date | Fri, 10 Nov 2006 01:41:53 +0000 |
parents | c8c591fe26f8 |
children | 311c95dddb16 |
rev | line source |
---|---|
781 | 1 /* |
2 * FFT/IFFT transforms | |
3 * Copyright (c) 2002 Fabrice Bellard. | |
4 * | |
3947
c8c591fe26f8
Change license headers to say 'FFmpeg' instead of 'this program/this library'
diego
parents:
3746
diff
changeset
|
5 * This file is part of FFmpeg. |
c8c591fe26f8
Change license headers to say 'FFmpeg' instead of 'this program/this library'
diego
parents:
3746
diff
changeset
|
6 * |
c8c591fe26f8
Change license headers to say 'FFmpeg' instead of 'this program/this library'
diego
parents:
3746
diff
changeset
|
7 * FFmpeg is free software; you can redistribute it and/or |
781 | 8 * modify it under the terms of the GNU Lesser General Public |
9 * 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
|
10 * version 2.1 of the License, or (at your option) any later version. |
781 | 11 * |
3947
c8c591fe26f8
Change license headers to say 'FFmpeg' instead of 'this program/this library'
diego
parents:
3746
diff
changeset
|
12 * FFmpeg is distributed in the hope that it will be useful, |
781 | 13 * but WITHOUT ANY WARRANTY; without even the implied warranty of |
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
15 * Lesser General Public License for more details. | |
16 * | |
17 * 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
|
18 * 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
|
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
781 | 20 */ |
1106 | 21 |
22 /** | |
23 * @file fft.c | |
24 * FFT/IFFT transforms. | |
25 */ | |
26 | |
781 | 27 #include "dsputil.h" |
28 | |
29 /** | |
30 * The size of the FFT is 2^nbits. If inverse is TRUE, inverse FFT is | |
2967 | 31 * done |
781 | 32 */ |
1879
dd63cb7e5080
fft_*() renamed into ff_fft_*() patch by (Gildas Bazin <gbazin at altern dot org>)
michael
parents:
1106
diff
changeset
|
33 int ff_fft_init(FFTContext *s, int nbits, int inverse) |
781 | 34 { |
35 int i, j, m, n; | |
36 float alpha, c1, s1, s2; | |
2967 | 37 |
781 | 38 s->nbits = nbits; |
39 n = 1 << nbits; | |
40 | |
41 s->exptab = av_malloc((n / 2) * sizeof(FFTComplex)); | |
42 if (!s->exptab) | |
43 goto fail; | |
44 s->revtab = av_malloc(n * sizeof(uint16_t)); | |
45 if (!s->revtab) | |
46 goto fail; | |
47 s->inverse = inverse; | |
48 | |
49 s2 = inverse ? 1.0 : -1.0; | |
2967 | 50 |
781 | 51 for(i=0;i<(n/2);i++) { |
52 alpha = 2 * M_PI * (float)i / (float)n; | |
53 c1 = cos(alpha); | |
54 s1 = sin(alpha) * s2; | |
55 s->exptab[i].re = c1; | |
56 s->exptab[i].im = s1; | |
57 } | |
1879
dd63cb7e5080
fft_*() renamed into ff_fft_*() patch by (Gildas Bazin <gbazin at altern dot org>)
michael
parents:
1106
diff
changeset
|
58 s->fft_calc = ff_fft_calc_c; |
3555 | 59 s->imdct_calc = ff_imdct_calc; |
781 | 60 s->exptab1 = NULL; |
61 | |
62 /* compute constant table for HAVE_SSE version */ | |
3590
a3d97c60ea07
ff_fft_calc_3dn/3dn2/sse: convert intrinsics to inline asm.
lorenm
parents:
3567
diff
changeset
|
63 #if defined(HAVE_MMX) \ |
3567
1f8730f62765
r5954 broke fft on cpus with 3dnow but without mm3dnow.h
lorenm
parents:
3555
diff
changeset
|
64 || (defined(HAVE_ALTIVEC) && !defined(ALTIVEC_USE_REFERENCE_C_CODE)) |
975
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
65 { |
3567
1f8730f62765
r5954 broke fft on cpus with 3dnow but without mm3dnow.h
lorenm
parents:
3555
diff
changeset
|
66 int has_vectors = mm_support(); |
781 | 67 |
3567
1f8730f62765
r5954 broke fft on cpus with 3dnow but without mm3dnow.h
lorenm
parents:
3555
diff
changeset
|
68 if (has_vectors) { |
975
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
69 #if defined(HAVE_MMX) |
3746 | 70 if (has_vectors & MM_3DNOWEXT) { |
71 /* 3DNowEx for K7/K8 */ | |
3567
1f8730f62765
r5954 broke fft on cpus with 3dnow but without mm3dnow.h
lorenm
parents:
3555
diff
changeset
|
72 s->imdct_calc = ff_imdct_calc_3dn2; |
1f8730f62765
r5954 broke fft on cpus with 3dnow but without mm3dnow.h
lorenm
parents:
3555
diff
changeset
|
73 s->fft_calc = ff_fft_calc_3dn2; |
3746 | 74 } else if (has_vectors & MM_3DNOW) { |
3567
1f8730f62765
r5954 broke fft on cpus with 3dnow but without mm3dnow.h
lorenm
parents:
3555
diff
changeset
|
75 /* 3DNow! for K6-2/3 */ |
1f8730f62765
r5954 broke fft on cpus with 3dnow but without mm3dnow.h
lorenm
parents:
3555
diff
changeset
|
76 s->fft_calc = ff_fft_calc_3dn; |
3746 | 77 } else if (has_vectors & MM_SSE) { |
3591 | 78 /* SSE for P3/P4 */ |
3746 | 79 s->imdct_calc = ff_imdct_calc_sse; |
3567
1f8730f62765
r5954 broke fft on cpus with 3dnow but without mm3dnow.h
lorenm
parents:
3555
diff
changeset
|
80 s->fft_calc = ff_fft_calc_sse; |
3746 | 81 } |
3567
1f8730f62765
r5954 broke fft on cpus with 3dnow but without mm3dnow.h
lorenm
parents:
3555
diff
changeset
|
82 #else /* HAVE_MMX */ |
1f8730f62765
r5954 broke fft on cpus with 3dnow but without mm3dnow.h
lorenm
parents:
3555
diff
changeset
|
83 if (has_vectors & MM_ALTIVEC) |
1f8730f62765
r5954 broke fft on cpus with 3dnow but without mm3dnow.h
lorenm
parents:
3555
diff
changeset
|
84 s->fft_calc = ff_fft_calc_altivec; |
1f8730f62765
r5954 broke fft on cpus with 3dnow but without mm3dnow.h
lorenm
parents:
3555
diff
changeset
|
85 #endif |
1f8730f62765
r5954 broke fft on cpus with 3dnow but without mm3dnow.h
lorenm
parents:
3555
diff
changeset
|
86 } |
1f8730f62765
r5954 broke fft on cpus with 3dnow but without mm3dnow.h
lorenm
parents:
3555
diff
changeset
|
87 if (s->fft_calc != ff_fft_calc_c) { |
975
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
88 int np, nblocks, np2, l; |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
89 FFTComplex *q; |
2967 | 90 |
975
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
91 np = 1 << nbits; |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
92 nblocks = np >> 3; |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
93 np2 = np >> 1; |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
94 s->exptab1 = av_malloc(np * 2 * sizeof(FFTComplex)); |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
95 if (!s->exptab1) |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
96 goto fail; |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
97 q = s->exptab1; |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
98 do { |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
99 for(l = 0; l < np2; l += 2 * nblocks) { |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
100 *q++ = s->exptab[l]; |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
101 *q++ = s->exptab[l + nblocks]; |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
102 |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
103 q->re = -s->exptab[l].im; |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
104 q->im = s->exptab[l].re; |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
105 q++; |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
106 q->re = -s->exptab[l + nblocks].im; |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
107 q->im = s->exptab[l + nblocks].re; |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
108 q++; |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
109 } |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
110 nblocks = nblocks >> 1; |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
111 } while (nblocks != 0); |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
112 av_freep(&s->exptab); |
e05d525505c5
fft altivec by Romain Dolbeau - simplified patch, test it on PPC with fft-test and wma decoding
bellard
parents:
971
diff
changeset
|
113 } |
781 | 114 } |
115 #endif | |
116 | |
117 /* compute bit reverse table */ | |
118 | |
119 for(i=0;i<n;i++) { | |
120 m=0; | |
121 for(j=0;j<nbits;j++) { | |
122 m |= ((i >> j) & 1) << (nbits-j-1); | |
123 } | |
124 s->revtab[i]=m; | |
125 } | |
126 return 0; | |
127 fail: | |
128 av_freep(&s->revtab); | |
129 av_freep(&s->exptab); | |
130 av_freep(&s->exptab1); | |
131 return -1; | |
132 } | |
133 | |
134 /* butter fly op */ | |
135 #define BF(pre, pim, qre, qim, pre1, pim1, qre1, qim1) \ | |
136 {\ | |
137 FFTSample ax, ay, bx, by;\ | |
138 bx=pre1;\ | |
139 by=pim1;\ | |
140 ax=qre1;\ | |
141 ay=qim1;\ | |
142 pre = (bx + ax);\ | |
143 pim = (by + ay);\ | |
144 qre = (bx - ax);\ | |
145 qim = (by - ay);\ | |
146 } | |
147 | |
148 #define MUL16(a,b) ((a) * (b)) | |
149 | |
150 #define CMUL(pre, pim, are, aim, bre, bim) \ | |
151 {\ | |
152 pre = (MUL16(are, bre) - MUL16(aim, bim));\ | |
153 pim = (MUL16(are, bim) + MUL16(bre, aim));\ | |
154 } | |
155 | |
156 /** | |
1879
dd63cb7e5080
fft_*() renamed into ff_fft_*() patch by (Gildas Bazin <gbazin at altern dot org>)
michael
parents:
1106
diff
changeset
|
157 * Do a complex FFT with the parameters defined in ff_fft_init(). The |
781 | 158 * input data must be permuted before with s->revtab table. No |
2967 | 159 * 1.0/sqrt(n) normalization is done. |
781 | 160 */ |
1879
dd63cb7e5080
fft_*() renamed into ff_fft_*() patch by (Gildas Bazin <gbazin at altern dot org>)
michael
parents:
1106
diff
changeset
|
161 void ff_fft_calc_c(FFTContext *s, FFTComplex *z) |
781 | 162 { |
163 int ln = s->nbits; | |
2979 | 164 int j, np, np2; |
165 int nblocks, nloops; | |
781 | 166 register FFTComplex *p, *q; |
167 FFTComplex *exptab = s->exptab; | |
168 int l; | |
169 FFTSample tmp_re, tmp_im; | |
170 | |
171 np = 1 << ln; | |
172 | |
173 /* pass 0 */ | |
174 | |
175 p=&z[0]; | |
176 j=(np >> 1); | |
177 do { | |
2967 | 178 BF(p[0].re, p[0].im, p[1].re, p[1].im, |
781 | 179 p[0].re, p[0].im, p[1].re, p[1].im); |
180 p+=2; | |
181 } while (--j != 0); | |
182 | |
183 /* pass 1 */ | |
184 | |
2967 | 185 |
781 | 186 p=&z[0]; |
187 j=np >> 2; | |
188 if (s->inverse) { | |
189 do { | |
2967 | 190 BF(p[0].re, p[0].im, p[2].re, p[2].im, |
781 | 191 p[0].re, p[0].im, p[2].re, p[2].im); |
2967 | 192 BF(p[1].re, p[1].im, p[3].re, p[3].im, |
781 | 193 p[1].re, p[1].im, -p[3].im, p[3].re); |
194 p+=4; | |
195 } while (--j != 0); | |
196 } else { | |
197 do { | |
2967 | 198 BF(p[0].re, p[0].im, p[2].re, p[2].im, |
781 | 199 p[0].re, p[0].im, p[2].re, p[2].im); |
2967 | 200 BF(p[1].re, p[1].im, p[3].re, p[3].im, |
781 | 201 p[1].re, p[1].im, p[3].im, -p[3].re); |
202 p+=4; | |
203 } while (--j != 0); | |
204 } | |
205 /* pass 2 .. ln-1 */ | |
206 | |
207 nblocks = np >> 3; | |
208 nloops = 1 << 2; | |
209 np2 = np >> 1; | |
210 do { | |
211 p = z; | |
212 q = z + nloops; | |
213 for (j = 0; j < nblocks; ++j) { | |
214 BF(p->re, p->im, q->re, q->im, | |
215 p->re, p->im, q->re, q->im); | |
2967 | 216 |
781 | 217 p++; |
218 q++; | |
219 for(l = nblocks; l < np2; l += nblocks) { | |
220 CMUL(tmp_re, tmp_im, exptab[l].re, exptab[l].im, q->re, q->im); | |
221 BF(p->re, p->im, q->re, q->im, | |
222 p->re, p->im, tmp_re, tmp_im); | |
223 p++; | |
224 q++; | |
225 } | |
226 | |
227 p += nloops; | |
228 q += nloops; | |
229 } | |
230 nblocks = nblocks >> 1; | |
231 nloops = nloops << 1; | |
232 } while (nblocks != 0); | |
233 } | |
234 | |
235 /** | |
1879
dd63cb7e5080
fft_*() renamed into ff_fft_*() patch by (Gildas Bazin <gbazin at altern dot org>)
michael
parents:
1106
diff
changeset
|
236 * Do the permutation needed BEFORE calling ff_fft_calc() |
781 | 237 */ |
1879
dd63cb7e5080
fft_*() renamed into ff_fft_*() patch by (Gildas Bazin <gbazin at altern dot org>)
michael
parents:
1106
diff
changeset
|
238 void ff_fft_permute(FFTContext *s, FFTComplex *z) |
781 | 239 { |
240 int j, k, np; | |
241 FFTComplex tmp; | |
242 const uint16_t *revtab = s->revtab; | |
2967 | 243 |
781 | 244 /* reverse */ |
245 np = 1 << s->nbits; | |
246 for(j=0;j<np;j++) { | |
247 k = revtab[j]; | |
248 if (k < j) { | |
249 tmp = z[k]; | |
250 z[k] = z[j]; | |
251 z[j] = tmp; | |
252 } | |
253 } | |
254 } | |
255 | |
1879
dd63cb7e5080
fft_*() renamed into ff_fft_*() patch by (Gildas Bazin <gbazin at altern dot org>)
michael
parents:
1106
diff
changeset
|
256 void ff_fft_end(FFTContext *s) |
781 | 257 { |
258 av_freep(&s->revtab); | |
259 av_freep(&s->exptab); | |
260 av_freep(&s->exptab1); | |
261 } | |
262 |