libavcodec/dct-test.c

   1 /*
   2  * (c) 2001 Fabrice Bellard
   3  *     2007 Marc Hoffman <marc.hoffman@analog.com>
   4  *
   5  * This file is part of FFmpeg.
   6  *
   7  * FFmpeg is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU Lesser General Public
   9  * License as published by the Free Software Foundation; either
  10  * version 2.1 of the License, or (at your option) any later version.
  11  *
  12  * FFmpeg is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15  * Lesser General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU Lesser General Public
  18  * License along with FFmpeg; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  20  */
  21
  22 /**
  23  * @file libavcodec/dct-test.c
  24  * DCT test (c) 2001 Fabrice Bellard
  25  * Started from sample code by Juan J. Sierralta P.
  26  */
  27
  28 #include <stdlib.h>
  29 #include <stdio.h>
  30 #include <string.h>
  31 #include <sys/time.h>
  32 #include <unistd.h>
  33 #include <math.h>
  34
  35 #include "libavutil/common.h"
  36 #include "libavutil/lfg.h"
  37
  38 #include "simple_idct.h"
  39 #include "aandcttab.h"
  40 #include "faandct.h"
  41 #include "faanidct.h"
  42 #include "x86/idct_xvid.h"
  43
  44 #undef printf
  45
  46 void *fast_memcpy(void *a, const void *b, size_t c){return memcpy(a,b,c);};
  47
  48 /* reference fdct/idct */
  49 void ff_ref_fdct(DCTELEM *block);
  50 void ff_ref_idct(DCTELEM *block);
  51 void ff_ref_dct_init(void);
  52
  53 void ff_mmx_idct(DCTELEM *data);
  54 void ff_mmxext_idct(DCTELEM *data);
  55
  56 void odivx_idct_c(short *block);
  57
  58 // BFIN
  59 void ff_bfin_idct(DCTELEM *block);
  60 void ff_bfin_fdct(DCTELEM *block);
  61
  62 // ALTIVEC
  63 void fdct_altivec(DCTELEM *block);
  64 //void idct_altivec(DCTELEM *block);?? no routine
  65
  66 // ARM
  67 void j_rev_dct_ARM(DCTELEM *data);
  68 void simple_idct_ARM(DCTELEM *data);
  69 void simple_idct_armv5te(DCTELEM *data);
  70 void ff_simple_idct_armv6(DCTELEM *data);
  71 void ff_simple_idct_neon(DCTELEM *data);
  72
  73 void ff_simple_idct_axp(DCTELEM *data);
  74
  75 struct algo {
  76   const char *name;
  77   enum { FDCT, IDCT } is_idct;
  78   void (* func) (DCTELEM *block);
  79   void (* ref)  (DCTELEM *block);
  80   enum formattag { NO_PERM,MMX_PERM, MMX_SIMPLE_PERM, SCALE_PERM, SSE2_PERM, PARTTRANS_PERM } format;
  81   int  mm_support;
  82 };
  83
  84 #ifndef FAAN_POSTSCALE
  85 #define FAAN_SCALE SCALE_PERM
  86 #else
  87 #define FAAN_SCALE NO_PERM
  88 #endif
  89
  90 static int cpu_flags;
  91
  92 struct algo algos[] = {
  93   {"REF-DBL",         0, ff_ref_fdct,        ff_ref_fdct, NO_PERM},
  94   {"FAAN",            0, ff_faandct,         ff_ref_fdct, FAAN_SCALE},
  95   {"FAANI",           1, ff_faanidct,        ff_ref_idct, NO_PERM},
  96   {"IJG-AAN-INT",     0, fdct_ifast,         ff_ref_fdct, SCALE_PERM},
  97   {"IJG-LLM-INT",     0, ff_jpeg_fdct_islow, ff_ref_fdct, NO_PERM},
  98   {"REF-DBL",         1, ff_ref_idct,        ff_ref_idct, NO_PERM},
  99   {"INT",             1, j_rev_dct,          ff_ref_idct, MMX_PERM},
 100   {"SIMPLE-C",        1, ff_simple_idct,     ff_ref_idct, NO_PERM},
 101
 102 #if HAVE_MMX
 103   {"MMX",             0, ff_fdct_mmx,        ff_ref_fdct, NO_PERM, FF_MM_MMX},
 104 #if HAVE_MMX2
 105   {"MMX2",            0, ff_fdct_mmx2,       ff_ref_fdct, NO_PERM, FF_MM_MMX2},
 106   {"SSE2",            0, ff_fdct_sse2,       ff_ref_fdct, NO_PERM, FF_MM_SSE2},
 107 #endif
 108
 109 #if CONFIG_GPL
 110   {"LIBMPEG2-MMX",    1, ff_mmx_idct,        ff_ref_idct, MMX_PERM, FF_MM_MMX},
 111   {"LIBMPEG2-MMX2",   1, ff_mmxext_idct,     ff_ref_idct, MMX_PERM, FF_MM_MMX2},
 112 #endif
 113   {"SIMPLE-MMX",      1, ff_simple_idct_mmx, ff_ref_idct, MMX_SIMPLE_PERM, FF_MM_MMX},
 114   {"XVID-MMX",        1, ff_idct_xvid_mmx,   ff_ref_idct, NO_PERM, FF_MM_MMX},
 115   {"XVID-MMX2",       1, ff_idct_xvid_mmx2,  ff_ref_idct, NO_PERM, FF_MM_MMX2},
 116   {"XVID-SSE2",       1, ff_idct_xvid_sse2,  ff_ref_idct, SSE2_PERM, FF_MM_SSE2},
 117 #endif
 118
 119 #if HAVE_ALTIVEC
 120   {"altivecfdct",     0, fdct_altivec,       ff_ref_fdct, NO_PERM, FF_MM_ALTIVEC},
 121 #endif
 122
 123 #if ARCH_BFIN
 124   {"BFINfdct",        0, ff_bfin_fdct,       ff_ref_fdct, NO_PERM},
 125   {"BFINidct",        1, ff_bfin_idct,       ff_ref_idct, NO_PERM},
 126 #endif
 127
 128 #if ARCH_ARM
 129   {"SIMPLE-ARM",      1, simple_idct_ARM,    ff_ref_idct, NO_PERM },
 130   {"INT-ARM",         1, j_rev_dct_ARM,      ff_ref_idct, MMX_PERM },
 131 #if HAVE_ARMV5TE
 132   {"SIMPLE-ARMV5TE",  1, simple_idct_armv5te, ff_ref_idct, NO_PERM },
 133 #endif
 134 #if HAVE_ARMV6
 135   {"SIMPLE-ARMV6",    1, ff_simple_idct_armv6, ff_ref_idct, MMX_PERM },
 136 #endif
 137 #if HAVE_NEON
 138   {"SIMPLE-NEON",     1, ff_simple_idct_neon, ff_ref_idct, PARTTRANS_PERM },
 139 #endif
 140 #endif /* ARCH_ARM */
 141
 142 #if ARCH_ALPHA
 143   {"SIMPLE-ALPHA",    1, ff_simple_idct_axp,  ff_ref_idct, NO_PERM },
 144 #endif
 145
 146   { 0 }
 147 };
 148
 149 #define AANSCALE_BITS 12
 150
 151 uint8_t cropTbl[256 + 2 * MAX_NEG_CROP];
 152
 153 static int64_t gettime(void)
 154 {
 155     struct timeval tv;
 156     gettimeofday(&tv,NULL);
 157     return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
 158 }
 159
 160 #define NB_ITS 20000
 161 #define NB_ITS_SPEED 50000
 162
 163 static short idct_mmx_perm[64];
 164
 165 static short idct_simple_mmx_perm[64]={
 166         0x00, 0x08, 0x04, 0x09, 0x01, 0x0C, 0x05, 0x0D,
 167         0x10, 0x18, 0x14, 0x19, 0x11, 0x1C, 0x15, 0x1D,
 168         0x20, 0x28, 0x24, 0x29, 0x21, 0x2C, 0x25, 0x2D,
 169         0x12, 0x1A, 0x16, 0x1B, 0x13, 0x1E, 0x17, 0x1F,
 170         0x02, 0x0A, 0x06, 0x0B, 0x03, 0x0E, 0x07, 0x0F,
 171         0x30, 0x38, 0x34, 0x39, 0x31, 0x3C, 0x35, 0x3D,
 172         0x22, 0x2A, 0x26, 0x2B, 0x23, 0x2E, 0x27, 0x2F,
 173         0x32, 0x3A, 0x36, 0x3B, 0x33, 0x3E, 0x37, 0x3F,
 174 };
 175
 176 static const uint8_t idct_sse2_row_perm[8] = {0, 4, 1, 5, 2, 6, 3, 7};
 177
 178 static void idct_mmx_init(void)
 179 {
 180     int i;
 181
 182     /* the mmx/mmxext idct uses a reordered input, so we patch scan tables */
 183     for (i = 0; i < 64; i++) {
 184         idct_mmx_perm[i] = (i & 0x38) | ((i & 6) >> 1) | ((i & 1) << 2);
 185 //        idct_simple_mmx_perm[i] = simple_block_permute_op(i);
 186     }
 187 }
 188
 189 static DCTELEM block[64] __attribute__ ((aligned (16)));
 190 static DCTELEM block1[64] __attribute__ ((aligned (8)));
 191 static DCTELEM block_org[64] __attribute__ ((aligned (8)));
 192
 193 static inline void mmx_emms(void)
 194 {
 195 #if HAVE_MMX
 196     if (cpu_flags & FF_MM_MMX)
 197         __asm__ volatile ("emms\n\t");
 198 #endif
 199 }
 200
 201 static void dct_error(const char *name, int is_idct,
 202                void (*fdct_func)(DCTELEM *block),
 203                void (*fdct_ref)(DCTELEM *block), int form, int test)
 204 {
 205     int it, i, scale;
 206     int err_inf, v;
 207     int64_t err2, ti, ti1, it1;
 208     int64_t sysErr[64], sysErrMax=0;
 209     int maxout=0;
 210     int blockSumErrMax=0, blockSumErr;
 211     AVLFG prng;
 212
 213     av_lfg_init(&prng, 1);
 214
 215     err_inf = 0;
 216     err2 = 0;
 217     for(i=0; i<64; i++) sysErr[i]=0;
 218     for(it=0;it<NB_ITS;it++) {
 219         for(i=0;i<64;i++)
 220             block1[i] = 0;
 221         switch(test){
 222         case 0:
 223             for(i=0;i<64;i++)
 224                 block1[i] = (av_lfg_get(&prng) % 512) -256;
 225             if (is_idct){
 226                 ff_ref_fdct(block1);
 227
 228                 for(i=0;i<64;i++)
 229                     block1[i]>>=3;
 230             }
 231         break;
 232         case 1:{
 233             int num = av_lfg_get(&prng) % 10 + 1;
 234             for(i=0;i<num;i++)
 235                 block1[av_lfg_get(&prng) % 64] = av_lfg_get(&prng) % 512 -256;
 236         }break;
 237         case 2:
 238             block1[0] = av_lfg_get(&prng) % 4096 - 2048;
 239             block1[63]= (block1[0]&1)^1;
 240         break;
 241         }
 242
 243 #if 0 // simulate mismatch control
 244 { int sum=0;
 245         for(i=0;i<64;i++)
 246            sum+=block1[i];
 247
 248         if((sum&1)==0) block1[63]^=1;
 249 }
 250 #endif
 251
 252         for(i=0; i<64; i++)
 253             block_org[i]= block1[i];
 254
 255         if (form == MMX_PERM) {
 256             for(i=0;i<64;i++)
 257                 block[idct_mmx_perm[i]] = block1[i];
 258             } else if (form == MMX_SIMPLE_PERM) {
 259             for(i=0;i<64;i++)
 260                 block[idct_simple_mmx_perm[i]] = block1[i];
 261
 262         } else if (form == SSE2_PERM) {
 263             for(i=0; i<64; i++)
 264                 block[(i&0x38) | idct_sse2_row_perm[i&7]] = block1[i];
 265         } else if (form == PARTTRANS_PERM) {
 266             for(i=0; i<64; i++)
 267                 block[(i&0x24) | ((i&3)<<3) | ((i>>3)&3)] = block1[i];
 268         } else {
 269             for(i=0; i<64; i++)
 270                 block[i]= block1[i];
 271         }
 272 #if 0 // simulate mismatch control for tested IDCT but not the ref
 273 { int sum=0;
 274         for(i=0;i<64;i++)
 275            sum+=block[i];
 276
 277         if((sum&1)==0) block[63]^=1;
 278 }
 279 #endif
 280
 281         fdct_func(block);
 282         mmx_emms();
 283
 284         if (form == SCALE_PERM) {
 285             for(i=0; i<64; i++) {
 286                 scale = 8*(1 << (AANSCALE_BITS + 11)) / ff_aanscales[i];
 287                 block[i] = (block[i] * scale /*+ (1<<(AANSCALE_BITS-1))*/) >> AANSCALE_BITS;
 288             }
 289         }
 290
 291         fdct_ref(block1);
 292
 293         blockSumErr=0;
 294         for(i=0;i<64;i++) {
 295             v = abs(block[i] - block1[i]);
 296             if (v > err_inf)
 297                 err_inf = v;
 298             err2 += v * v;
 299             sysErr[i] += block[i] - block1[i];
 300             blockSumErr += v;
 301             if( abs(block[i])>maxout) maxout=abs(block[i]);
 302         }
 303         if(blockSumErrMax < blockSumErr) blockSumErrMax= blockSumErr;
 304 #if 0 // print different matrix pairs
 305         if(blockSumErr){
 306             printf("\n");
 307             for(i=0; i<64; i++){
 308                 if((i&7)==0) printf("\n");
 309                 printf("%4d ", block_org[i]);
 310             }
 311             for(i=0; i<64; i++){
 312                 if((i&7)==0) printf("\n");
 313                 printf("%4d ", block[i] - block1[i]);
 314             }
 315         }
 316 #endif
 317     }
 318     for(i=0; i<64; i++) sysErrMax= FFMAX(sysErrMax, FFABS(sysErr[i]));
 319
 320 #if 1 // dump systematic errors
 321     for(i=0; i<64; i++){
 322         if(i%8==0) printf("\n");
 323         printf("%7d ", (int)sysErr[i]);
 324     }
 325     printf("\n");
 326 #endif
 327
 328     printf("%s %s: err_inf=%d err2=%0.8f syserr=%0.8f maxout=%d blockSumErr=%d\n",
 329            is_idct ? "IDCT" : "DCT",
 330            name, err_inf, (double)err2 / NB_ITS / 64.0, (double)sysErrMax / NB_ITS, maxout, blockSumErrMax);
 331 #if 1 //Speed test
 332     /* speed test */
 333     for(i=0;i<64;i++)
 334         block1[i] = 0;
 335     switch(test){
 336     case 0:
 337         for(i=0;i<64;i++)
 338             block1[i] = av_lfg_get(&prng) % 512 -256;
 339         if (is_idct){
 340             ff_ref_fdct(block1);
 341
 342             for(i=0;i<64;i++)
 343                 block1[i]>>=3;
 344         }
 345     break;
 346     case 1:{
 347     case 2:
 348         block1[0] = av_lfg_get(&prng) % 512 -256;
 349         block1[1] = av_lfg_get(&prng) % 512 -256;
 350         block1[2] = av_lfg_get(&prng) % 512 -256;
 351         block1[3] = av_lfg_get(&prng) % 512 -256;
 352     }break;
 353     }
 354
 355     if (form == MMX_PERM) {
 356         for(i=0;i<64;i++)
 357             block[idct_mmx_perm[i]] = block1[i];
 358     } else if(form == MMX_SIMPLE_PERM) {
 359         for(i=0;i<64;i++)
 360             block[idct_simple_mmx_perm[i]] = block1[i];
 361     } else {
 362         for(i=0; i<64; i++)
 363             block[i]= block1[i];
 364     }
 365
 366     ti = gettime();
 367     it1 = 0;
 368     do {
 369         for(it=0;it<NB_ITS_SPEED;it++) {
 370             for(i=0; i<64; i++)
 371                 block[i]= block1[i];
 372 //            memcpy(block, block1, sizeof(DCTELEM) * 64);
 373 // do not memcpy especially not fastmemcpy because it does movntq !!!
 374             fdct_func(block);
 375         }
 376         it1 += NB_ITS_SPEED;
 377         ti1 = gettime() - ti;
 378     } while (ti1 < 1000000);
 379     mmx_emms();
 380
 381     printf("%s %s: %0.1f kdct/s\n",
 382            is_idct ? "IDCT" : "DCT",
 383            name, (double)it1 * 1000.0 / (double)ti1);
 384 #endif
 385 }
 386
 387 static uint8_t img_dest[64] __attribute__ ((aligned (8)));
 388 static uint8_t img_dest1[64] __attribute__ ((aligned (8)));
 389
 390 static void idct248_ref(uint8_t *dest, int linesize, int16_t *block)
 391 {
 392     static int init;
 393     static double c8[8][8];
 394     static double c4[4][4];
 395     double block1[64], block2[64], block3[64];
 396     double s, sum, v;
 397     int i, j, k;
 398
 399     if (!init) {
 400         init = 1;
 401
 402         for(i=0;i<8;i++) {
 403             sum = 0;
 404             for(j=0;j<8;j++) {
 405                 s = (i==0) ? sqrt(1.0/8.0) : sqrt(1.0/4.0);
 406                 c8[i][j] = s * cos(M_PI * i * (j + 0.5) / 8.0);
 407                 sum += c8[i][j] * c8[i][j];
 408             }
 409         }
 410
 411         for(i=0;i<4;i++) {
 412             sum = 0;
 413             for(j=0;j<4;j++) {
 414                 s = (i==0) ? sqrt(1.0/4.0) : sqrt(1.0/2.0);
 415                 c4[i][j] = s * cos(M_PI * i * (j + 0.5) / 4.0);
 416                 sum += c4[i][j] * c4[i][j];
 417             }
 418         }
 419     }
 420
 421     /* butterfly */
 422     s = 0.5 * sqrt(2.0);
 423     for(i=0;i<4;i++) {
 424         for(j=0;j<8;j++) {
 425             block1[8*(2*i)+j] = (block[8*(2*i)+j] + block[8*(2*i+1)+j]) * s;
 426             block1[8*(2*i+1)+j] = (block[8*(2*i)+j] - block[8*(2*i+1)+j]) * s;
 427         }
 428     }
 429
 430     /* idct8 on lines */
 431     for(i=0;i<8;i++) {
 432         for(j=0;j<8;j++) {
 433             sum = 0;
 434             for(k=0;k<8;k++)
 435                 sum += c8[k][j] * block1[8*i+k];
 436             block2[8*i+j] = sum;
 437         }
 438     }
 439
 440     /* idct4 */
 441     for(i=0;i<8;i++) {
 442         for(j=0;j<4;j++) {
 443             /* top */
 444             sum = 0;
 445             for(k=0;k<4;k++)
 446                 sum += c4[k][j] * block2[8*(2*k)+i];
 447             block3[8*(2*j)+i] = sum;
 448
 449             /* bottom */
 450             sum = 0;
 451             for(k=0;k<4;k++)
 452                 sum += c4[k][j] * block2[8*(2*k+1)+i];
 453             block3[8*(2*j+1)+i] = sum;
 454         }
 455     }
 456
 457     /* clamp and store the result */
 458     for(i=0;i<8;i++) {
 459         for(j=0;j<8;j++) {
 460             v = block3[8*i+j];
 461             if (v < 0)
 462                 v = 0;
 463             else if (v > 255)
 464                 v = 255;
 465             dest[i * linesize + j] = (int)rint(v);
 466         }
 467     }
 468 }
 469
 470 static void idct248_error(const char *name,
 471                     void (*idct248_put)(uint8_t *dest, int line_size, int16_t *block))
 472 {
 473     int it, i, it1, ti, ti1, err_max, v;
 474
 475     AVLFG prng;
 476
 477     av_lfg_init(&prng, 1);
 478
 479     /* just one test to see if code is correct (precision is less
 480        important here) */
 481     err_max = 0;
 482     for(it=0;it<NB_ITS;it++) {
 483
 484         /* XXX: use forward transform to generate values */
 485         for(i=0;i<64;i++)
 486             block1[i] = av_lfg_get(&prng) % 256 - 128;
 487         block1[0] += 1024;
 488
 489         for(i=0; i<64; i++)
 490             block[i]= block1[i];
 491         idct248_ref(img_dest1, 8, block);
 492
 493         for(i=0; i<64; i++)
 494             block[i]= block1[i];
 495         idct248_put(img_dest, 8, block);
 496
 497         for(i=0;i<64;i++) {
 498             v = abs((int)img_dest[i] - (int)img_dest1[i]);
 499             if (v == 255)
 500                 printf("%d %d\n", img_dest[i], img_dest1[i]);
 501             if (v > err_max)
 502                 err_max = v;
 503         }
 504 #if 0
 505         printf("ref=\n");
 506         for(i=0;i<8;i++) {
 507             int j;
 508             for(j=0;j<8;j++) {
 509                 printf(" %3d", img_dest1[i*8+j]);
 510             }
 511             printf("\n");
 512         }
 513
 514         printf("out=\n");
 515         for(i=0;i<8;i++) {
 516             int j;
 517             for(j=0;j<8;j++) {
 518                 printf(" %3d", img_dest[i*8+j]);
 519             }
 520             printf("\n");
 521         }
 522 #endif
 523     }
 524     printf("%s %s: err_inf=%d\n",
 525            1 ? "IDCT248" : "DCT248",
 526            name, err_max);
 527
 528     ti = gettime();
 529     it1 = 0;
 530     do {
 531         for(it=0;it<NB_ITS_SPEED;it++) {
 532             for(i=0; i<64; i++)
 533                 block[i]= block1[i];
 534 //            memcpy(block, block1, sizeof(DCTELEM) * 64);
 535 // do not memcpy especially not fastmemcpy because it does movntq !!!
 536             idct248_put(img_dest, 8, block);
 537         }
 538         it1 += NB_ITS_SPEED;
 539         ti1 = gettime() - ti;
 540     } while (ti1 < 1000000);
 541     mmx_emms();
 542
 543     printf("%s %s: %0.1f kdct/s\n",
 544            1 ? "IDCT248" : "DCT248",
 545            name, (double)it1 * 1000.0 / (double)ti1);
 546 }
 547
 548 static void help(void)
 549 {
 550     printf("dct-test [-i] [<test-number>]\n"
 551            "test-number 0 -> test with random matrixes\n"
 552            "            1 -> test with random sparse matrixes\n"
 553            "            2 -> do 3. test from mpeg4 std\n"
 554            "-i          test IDCT implementations\n"
 555            "-4          test IDCT248 implementations\n");
 556 }
 557
 558 int main(int argc, char **argv)
 559 {
 560     int test_idct = 0, test_248_dct = 0;
 561     int c,i;
 562     int test=1;
 563     cpu_flags = mm_support();
 564
 565     ff_ref_dct_init();
 566     idct_mmx_init();
 567
 568     for(i=0;i<256;i++) cropTbl[i + MAX_NEG_CROP] = i;
 569     for(i=0;i<MAX_NEG_CROP;i++) {
 570         cropTbl[i] = 0;
 571         cropTbl[i + MAX_NEG_CROP + 256] = 255;
 572     }
 573
 574     for(;;) {
 575         c = getopt(argc, argv, "ih4");
 576         if (c == -1)
 577             break;
 578         switch(c) {
 579         case 'i':
 580             test_idct = 1;
 581             break;
 582         case '4':
 583             test_248_dct = 1;
 584             break;
 585         default :
 586         case 'h':
 587             help();
 588             return 0;
 589         }
 590     }
 591
 592     if(optind <argc) test= atoi(argv[optind]);
 593
 594     printf("ffmpeg DCT/IDCT test\n");
 595
 596     if (test_248_dct) {
 597         idct248_error("SIMPLE-C", ff_simple_idct248_put);
 598     } else {
 599       for (i=0;algos[i].name;i++)
 600         if (algos[i].is_idct == test_idct && !(~cpu_flags & algos[i].mm_support)) {
 601           dct_error (algos[i].name, algos[i].is_idct, algos[i].func, algos[i].ref, algos[i].format, test);
 602         }
 603     }
 604     return 0;
 605 }