git.sesse.net Git - ffmpeg/blob - libavcodec/tests/fft.c

   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
  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
 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
 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 */
 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
 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
 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
 198 static FFTSample frandom(AVLFG *prng)
 199 {
 200     return (int16_t) av_lfg_get(prng) / 32768.0 * RANGE;
 201 }
 202
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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 }