1 /*
2 * (c) 2002 Fabrice Bellard
3 *
4 * This file is part of FFmpeg.
5 *
6 * FFmpeg 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 * FFmpeg 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 FFmpeg; 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 "config.h"
27
28 #ifndef AVFFT
29 #define AVFFT 0
30 #endif
31
32 #include <math.h>
33 #if HAVE_UNISTD_H
34 #include <unistd.h>
35 #endif
36 #include <stdio.h>
37 #include <stdlib.h>
38 #include <string.h>
39
40 #include "libavutil/cpu.h"
41 #include "libavutil/error.h"
42 #include "libavutil/lfg.h"
43 #include "libavutil/log.h"
44 #include "libavutil/mathematics.h"
45 #include "libavutil/time.h"
46
47 #if AVFFT
48 #include "libavcodec/avfft.h"
49 #else
50 #include "libavcodec/fft.h"
51 #endif
52
53 #if FFT_FLOAT
54 #include "libavcodec/dct.h"
55 #include "libavcodec/rdft.h"
56 #endif
57
58 /* reference fft */
59
60 #define MUL16(a, b) ((a) * (b))
61
62 #define CMAC(pre, pim, are, aim, bre, bim) \
63 { \
64 pre += (MUL16(are, bre) - MUL16(aim, bim)); \
65 pim += (MUL16(are, bim) + MUL16(bre, aim)); \
66 }
67
68 #if FFT_FLOAT || AVFFT
69 #define RANGE 1.0
70 #define REF_SCALE(x, bits) (x)
71 #define FMT "%10.6f"
72 #else
73 #define RANGE 8388608
74 #define REF_SCALE(x, bits) (x)
75 #define FMT "%6d"
76 #endif
77
78 static struct {
79 float re, im;
80 } *exptab;
81
fft_ref_init(int nbits, int inverse)82 static int fft_ref_init(int nbits, int inverse)
83 {
84 int i, n = 1 << nbits;
85
86 exptab = av_malloc_array((n / 2), sizeof(*exptab));
87 if (!exptab)
88 return AVERROR(ENOMEM);
89
90 for (i = 0; i < (n / 2); i++) {
91 double alpha = 2 * M_PI * (float) i / (float) n;
92 double c1 = cos(alpha), s1 = sin(alpha);
93 if (!inverse)
94 s1 = -s1;
95 exptab[i].re = c1;
96 exptab[i].im = s1;
97 }
98 return 0;
99 }
100
fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits)101 static void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits)
102 {
103 int i, j;
104 int n = 1 << nbits;
105 int n2 = n >> 1;
106
107 for (i = 0; i < n; i++) {
108 double tmp_re = 0, tmp_im = 0;
109 FFTComplex *q = tab;
110 for (j = 0; j < n; j++) {
111 double s, c;
112 int k = (i * j) & (n - 1);
113 if (k >= n2) {
114 c = -exptab[k - n2].re;
115 s = -exptab[k - n2].im;
116 } else {
117 c = exptab[k].re;
118 s = exptab[k].im;
119 }
120 CMAC(tmp_re, tmp_im, c, s, q->re, q->im);
121 q++;
122 }
123 tabr[i].re = REF_SCALE(tmp_re, nbits);
124 tabr[i].im = REF_SCALE(tmp_im, nbits);
125 }
126 }
127
128 #if CONFIG_MDCT
imdct_ref(FFTSample *out, FFTSample *in, int nbits)129 static void imdct_ref(FFTSample *out, FFTSample *in, int nbits)
130 {
131 int i, k, n = 1 << nbits;
132
133 for (i = 0; i < n; i++) {
134 double sum = 0;
135 for (k = 0; k < n / 2; k++) {
136 int a = (2 * i + 1 + (n / 2)) * (2 * k + 1);
137 double f = cos(M_PI * a / (double) (2 * n));
138 sum += f * in[k];
139 }
140 out[i] = REF_SCALE(-sum, nbits - 2);
141 }
142 }
143
144 /* NOTE: no normalisation by 1 / N is done */
mdct_ref(FFTSample *output, FFTSample *input, int nbits)145 static void mdct_ref(FFTSample *output, FFTSample *input, int nbits)
146 {
147 int i, k, n = 1 << nbits;
148
149 /* do it by hand */
150 for (k = 0; k < n / 2; k++) {
151 double s = 0;
152 for (i = 0; i < n; i++) {
153 double a = (2 * M_PI * (2 * i + 1 + n / 2) * (2 * k + 1) / (4 * n));
154 s += input[i] * cos(a);
155 }
156 output[k] = REF_SCALE(s, nbits - 1);
157 }
158 }
159 #endif /* CONFIG_MDCT */
160
161 #if FFT_FLOAT
162 #if CONFIG_DCT
idct_ref(FFTSample *output, FFTSample *input, int nbits)163 static void idct_ref(FFTSample *output, FFTSample *input, int nbits)
164 {
165 int i, k, n = 1 << nbits;
166
167 /* do it by hand */
168 for (i = 0; i < n; i++) {
169 double s = 0.5 * input[0];
170 for (k = 1; k < n; k++) {
171 double a = M_PI * k * (i + 0.5) / n;
172 s += input[k] * cos(a);
173 }
174 output[i] = 2 * s / n;
175 }
176 }
177
dct_ref(FFTSample *output, FFTSample *input, int nbits)178 static void dct_ref(FFTSample *output, FFTSample *input, int nbits)
179 {
180 int i, k, n = 1 << nbits;
181
182 /* do it by hand */
183 for (k = 0; k < n; k++) {
184 double s = 0;
185 for (i = 0; i < n; i++) {
186 double a = M_PI * k * (i + 0.5) / n;
187 s += input[i] * cos(a);
188 }
189 output[k] = s;
190 }
191 }
192 #endif /* CONFIG_DCT */
193 #endif /* FFT_FLOAT */
194
frandom(AVLFG *prng)195 static FFTSample frandom(AVLFG *prng)
196 {
197 return (int16_t) av_lfg_get(prng) / 32768.0 * RANGE;
198 }
199
check_diff(FFTSample *tab1, FFTSample *tab2, int n, double scale)200 static int check_diff(FFTSample *tab1, FFTSample *tab2, int n, double scale)
201 {
202 int i, err = 0;
203 double error = 0, max = 0;
204
205 for (i = 0; i < n; i++) {
206 double e = fabs(tab1[i] - (tab2[i] / scale)) / RANGE;
207 if (e >= 1e-3) {
208 av_log(NULL, AV_LOG_ERROR, "ERROR %5d: "FMT" "FMT"\n",
209 i, tab1[i], tab2[i]);
210 err = 1;
211 }
212 error += e * e;
213 if (e > max)
214 max = e;
215 }
216 av_log(NULL, AV_LOG_INFO, "max:%f e:%g\n", max, sqrt(error / n));
217 return err;
218 }
219
fft_init(FFTContext **s, int nbits, int inverse)220 static inline void fft_init(FFTContext **s, int nbits, int inverse)
221 {
222 #if AVFFT
223 *s = av_fft_init(nbits, inverse);
224 #else
225 ff_fft_init(*s, nbits, inverse);
226 #endif
227 }
228
mdct_init(FFTContext **s, int nbits, int inverse, double scale)229 static inline void mdct_init(FFTContext **s, int nbits, int inverse, double scale)
230 {
231 #if AVFFT
232 *s = av_mdct_init(nbits, inverse, scale);
233 #else
234 ff_mdct_init(*s, nbits, inverse, scale);
235 #endif
236 }
237
mdct_calc(FFTContext *s, FFTSample *output, const FFTSample *input)238 static inline void mdct_calc(FFTContext *s, FFTSample *output, const FFTSample *input)
239 {
240 #if AVFFT
241 av_mdct_calc(s, output, input);
242 #else
243 s->mdct_calc(s, output, input);
244 #endif
245 }
246
imdct_calc(struct FFTContext *s, FFTSample *output, const FFTSample *input)247 static inline void imdct_calc(struct FFTContext *s, FFTSample *output, const FFTSample *input)
248 {
249 #if AVFFT
250 av_imdct_calc(s, output, input);
251 #else
252 s->imdct_calc(s, output, input);
253 #endif
254 }
255
fft_permute(FFTContext *s, FFTComplex *z)256 static inline void fft_permute(FFTContext *s, FFTComplex *z)
257 {
258 #if AVFFT
259 av_fft_permute(s, z);
260 #else
261 s->fft_permute(s, z);
262 #endif
263 }
264
fft_calc(FFTContext *s, FFTComplex *z)265 static inline void fft_calc(FFTContext *s, FFTComplex *z)
266 {
267 #if AVFFT
268 av_fft_calc(s, z);
269 #else
270 s->fft_calc(s, z);
271 #endif
272 }
273
mdct_end(FFTContext *s)274 static inline void mdct_end(FFTContext *s)
275 {
276 #if AVFFT
277 av_mdct_end(s);
278 #else
279 ff_mdct_end(s);
280 #endif
281 }
282
fft_end(FFTContext *s)283 static inline void fft_end(FFTContext *s)
284 {
285 #if AVFFT
286 av_fft_end(s);
287 #else
288 ff_fft_end(s);
289 #endif
290 }
291
292 #if FFT_FLOAT
rdft_init(RDFTContext **r, int nbits, enum RDFTransformType trans)293 static inline void rdft_init(RDFTContext **r, int nbits, enum RDFTransformType trans)
294 {
295 #if AVFFT
296 *r = av_rdft_init(nbits, trans);
297 #else
298 ff_rdft_init(*r, nbits, trans);
299 #endif
300 }
301
dct_init(DCTContext **d, int nbits, enum DCTTransformType trans)302 static inline void dct_init(DCTContext **d, int nbits, enum DCTTransformType trans)
303 {
304 #if AVFFT
305 *d = av_dct_init(nbits, trans);
306 #else
307 ff_dct_init(*d, nbits, trans);
308 #endif
309 }
310
rdft_calc(RDFTContext *r, FFTSample *tab)311 static inline void rdft_calc(RDFTContext *r, FFTSample *tab)
312 {
313 #if AVFFT
314 av_rdft_calc(r, tab);
315 #else
316 r->rdft_calc(r, tab);
317 #endif
318 }
319
dct_calc(DCTContext *d, FFTSample *data)320 static inline void dct_calc(DCTContext *d, FFTSample *data)
321 {
322 #if AVFFT
323 av_dct_calc(d, data);
324 #else
325 d->dct_calc(d, data);
326 #endif
327 }
328
rdft_end(RDFTContext *r)329 static inline void rdft_end(RDFTContext *r)
330 {
331 #if AVFFT
332 av_rdft_end(r);
333 #else
334 ff_rdft_end(r);
335 #endif
336 }
337
dct_end(DCTContext *d)338 static inline void dct_end(DCTContext *d)
339 {
340 #if AVFFT
341 av_dct_end(d);
342 #else
343 ff_dct_end(d);
344 #endif
345 }
346 #endif /* FFT_FLOAT */
347
help(void)348 static void help(void)
349 {
350 av_log(NULL, AV_LOG_INFO,
351 "usage: fft-test [-h] [-s] [-i] [-n b]\n"
352 "-h print this help\n"
353 "-s speed test\n"
354 "-m (I)MDCT test\n"
355 "-d (I)DCT test\n"
356 "-r (I)RDFT test\n"
357 "-i inverse transform test\n"
358 "-n b set the transform size to 2^b\n"
359 "-f x set scale factor for output data of (I)MDCT to x\n");
360 }
361
362 enum tf_transform {
363 TRANSFORM_FFT,
364 TRANSFORM_MDCT,
365 TRANSFORM_RDFT,
366 TRANSFORM_DCT,
367 };
368
369 #if !HAVE_GETOPT
370 #include "compat/getopt.c"
371 #endif
372
main(int argc, char **argv)373 int main(int argc, char **argv)
374 {
375 FFTComplex *tab, *tab1, *tab_ref;
376 FFTSample *tab2;
377 enum tf_transform transform = TRANSFORM_FFT;
378 FFTContext *m, *s;
379 #if FFT_FLOAT
380 RDFTContext *r;
381 DCTContext *d;
382 #endif /* FFT_FLOAT */
383 int it, i, err = 1;
384 int do_speed = 0, do_inverse = 0;
385 int fft_nbits = 9, fft_size;
386 double scale = 1.0;
387 AVLFG prng;
388
389 #if !AVFFT
390 s = av_mallocz(sizeof(*s));
391 m = av_mallocz(sizeof(*m));
392 #endif
393
394 #if !AVFFT && FFT_FLOAT
395 r = av_mallocz(sizeof(*r));
396 d = av_mallocz(sizeof(*d));
397 #endif
398
399 av_lfg_init(&prng, 1);
400
401 for (;;) {
402 int c = getopt(argc, argv, "hsimrdn:f:c:");
403 if (c == -1)
404 break;
405 switch (c) {
406 case 'h':
407 help();
408 return 1;
409 case 's':
410 do_speed = 1;
411 break;
412 case 'i':
413 do_inverse = 1;
414 break;
415 case 'm':
416 transform = TRANSFORM_MDCT;
417 break;
418 case 'r':
419 transform = TRANSFORM_RDFT;
420 break;
421 case 'd':
422 transform = TRANSFORM_DCT;
423 break;
424 case 'n':
425 fft_nbits = atoi(optarg);
426 break;
427 case 'f':
428 scale = atof(optarg);
429 break;
430 case 'c':
431 {
432 unsigned cpuflags = av_get_cpu_flags();
433
434 if (av_parse_cpu_caps(&cpuflags, optarg) < 0)
435 return 1;
436
437 av_force_cpu_flags(cpuflags);
438 break;
439 }
440 }
441 }
442
443 fft_size = 1 << fft_nbits;
444 tab = av_malloc_array(fft_size, sizeof(FFTComplex));
445 tab1 = av_malloc_array(fft_size, sizeof(FFTComplex));
446 tab_ref = av_malloc_array(fft_size, sizeof(FFTComplex));
447 tab2 = av_malloc_array(fft_size, sizeof(FFTSample));
448
449 if (!(tab && tab1 && tab_ref && tab2))
450 goto cleanup;
451
452 switch (transform) {
453 #if CONFIG_MDCT
454 case TRANSFORM_MDCT:
455 av_log(NULL, AV_LOG_INFO, "Scale factor is set to %f\n", scale);
456 if (do_inverse)
457 av_log(NULL, AV_LOG_INFO, "IMDCT");
458 else
459 av_log(NULL, AV_LOG_INFO, "MDCT");
460 mdct_init(&m, fft_nbits, do_inverse, scale);
461 break;
462 #endif /* CONFIG_MDCT */
463 case TRANSFORM_FFT:
464 if (do_inverse)
465 av_log(NULL, AV_LOG_INFO, "IFFT");
466 else
467 av_log(NULL, AV_LOG_INFO, "FFT");
468 fft_init(&s, fft_nbits, do_inverse);
469 if ((err = fft_ref_init(fft_nbits, do_inverse)) < 0)
470 goto cleanup;
471 break;
472 #if FFT_FLOAT
473 # if CONFIG_RDFT
474 case TRANSFORM_RDFT:
475 if (do_inverse)
476 av_log(NULL, AV_LOG_INFO, "IDFT_C2R");
477 else
478 av_log(NULL, AV_LOG_INFO, "DFT_R2C");
479 rdft_init(&r, fft_nbits, do_inverse ? IDFT_C2R : DFT_R2C);
480 if ((err = fft_ref_init(fft_nbits, do_inverse)) < 0)
481 goto cleanup;
482 break;
483 # endif /* CONFIG_RDFT */
484 # if CONFIG_DCT
485 case TRANSFORM_DCT:
486 if (do_inverse)
487 av_log(NULL, AV_LOG_INFO, "DCT_III");
488 else
489 av_log(NULL, AV_LOG_INFO, "DCT_II");
490 dct_init(&d, fft_nbits, do_inverse ? DCT_III : DCT_II);
491 break;
492 # endif /* CONFIG_DCT */
493 #endif /* FFT_FLOAT */
494 default:
495 av_log(NULL, AV_LOG_ERROR, "Requested transform not supported\n");
496 goto cleanup;
497 }
498 av_log(NULL, AV_LOG_INFO, " %d test\n", fft_size);
499
500 /* generate random data */
501
502 for (i = 0; i < fft_size; i++) {
503 tab1[i].re = frandom(&prng);
504 tab1[i].im = frandom(&prng);
505 }
506
507 /* checking result */
508 av_log(NULL, AV_LOG_INFO, "Checking...\n");
509
510 switch (transform) {
511 #if CONFIG_MDCT
512 case TRANSFORM_MDCT:
513 if (do_inverse) {
514 imdct_ref(&tab_ref->re, &tab1->re, fft_nbits);
515 imdct_calc(m, tab2, &tab1->re);
516 err = check_diff(&tab_ref->re, tab2, fft_size, scale);
517 } else {
518 mdct_ref(&tab_ref->re, &tab1->re, fft_nbits);
519 mdct_calc(m, tab2, &tab1->re);
520 err = check_diff(&tab_ref->re, tab2, fft_size / 2, scale);
521 }
522 break;
523 #endif /* CONFIG_MDCT */
524 case TRANSFORM_FFT:
525 memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
526 fft_permute(s, tab);
527 fft_calc(s, tab);
528
529 fft_ref(tab_ref, tab1, fft_nbits);
530 err = check_diff(&tab_ref->re, &tab->re, fft_size * 2, 1.0);
531 break;
532 #if FFT_FLOAT
533 #if CONFIG_RDFT
534 case TRANSFORM_RDFT:
535 {
536 int fft_size_2 = fft_size >> 1;
537 if (do_inverse) {
538 tab1[0].im = 0;
539 tab1[fft_size_2].im = 0;
540 for (i = 1; i < fft_size_2; i++) {
541 tab1[fft_size_2 + i].re = tab1[fft_size_2 - i].re;
542 tab1[fft_size_2 + i].im = -tab1[fft_size_2 - i].im;
543 }
544
545 memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
546 tab2[1] = tab1[fft_size_2].re;
547
548 rdft_calc(r, tab2);
549 fft_ref(tab_ref, tab1, fft_nbits);
550 for (i = 0; i < fft_size; i++) {
551 tab[i].re = tab2[i];
552 tab[i].im = 0;
553 }
554 err = check_diff(&tab_ref->re, &tab->re, fft_size * 2, 0.5);
555 } else {
556 for (i = 0; i < fft_size; i++) {
557 tab2[i] = tab1[i].re;
558 tab1[i].im = 0;
559 }
560 rdft_calc(r, tab2);
561 fft_ref(tab_ref, tab1, fft_nbits);
562 tab_ref[0].im = tab_ref[fft_size_2].re;
563 err = check_diff(&tab_ref->re, tab2, fft_size, 1.0);
564 }
565 break;
566 }
567 #endif /* CONFIG_RDFT */
568 #if CONFIG_DCT
569 case TRANSFORM_DCT:
570 memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
571 dct_calc(d, &tab->re);
572 if (do_inverse)
573 idct_ref(&tab_ref->re, &tab1->re, fft_nbits);
574 else
575 dct_ref(&tab_ref->re, &tab1->re, fft_nbits);
576 err = check_diff(&tab_ref->re, &tab->re, fft_size, 1.0);
577 break;
578 #endif /* CONFIG_DCT */
579 #endif /* FFT_FLOAT */
580 }
581
582 /* do a speed test */
583
584 if (do_speed) {
585 int64_t time_start, duration;
586 int nb_its;
587
588 av_log(NULL, AV_LOG_INFO, "Speed test...\n");
589 /* we measure during about 1 seconds */
590 nb_its = 1;
591 for (;;) {
592 time_start = av_gettime_relative();
593 for (it = 0; it < nb_its; it++) {
594 switch (transform) {
595 case TRANSFORM_MDCT:
596 if (do_inverse)
597 imdct_calc(m, &tab->re, &tab1->re);
598 else
599 mdct_calc(m, &tab->re, &tab1->re);
600 break;
601 case TRANSFORM_FFT:
602 memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
603 fft_calc(s, tab);
604 break;
605 #if FFT_FLOAT
606 case TRANSFORM_RDFT:
607 memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
608 rdft_calc(r, tab2);
609 break;
610 case TRANSFORM_DCT:
611 memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
612 dct_calc(d, tab2);
613 break;
614 #endif /* FFT_FLOAT */
615 }
616 }
617 duration = av_gettime_relative() - time_start;
618 if (duration >= 1000000)
619 break;
620 nb_its *= 2;
621 }
622 av_log(NULL, AV_LOG_INFO,
623 "time: %0.1f us/transform [total time=%0.2f s its=%d]\n",
624 (double) duration / nb_its,
625 (double) duration / 1000000.0,
626 nb_its);
627 }
628
629 switch (transform) {
630 #if CONFIG_MDCT
631 case TRANSFORM_MDCT:
632 mdct_end(m);
633 break;
634 #endif /* CONFIG_MDCT */
635 case TRANSFORM_FFT:
636 fft_end(s);
637 break;
638 #if FFT_FLOAT
639 # if CONFIG_RDFT
640 case TRANSFORM_RDFT:
641 rdft_end(r);
642 break;
643 # endif /* CONFIG_RDFT */
644 # if CONFIG_DCT
645 case TRANSFORM_DCT:
646 dct_end(d);
647 break;
648 # endif /* CONFIG_DCT */
649 #endif /* FFT_FLOAT */
650 }
651
652 cleanup:
653 av_free(tab);
654 av_free(tab1);
655 av_free(tab2);
656 av_free(tab_ref);
657 av_free(exptab);
658
659 #if !AVFFT
660 av_free(s);
661 av_free(m);
662 #endif
663
664 #if !AVFFT && FFT_FLOAT
665 av_free(r);
666 av_free(d);
667 #endif
668
669 if (err)
670 printf("Error: %d.\n", err);
671
672 return !!err;
673 }
674