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