diff options
Diffstat (limited to 'lib/rbcodec/codecs/libwmavoice/fft.c')
-rw-r--r-- | lib/rbcodec/codecs/libwmavoice/fft.c | 298 |
1 files changed, 298 insertions, 0 deletions
diff --git a/lib/rbcodec/codecs/libwmavoice/fft.c b/lib/rbcodec/codecs/libwmavoice/fft.c new file mode 100644 index 0000000000..a030534a26 --- /dev/null +++ b/lib/rbcodec/codecs/libwmavoice/fft.c | |||
@@ -0,0 +1,298 @@ | |||
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 | ||
26 | * FFT/IFFT transforms. | ||
27 | */ | ||
28 | |||
29 | #include <stdlib.h> | ||
30 | #include <string.h> | ||
31 | #include "libavutil/mathematics.h" | ||
32 | #include "fft.h" | ||
33 | |||
34 | /* cos(2*pi*x/n) for 0<=x<=n/4, followed by its reverse */ | ||
35 | #if !CONFIG_HARDCODED_TABLES | ||
36 | COSTABLE(16); | ||
37 | COSTABLE(32); | ||
38 | COSTABLE(64); | ||
39 | COSTABLE(128); | ||
40 | COSTABLE(256); | ||
41 | COSTABLE(512); | ||
42 | COSTABLE(1024); | ||
43 | COSTABLE(2048); | ||
44 | COSTABLE(4096); | ||
45 | COSTABLE(8192); | ||
46 | COSTABLE(16384); | ||
47 | COSTABLE(32768); | ||
48 | COSTABLE(65536); | ||
49 | #endif | ||
50 | COSTABLE_CONST FFTSample * const ff_cos_tabs[] = { | ||
51 | NULL, NULL, NULL, NULL, | ||
52 | ff_cos_16, ff_cos_32, ff_cos_64, ff_cos_128, ff_cos_256, ff_cos_512, ff_cos_1024, | ||
53 | ff_cos_2048, ff_cos_4096, ff_cos_8192, ff_cos_16384, ff_cos_32768, ff_cos_65536, | ||
54 | }; | ||
55 | |||
56 | static int split_radix_permutation(int i, int n, int inverse) | ||
57 | { | ||
58 | int m; | ||
59 | if(n <= 2) return i&1; | ||
60 | m = n >> 1; | ||
61 | if(!(i&m)) return split_radix_permutation(i, m, inverse)*2; | ||
62 | m >>= 1; | ||
63 | if(inverse == !(i&m)) return split_radix_permutation(i, m, inverse)*4 + 1; | ||
64 | else return split_radix_permutation(i, m, inverse)*4 - 1; | ||
65 | } | ||
66 | |||
67 | av_cold void ff_init_ff_cos_tabs(int index) | ||
68 | { | ||
69 | #if !CONFIG_HARDCODED_TABLES | ||
70 | int i; | ||
71 | int m = 1<<index; | ||
72 | double freq = 2*M_PI/m; | ||
73 | FFTSample *tab = ff_cos_tabs[index]; | ||
74 | for(i=0; i<=m/4; i++) | ||
75 | tab[i] = cos(i*freq); | ||
76 | for(i=1; i<m/4; i++) | ||
77 | tab[m/2-i] = tab[i]; | ||
78 | #endif | ||
79 | } | ||
80 | |||
81 | av_cold int ff_fft_init(FFTContext *s, int nbits, int inverse) | ||
82 | { | ||
83 | int i, j, n; | ||
84 | |||
85 | if (nbits < 2 || nbits > 16) | ||
86 | goto fail; | ||
87 | s->nbits = nbits; | ||
88 | n = 1 << nbits; | ||
89 | |||
90 | s->revtab = av_malloc(n * sizeof(uint16_t)); | ||
91 | if (!s->revtab) | ||
92 | goto fail; | ||
93 | s->tmp_buf = av_malloc(n * sizeof(FFTComplex)); | ||
94 | if (!s->tmp_buf) | ||
95 | goto fail; | ||
96 | s->inverse = inverse; | ||
97 | |||
98 | s->fft_permute = ff_fft_permute_c; | ||
99 | s->fft_calc = ff_fft_calc_c; | ||
100 | #if CONFIG_MDCT | ||
101 | s->imdct_calc = ff_imdct_calc_c; | ||
102 | s->imdct_half = ff_imdct_half_c; | ||
103 | s->mdct_calc = ff_mdct_calc_c; | ||
104 | #endif | ||
105 | |||
106 | #if 0 | ||
107 | if (ARCH_ARM) ff_fft_init_arm(s); | ||
108 | if (HAVE_ALTIVEC) ff_fft_init_altivec(s); | ||
109 | if (HAVE_MMX) ff_fft_init_mmx(s); | ||
110 | #endif | ||
111 | |||
112 | for(j=4; j<=nbits; j++) { | ||
113 | ff_init_ff_cos_tabs(j); | ||
114 | } | ||
115 | for(i=0; i<n; i++) | ||
116 | s->revtab[-split_radix_permutation(i, n, s->inverse) & (n-1)] = i; | ||
117 | |||
118 | return 0; | ||
119 | fail: | ||
120 | av_freep(&s->revtab); | ||
121 | av_freep(&s->tmp_buf); | ||
122 | return -1; | ||
123 | } | ||
124 | |||
125 | void ff_fft_permute_c(FFTContext *s, FFTComplex *z) | ||
126 | { | ||
127 | int j, np; | ||
128 | const uint16_t *revtab = s->revtab; | ||
129 | np = 1 << s->nbits; | ||
130 | /* TODO: handle split-radix permute in a more optimal way, probably in-place */ | ||
131 | for(j=0;j<np;j++) s->tmp_buf[revtab[j]] = z[j]; | ||
132 | memcpy(z, s->tmp_buf, np * sizeof(FFTComplex)); | ||
133 | } | ||
134 | |||
135 | av_cold void ff_fft_end(FFTContext *s) | ||
136 | { | ||
137 | av_freep(&s->revtab); | ||
138 | av_freep(&s->tmp_buf); | ||
139 | } | ||
140 | |||
141 | #define sqrthalf (float)M_SQRT1_2 | ||
142 | |||
143 | #define BF(x,y,a,b) {\ | ||
144 | x = a - b;\ | ||
145 | y = a + b;\ | ||
146 | } | ||
147 | |||
148 | #define BUTTERFLIES(a0,a1,a2,a3) {\ | ||
149 | BF(t3, t5, t5, t1);\ | ||
150 | BF(a2.re, a0.re, a0.re, t5);\ | ||
151 | BF(a3.im, a1.im, a1.im, t3);\ | ||
152 | BF(t4, t6, t2, t6);\ | ||
153 | BF(a3.re, a1.re, a1.re, t4);\ | ||
154 | BF(a2.im, a0.im, a0.im, t6);\ | ||
155 | } | ||
156 | |||
157 | // force loading all the inputs before storing any. | ||
158 | // this is slightly slower for small data, but avoids store->load aliasing | ||
159 | // for addresses separated by large powers of 2. | ||
160 | #define BUTTERFLIES_BIG(a0,a1,a2,a3) {\ | ||
161 | FFTSample r0=a0.re, i0=a0.im, r1=a1.re, i1=a1.im;\ | ||
162 | BF(t3, t5, t5, t1);\ | ||
163 | BF(a2.re, a0.re, r0, t5);\ | ||
164 | BF(a3.im, a1.im, i1, t3);\ | ||
165 | BF(t4, t6, t2, t6);\ | ||
166 | BF(a3.re, a1.re, r1, t4);\ | ||
167 | BF(a2.im, a0.im, i0, t6);\ | ||
168 | } | ||
169 | |||
170 | #define TRANSFORM(a0,a1,a2,a3,wre,wim) {\ | ||
171 | t1 = a2.re * wre + a2.im * wim;\ | ||
172 | t2 = a2.im * wre - a2.re * wim;\ | ||
173 | t5 = a3.re * wre - a3.im * wim;\ | ||
174 | t6 = a3.im * wre + a3.re * wim;\ | ||
175 | BUTTERFLIES(a0,a1,a2,a3)\ | ||
176 | } | ||
177 | |||
178 | #define TRANSFORM_ZERO(a0,a1,a2,a3) {\ | ||
179 | t1 = a2.re;\ | ||
180 | t2 = a2.im;\ | ||
181 | t5 = a3.re;\ | ||
182 | t6 = a3.im;\ | ||
183 | BUTTERFLIES(a0,a1,a2,a3)\ | ||
184 | } | ||
185 | |||
186 | /* z[0...8n-1], w[1...2n-1] */ | ||
187 | #define PASS(name)\ | ||
188 | static void name(FFTComplex *z, const FFTSample *wre, unsigned int n)\ | ||
189 | {\ | ||
190 | FFTSample t1, t2, t3, t4, t5, t6;\ | ||
191 | int o1 = 2*n;\ | ||
192 | int o2 = 4*n;\ | ||
193 | int o3 = 6*n;\ | ||
194 | const FFTSample *wim = wre+o1;\ | ||
195 | n--;\ | ||
196 | \ | ||
197 | TRANSFORM_ZERO(z[0],z[o1],z[o2],z[o3]);\ | ||
198 | TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\ | ||
199 | do {\ | ||
200 | z += 2;\ | ||
201 | wre += 2;\ | ||
202 | wim -= 2;\ | ||
203 | TRANSFORM(z[0],z[o1],z[o2],z[o3],wre[0],wim[0]);\ | ||
204 | TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\ | ||
205 | } while(--n);\ | ||
206 | } | ||
207 | |||
208 | PASS(pass) | ||
209 | #undef BUTTERFLIES | ||
210 | #define BUTTERFLIES BUTTERFLIES_BIG | ||
211 | PASS(pass_big) | ||
212 | |||
213 | #define DECL_FFT(n,n2,n4)\ | ||
214 | static void fft##n(FFTComplex *z)\ | ||
215 | {\ | ||
216 | fft##n2(z);\ | ||
217 | fft##n4(z+n4*2);\ | ||
218 | fft##n4(z+n4*3);\ | ||
219 | pass(z,ff_cos_##n,n4/2);\ | ||
220 | } | ||
221 | |||
222 | static void fft4(FFTComplex *z) | ||
223 | { | ||
224 | FFTSample t1, t2, t3, t4, t5, t6, t7, t8; | ||
225 | |||
226 | BF(t3, t1, z[0].re, z[1].re); | ||
227 | BF(t8, t6, z[3].re, z[2].re); | ||
228 | BF(z[2].re, z[0].re, t1, t6); | ||
229 | BF(t4, t2, z[0].im, z[1].im); | ||
230 | BF(t7, t5, z[2].im, z[3].im); | ||
231 | BF(z[3].im, z[1].im, t4, t8); | ||
232 | BF(z[3].re, z[1].re, t3, t7); | ||
233 | BF(z[2].im, z[0].im, t2, t5); | ||
234 | } | ||
235 | |||
236 | static void fft8(FFTComplex *z) | ||
237 | { | ||
238 | FFTSample t1, t2, t3, t4, t5, t6, t7, t8; | ||
239 | |||
240 | fft4(z); | ||
241 | |||
242 | BF(t1, z[5].re, z[4].re, -z[5].re); | ||
243 | BF(t2, z[5].im, z[4].im, -z[5].im); | ||
244 | BF(t3, z[7].re, z[6].re, -z[7].re); | ||
245 | BF(t4, z[7].im, z[6].im, -z[7].im); | ||
246 | BF(t8, t1, t3, t1); | ||
247 | BF(t7, t2, t2, t4); | ||
248 | BF(z[4].re, z[0].re, z[0].re, t1); | ||
249 | BF(z[4].im, z[0].im, z[0].im, t2); | ||
250 | BF(z[6].re, z[2].re, z[2].re, t7); | ||
251 | BF(z[6].im, z[2].im, z[2].im, t8); | ||
252 | |||
253 | TRANSFORM(z[1],z[3],z[5],z[7],sqrthalf,sqrthalf); | ||
254 | } | ||
255 | |||
256 | #if !CONFIG_SMALL | ||
257 | static void fft16(FFTComplex *z) | ||
258 | { | ||
259 | FFTSample t1, t2, t3, t4, t5, t6; | ||
260 | |||
261 | fft8(z); | ||
262 | fft4(z+8); | ||
263 | fft4(z+12); | ||
264 | |||
265 | TRANSFORM_ZERO(z[0],z[4],z[8],z[12]); | ||
266 | TRANSFORM(z[2],z[6],z[10],z[14],sqrthalf,sqrthalf); | ||
267 | TRANSFORM(z[1],z[5],z[9],z[13],ff_cos_16[1],ff_cos_16[3]); | ||
268 | TRANSFORM(z[3],z[7],z[11],z[15],ff_cos_16[3],ff_cos_16[1]); | ||
269 | } | ||
270 | #else | ||
271 | DECL_FFT(16,8,4) | ||
272 | #endif | ||
273 | DECL_FFT(32,16,8) | ||
274 | DECL_FFT(64,32,16) | ||
275 | DECL_FFT(128,64,32) | ||
276 | DECL_FFT(256,128,64) | ||
277 | DECL_FFT(512,256,128) | ||
278 | #if !CONFIG_SMALL | ||
279 | #define pass pass_big | ||
280 | #endif | ||
281 | DECL_FFT(1024,512,256) | ||
282 | DECL_FFT(2048,1024,512) | ||
283 | DECL_FFT(4096,2048,1024) | ||
284 | DECL_FFT(8192,4096,2048) | ||
285 | DECL_FFT(16384,8192,4096) | ||
286 | DECL_FFT(32768,16384,8192) | ||
287 | DECL_FFT(65536,32768,16384) | ||
288 | |||
289 | static void (* const fft_dispatch[])(FFTComplex*) = { | ||
290 | fft4, fft8, fft16, fft32, fft64, fft128, fft256, fft512, fft1024, | ||
291 | fft2048, fft4096, fft8192, fft16384, fft32768, fft65536, | ||
292 | }; | ||
293 | |||
294 | void ff_fft_calc_c(FFTContext *s, FFTComplex *z) | ||
295 | { | ||
296 | fft_dispatch[s->nbits-2](z); | ||
297 | } | ||
298 | |||