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 68struct { 69 float re, im; 70} *exptab; 71 72static 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 91static 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 121static 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 */ 139static 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 159static 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} 175static 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 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; 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 222static 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 236enum 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 247int 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} 526