98bf6f583cc0f0c9ee127312bc4455f54b218759
[libav.git] / libavcodec / fft-test.c
1 /*
2 * (c) 2002 Fabrice Bellard
3 *
4 * This file is part of Libav.
5 *
6 * Libav is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU Lesser General Public
8 * License as published by the Free Software Foundation; either
9 * version 2.1 of the License, or (at your option) any later version.
10 *
11 * Libav is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with Libav; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19 */
20
21 /**
22 * @file
23 * FFT and MDCT tests.
24 */
25
26 #include "libavutil/cpu.h"
27 #include "libavutil/mathematics.h"
28 #include "libavutil/lfg.h"
29 #include "libavutil/log.h"
30 #include "libavutil/time.h"
31 #include "fft.h"
32 #include "dct.h"
33 #include "rdft.h"
34 #include <math.h>
35 #if HAVE_UNISTD_H
36 #include <unistd.h>
37 #endif
38 #include <stdio.h>
39 #include <stdlib.h>
40 #include <string.h>
41
42 /* reference fft */
43
44 #define MUL16(a,b) ((a) * (b))
45
46 #define CMAC(pre, pim, are, aim, bre, bim) \
47 {\
48 pre += (MUL16(are, bre) - MUL16(aim, bim));\
49 pim += (MUL16(are, bim) + MUL16(bre, aim));\
50 }
51
52 #if FFT_FLOAT
53 # define RANGE 1.0
54 # define REF_SCALE(x, bits) (x)
55 # define FMT "%10.6f"
56 #else
57 # define RANGE 16384
58 # define REF_SCALE(x, bits) ((x) / (1<<(bits)))
59 # define FMT "%6d"
60 #endif
61
62 static struct {
63 float re, im;
64 } *exptab;
65
66 static int fft_ref_init(int nbits, int inverse)
67 {
68 int i, n = 1 << nbits;
69
70 exptab = av_malloc((n / 2) * sizeof(*exptab));
71 if (!exptab)
72 return AVERROR(ENOMEM);
73
74 for (i = 0; i < (n/2); i++) {
75 double alpha = 2 * M_PI * (float)i / (float)n;
76 double c1 = cos(alpha), s1 = sin(alpha);
77 if (!inverse)
78 s1 = -s1;
79 exptab[i].re = c1;
80 exptab[i].im = s1;
81 }
82 return 0;
83 }
84
85 static void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits)
86 {
87 int i, j;
88 int n = 1 << nbits;
89 int n2 = n >> 1;
90
91 for (i = 0; i < n; i++) {
92 double tmp_re = 0, tmp_im = 0;
93 FFTComplex *q = tab;
94 for (j = 0; j < n; j++) {
95 double s, c;
96 int k = (i * j) & (n - 1);
97 if (k >= n2) {
98 c = -exptab[k - n2].re;
99 s = -exptab[k - n2].im;
100 } else {
101 c = exptab[k].re;
102 s = exptab[k].im;
103 }
104 CMAC(tmp_re, tmp_im, c, s, q->re, q->im);
105 q++;
106 }
107 tabr[i].re = REF_SCALE(tmp_re, nbits);
108 tabr[i].im = REF_SCALE(tmp_im, nbits);
109 }
110 }
111
112 #if CONFIG_MDCT
113 static void imdct_ref(FFTSample *out, FFTSample *in, int nbits)
114 {
115 int i, k, n = 1 << nbits;
116
117 for (i = 0; i < n; i++) {
118 double sum = 0;
119 for (k = 0; k < n/2; k++) {
120 int a = (2 * i + 1 + (n / 2)) * (2 * k + 1);
121 double f = cos(M_PI * a / (double)(2 * n));
122 sum += f * in[k];
123 }
124 out[i] = REF_SCALE(-sum, nbits - 2);
125 }
126 }
127
128 /* NOTE: no normalisation by 1 / N is done */
129 static void mdct_ref(FFTSample *output, FFTSample *input, int nbits)
130 {
131 int i, k, n = 1 << nbits;
132
133 /* do it by hand */
134 for (k = 0; k < n/2; k++) {
135 double s = 0;
136 for (i = 0; i < n; i++) {
137 double a = (2 * M_PI * (2 * i + 1 + n / 2) * (2 * k + 1) / (4 * n));
138 s += input[i] * cos(a);
139 }
140 output[k] = REF_SCALE(s, nbits - 1);
141 }
142 }
143 #endif /* CONFIG_MDCT */
144
145 #if FFT_FLOAT
146 #if CONFIG_DCT
147 static void idct_ref(float *output, float *input, int nbits)
148 {
149 int i, k, n = 1 << nbits;
150
151 /* do it by hand */
152 for (i = 0; i < n; i++) {
153 double s = 0.5 * input[0];
154 for (k = 1; k < n; k++) {
155 double a = M_PI * k * (i + 0.5) / n;
156 s += input[k] * cos(a);
157 }
158 output[i] = 2 * s / n;
159 }
160 }
161 static void dct_ref(float *output, float *input, int nbits)
162 {
163 int i, k, n = 1 << nbits;
164
165 /* do it by hand */
166 for (k = 0; k < n; k++) {
167 double s = 0;
168 for (i = 0; i < n; i++) {
169 double a = M_PI * k * (i + 0.5) / n;
170 s += input[i] * cos(a);
171 }
172 output[k] = s;
173 }
174 }
175 #endif /* CONFIG_DCT */
176 #endif /* FFT_FLOAT */
177
178
179 static FFTSample frandom(AVLFG *prng)
180 {
181 return (int16_t)av_lfg_get(prng) / 32768.0 * RANGE;
182 }
183
184 static int check_diff(FFTSample *tab1, FFTSample *tab2, int n, double scale)
185 {
186 int i, err = 0;
187 double error = 0, max = 0;
188
189 for (i = 0; i < n; i++) {
190 double e = fabsf(tab1[i] - (tab2[i] / scale)) / RANGE;
191 if (e >= 1e-3) {
192 av_log(NULL, AV_LOG_ERROR, "ERROR %5d: "FMT" "FMT"\n",
193 i, tab1[i], tab2[i]);
194 err = 1;
195 }
196 error+= e*e;
197 if(e>max) max= e;
198 }
199 av_log(NULL, AV_LOG_INFO, "max:%f e:%g\n", max, sqrt(error)/n);
200 return err;
201 }
202
203
204 static void help(void)
205 {
206 av_log(NULL, AV_LOG_INFO,"usage: fft-test [-h] [-s] [-i] [-n b]\n"
207 "-h print this help\n"
208 "-s speed test\n"
209 "-m (I)MDCT test\n"
210 "-d (I)DCT test\n"
211 "-r (I)RDFT test\n"
212 "-i inverse transform test\n"
213 "-n b set the transform size to 2^b\n"
214 "-f x set scale factor for output data of (I)MDCT to x\n"
215 );
216 }
217
218 enum tf_transform {
219 TRANSFORM_FFT,
220 TRANSFORM_MDCT,
221 TRANSFORM_RDFT,
222 TRANSFORM_DCT,
223 };
224
225 #if !HAVE_GETOPT
226 #include "compat/getopt.c"
227 #endif
228
229 int main(int argc, char **argv)
230 {
231 FFTComplex *tab, *tab1, *tab_ref;
232 FFTSample *tab2;
233 enum tf_transform transform = TRANSFORM_FFT;
234 FFTContext m, s;
235 #if FFT_FLOAT
236 RDFTContext r;
237 DCTContext d;
238 #endif /* FFT_FLOAT */
239 int it, i, err = 1;
240 int do_speed = 0, do_inverse = 0;
241 int fft_nbits = 9, fft_size;
242 double scale = 1.0;
243 AVLFG prng;
244
245 av_lfg_init(&prng, 1);
246
247 for(;;) {
248 int c = getopt(argc, argv, "hsimrdn:f:c:");
249 if (c == -1)
250 break;
251 switch(c) {
252 case 'h':
253 help();
254 return 1;
255 case 's':
256 do_speed = 1;
257 break;
258 case 'i':
259 do_inverse = 1;
260 break;
261 case 'm':
262 transform = TRANSFORM_MDCT;
263 break;
264 case 'r':
265 transform = TRANSFORM_RDFT;
266 break;
267 case 'd':
268 transform = TRANSFORM_DCT;
269 break;
270 case 'n':
271 fft_nbits = atoi(optarg);
272 break;
273 case 'f':
274 scale = atof(optarg);
275 break;
276 case 'c':
277 {
278 int cpuflags = av_parse_cpu_flags(optarg);
279 if (cpuflags < 0)
280 return 1;
281 av_set_cpu_flags_mask(cpuflags);
282 break;
283 }
284 }
285 }
286
287 fft_size = 1 << fft_nbits;
288 tab = av_malloc(fft_size * sizeof(FFTComplex));
289 tab1 = av_malloc(fft_size * sizeof(FFTComplex));
290 tab_ref = av_malloc(fft_size * sizeof(FFTComplex));
291 tab2 = av_malloc(fft_size * sizeof(FFTSample));
292
293 if (!(tab && tab1 && tab_ref && tab2))
294 goto cleanup;
295
296 switch (transform) {
297 #if CONFIG_MDCT
298 case TRANSFORM_MDCT:
299 av_log(NULL, AV_LOG_INFO,"Scale factor is set to %f\n", scale);
300 if (do_inverse)
301 av_log(NULL, AV_LOG_INFO,"IMDCT");
302 else
303 av_log(NULL, AV_LOG_INFO,"MDCT");
304 ff_mdct_init(&m, fft_nbits, do_inverse, scale);
305 break;
306 #endif /* CONFIG_MDCT */
307 case TRANSFORM_FFT:
308 if (do_inverse)
309 av_log(NULL, AV_LOG_INFO,"IFFT");
310 else
311 av_log(NULL, AV_LOG_INFO,"FFT");
312 ff_fft_init(&s, fft_nbits, do_inverse);
313 if (err = fft_ref_init(fft_nbits, do_inverse) < 0)
314 goto cleanup;
315 break;
316 #if FFT_FLOAT
317 #if CONFIG_RDFT
318 case TRANSFORM_RDFT:
319 if (do_inverse)
320 av_log(NULL, AV_LOG_INFO,"IDFT_C2R");
321 else
322 av_log(NULL, AV_LOG_INFO,"DFT_R2C");
323 ff_rdft_init(&r, fft_nbits, do_inverse ? IDFT_C2R : DFT_R2C);
324 if (err = fft_ref_init(fft_nbits, do_inverse) < 0)
325 goto cleanup;
326 break;
327 #endif /* CONFIG_RDFT */
328 #if CONFIG_DCT
329 case TRANSFORM_DCT:
330 if (do_inverse)
331 av_log(NULL, AV_LOG_INFO,"DCT_III");
332 else
333 av_log(NULL, AV_LOG_INFO,"DCT_II");
334 ff_dct_init(&d, fft_nbits, do_inverse ? DCT_III : DCT_II);
335 break;
336 #endif /* CONFIG_DCT */
337 #endif /* FFT_FLOAT */
338 default:
339 av_log(NULL, AV_LOG_ERROR, "Requested transform not supported\n");
340 goto cleanup;
341 }
342 av_log(NULL, AV_LOG_INFO," %d test\n", fft_size);
343
344 /* generate random data */
345
346 for (i = 0; i < fft_size; i++) {
347 tab1[i].re = frandom(&prng);
348 tab1[i].im = frandom(&prng);
349 }
350
351 /* checking result */
352 av_log(NULL, AV_LOG_INFO,"Checking...\n");
353
354 switch (transform) {
355 #if CONFIG_MDCT
356 case TRANSFORM_MDCT:
357 if (do_inverse) {
358 imdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits);
359 m.imdct_calc(&m, tab2, (FFTSample *)tab1);
360 err = check_diff((FFTSample *)tab_ref, tab2, fft_size, scale);
361 } else {
362 mdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits);
363
364 m.mdct_calc(&m, tab2, (FFTSample *)tab1);
365
366 err = check_diff((FFTSample *)tab_ref, tab2, fft_size / 2, scale);
367 }
368 break;
369 #endif /* CONFIG_MDCT */
370 case TRANSFORM_FFT:
371 memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
372 s.fft_permute(&s, tab);
373 s.fft_calc(&s, tab);
374
375 fft_ref(tab_ref, tab1, fft_nbits);
376 err = check_diff((FFTSample *)tab_ref, (FFTSample *)tab, fft_size * 2, 1.0);
377 break;
378 #if FFT_FLOAT
379 #if CONFIG_RDFT
380 case TRANSFORM_RDFT:
381 {
382 int fft_size_2 = fft_size >> 1;
383 if (do_inverse) {
384 tab1[ 0].im = 0;
385 tab1[fft_size_2].im = 0;
386 for (i = 1; i < fft_size_2; i++) {
387 tab1[fft_size_2+i].re = tab1[fft_size_2-i].re;
388 tab1[fft_size_2+i].im = -tab1[fft_size_2-i].im;
389 }
390
391 memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
392 tab2[1] = tab1[fft_size_2].re;
393
394 r.rdft_calc(&r, tab2);
395 fft_ref(tab_ref, tab1, fft_nbits);
396 for (i = 0; i < fft_size; i++) {
397 tab[i].re = tab2[i];
398 tab[i].im = 0;
399 }
400 err = check_diff((float *)tab_ref, (float *)tab, fft_size * 2, 0.5);
401 } else {
402 for (i = 0; i < fft_size; i++) {
403 tab2[i] = tab1[i].re;
404 tab1[i].im = 0;
405 }
406 r.rdft_calc(&r, tab2);
407 fft_ref(tab_ref, tab1, fft_nbits);
408 tab_ref[0].im = tab_ref[fft_size_2].re;
409 err = check_diff((float *)tab_ref, (float *)tab2, fft_size, 1.0);
410 }
411 break;
412 }
413 #endif /* CONFIG_RDFT */
414 #if CONFIG_DCT
415 case TRANSFORM_DCT:
416 memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
417 d.dct_calc(&d, tab);
418 if (do_inverse) {
419 idct_ref(tab_ref, tab1, fft_nbits);
420 } else {
421 dct_ref(tab_ref, tab1, fft_nbits);
422 }
423 err = check_diff((float *)tab_ref, (float *)tab, fft_size, 1.0);
424 break;
425 #endif /* CONFIG_DCT */
426 #endif /* FFT_FLOAT */
427 }
428
429 /* do a speed test */
430
431 if (do_speed) {
432 int64_t time_start, duration;
433 int nb_its;
434
435 av_log(NULL, AV_LOG_INFO,"Speed test...\n");
436 /* we measure during about 1 seconds */
437 nb_its = 1;
438 for(;;) {
439 time_start = av_gettime();
440 for (it = 0; it < nb_its; it++) {
441 switch (transform) {
442 case TRANSFORM_MDCT:
443 if (do_inverse) {
444 m.imdct_calc(&m, (FFTSample *)tab, (FFTSample *)tab1);
445 } else {
446 m.mdct_calc(&m, (FFTSample *)tab, (FFTSample *)tab1);
447 }
448 break;
449 case TRANSFORM_FFT:
450 memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
451 s.fft_calc(&s, tab);
452 break;
453 #if FFT_FLOAT
454 case TRANSFORM_RDFT:
455 memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
456 r.rdft_calc(&r, tab2);
457 break;
458 case TRANSFORM_DCT:
459 memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
460 d.dct_calc(&d, tab2);
461 break;
462 #endif /* FFT_FLOAT */
463 }
464 }
465 duration = av_gettime() - time_start;
466 if (duration >= 1000000)
467 break;
468 nb_its *= 2;
469 }
470 av_log(NULL, AV_LOG_INFO,"time: %0.1f us/transform [total time=%0.2f s its=%d]\n",
471 (double)duration / nb_its,
472 (double)duration / 1000000.0,
473 nb_its);
474 }
475
476 switch (transform) {
477 #if CONFIG_MDCT
478 case TRANSFORM_MDCT:
479 ff_mdct_end(&m);
480 break;
481 #endif /* CONFIG_MDCT */
482 case TRANSFORM_FFT:
483 ff_fft_end(&s);
484 break;
485 #if FFT_FLOAT
486 #if CONFIG_RDFT
487 case TRANSFORM_RDFT:
488 ff_rdft_end(&r);
489 break;
490 #endif /* CONFIG_RDFT */
491 #if CONFIG_DCT
492 case TRANSFORM_DCT:
493 ff_dct_end(&d);
494 break;
495 #endif /* CONFIG_DCT */
496 #endif /* FFT_FLOAT */
497 }
498
499 cleanup:
500 av_free(tab);
501 av_free(tab1);
502 av_free(tab2);
503 av_free(tab_ref);
504 av_free(exptab);
505
506 if (err)
507 printf("Error: %d.\n", err);
508
509 return !!err;
510 }