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 78static struct { 79 float re, im; 80} *exptab; 81 82static 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 101static 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 129static 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 */ 145static 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 163static 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 178static 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 195static FFTSample frandom(AVLFG *prng) 196{ 197 return (int16_t) av_lfg_get(prng) / 32768.0 * RANGE; 198} 199 200static 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 220static 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 229static 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 238static 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 247static 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 256static 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 265static 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 274static 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 283static 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 293static 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 302static 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 311static 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 320static 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 329static 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 338static 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 348static 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 362enum 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 373int 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 652cleanup: 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