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