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
  37 #include "simple_idct.h"
  38 #include "aandcttab.h"
  39 #include "faandct.h"
  40 #include "faanidct.h"
  41 #include "x86/idct_xvid.h"
  42
  43 #undef printf
  44 #undef random
  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 fdct(DCTELEM *block);
  50 void idct(DCTELEM *block);
  51 void init_fdct(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, fdct,               fdct, NO_PERM},
  94   {"FAAN",            0, ff_faandct,         fdct, FAAN_SCALE},
  95   {"FAANI",           1, ff_faanidct,        idct, NO_PERM},
  96   {"IJG-AAN-INT",     0, fdct_ifast,         fdct, SCALE_PERM},
  97   {"IJG-LLM-INT",     0, ff_jpeg_fdct_islow, fdct, NO_PERM},
  98   {"REF-DBL",         1, idct,               idct, NO_PERM},
  99   {"INT",             1, j_rev_dct,          idct, MMX_PERM},
 100   {"SIMPLE-C",        1, ff_simple_idct,     idct, NO_PERM},
 101
 102 #if HAVE_MMX
 103   {"MMX",             0, ff_fdct_mmx,        fdct, NO_PERM, FF_MM_MMX},
 104 #if HAVE_MMX2
 105   {"MMX2",            0, ff_fdct_mmx2,       fdct, NO_PERM, FF_MM_MMXEXT},
 106   {"SSE2",            0, ff_fdct_sse2,       fdct, NO_PERM, FF_MM_SSE2},
 107 #endif
 108
 109 #if CONFIG_GPL
 110   {"LIBMPEG2-MMX",    1, ff_mmx_idct,        idct, MMX_PERM, FF_MM_MMX},
 111   {"LIBMPEG2-MMXEXT", 1, ff_mmxext_idct,     idct, MMX_PERM, FF_MM_MMXEXT},
 112 #endif
 113   {"SIMPLE-MMX",      1, ff_simple_idct_mmx, idct, MMX_SIMPLE_PERM, FF_MM_MMX},
 114   {"XVID-MMX",        1, ff_idct_xvid_mmx,   idct, NO_PERM, FF_MM_MMX},
 115   {"XVID-MMX2",       1, ff_idct_xvid_mmx2,  idct, NO_PERM, FF_MM_MMXEXT},
 116   {"XVID-SSE2",       1, ff_idct_xvid_sse2,  idct, SSE2_PERM, FF_MM_SSE2},
 117 #endif
 118
 119 #if HAVE_ALTIVEC
 120   {"altivecfdct",     0, fdct_altivec,       fdct, NO_PERM, FF_MM_ALTIVEC},
 121 #endif
 122
 123 #if ARCH_BFIN
 124   {"BFINfdct",        0, ff_bfin_fdct,       fdct, NO_PERM},
 125   {"BFINidct",        1, ff_bfin_idct,       idct, NO_PERM},
 126 #endif
 127
 128 #if ARCH_ARM
 129   {"SIMPLE-ARM",      1, simple_idct_ARM,    idct, NO_PERM },
 130   {"INT-ARM",         1, j_rev_dct_ARM,      idct, MMX_PERM },
 131 #if HAVE_ARMV5TE
 132   {"SIMPLE-ARMV5TE",  1, simple_idct_armv5te, idct, NO_PERM },
 133 #endif
 134 #if HAVE_ARMV6
 135   {"SIMPLE-ARMV6",    1, ff_simple_idct_armv6, idct, MMX_PERM },
 136 #endif
 137 #if HAVE_NEON
 138   {"SIMPLE-NEON",     1, ff_simple_idct_neon, idct, PARTTRANS_PERM },
 139 #endif
 140 #endif /* ARCH_ARM */
 141
 142 #if ARCH_ALPHA
 143   {"SIMPLE-ALPHA",    1, ff_simple_idct_axp,  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 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 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 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
 212     srandom(0);
 213
 214     err_inf = 0;
 215     err2 = 0;
 216     for(i=0; i<64; i++) sysErr[i]=0;
 217     for(it=0;it<NB_ITS;it++) {
 218         for(i=0;i<64;i++)
 219             block1[i] = 0;
 220         switch(test){
 221         case 0:
 222             for(i=0;i<64;i++)
 223                 block1[i] = (random() % 512) -256;
 224             if (is_idct){
 225                 fdct(block1);
 226
 227                 for(i=0;i<64;i++)
 228                     block1[i]>>=3;
 229             }
 230         break;
 231         case 1:{
 232             int num= (random()%10)+1;
 233             for(i=0;i<num;i++)
 234                 block1[random()%64] = (random() % 512) -256;
 235         }break;
 236         case 2:
 237             block1[0]= (random()%4096)-2048;
 238             block1[63]= (block1[0]&1)^1;
 239         break;
 240         }
 241
 242 #if 0 // simulate mismatch control
 243 { int sum=0;
 244         for(i=0;i<64;i++)
 245            sum+=block1[i];
 246
 247         if((sum&1)==0) block1[63]^=1;
 248 }
 249 #endif
 250
 251         for(i=0; i<64; i++)
 252             block_org[i]= block1[i];
 253
 254         if (form == MMX_PERM) {
 255             for(i=0;i<64;i++)
 256                 block[idct_mmx_perm[i]] = block1[i];
 257             } else if (form == MMX_SIMPLE_PERM) {
 258             for(i=0;i<64;i++)
 259                 block[idct_simple_mmx_perm[i]] = block1[i];
 260
 261         } else if (form == SSE2_PERM) {
 262             for(i=0; i<64; i++)
 263                 block[(i&0x38) | idct_sse2_row_perm[i&7]] = block1[i];
 264         } else if (form == PARTTRANS_PERM) {
 265             for(i=0; i<64; i++)
 266                 block[(i&0x24) | ((i&3)<<3) | ((i>>3)&3)] = block1[i];
 267         } else {
 268             for(i=0; i<64; i++)
 269                 block[i]= block1[i];
 270         }
 271 #if 0 // simulate mismatch control for tested IDCT but not the ref
 272 { int sum=0;
 273         for(i=0;i<64;i++)
 274            sum+=block[i];
 275
 276         if((sum&1)==0) block[63]^=1;
 277 }
 278 #endif
 279
 280         fdct_func(block);
 281         mmx_emms();
 282
 283         if (form == SCALE_PERM) {
 284             for(i=0; i<64; i++) {
 285                 scale = 8*(1 << (AANSCALE_BITS + 11)) / ff_aanscales[i];
 286                 block[i] = (block[i] * scale /*+ (1<<(AANSCALE_BITS-1))*/) >> AANSCALE_BITS;
 287             }
 288         }
 289
 290         fdct_ref(block1);
 291
 292         blockSumErr=0;
 293         for(i=0;i<64;i++) {
 294             v = abs(block[i] - block1[i]);
 295             if (v > err_inf)
 296                 err_inf = v;
 297             err2 += v * v;
 298             sysErr[i] += block[i] - block1[i];
 299             blockSumErr += v;
 300             if( abs(block[i])>maxout) maxout=abs(block[i]);
 301         }
 302         if(blockSumErrMax < blockSumErr) blockSumErrMax= blockSumErr;
 303 #if 0 // print different matrix pairs
 304         if(blockSumErr){
 305             printf("\n");
 306             for(i=0; i<64; i++){
 307                 if((i&7)==0) printf("\n");
 308                 printf("%4d ", block_org[i]);
 309             }
 310             for(i=0; i<64; i++){
 311                 if((i&7)==0) printf("\n");
 312                 printf("%4d ", block[i] - block1[i]);
 313             }
 314         }
 315 #endif
 316     }
 317     for(i=0; i<64; i++) sysErrMax= FFMAX(sysErrMax, FFABS(sysErr[i]));
 318
 319 #if 1 // dump systematic errors
 320     for(i=0; i<64; i++){
 321         if(i%8==0) printf("\n");
 322         printf("%5d ", (int)sysErr[i]);
 323     }
 324     printf("\n");
 325 #endif
 326
 327     printf("%s %s: err_inf=%d err2=%0.8f syserr=%0.8f maxout=%d blockSumErr=%d\n",
 328            is_idct ? "IDCT" : "DCT",
 329            name, err_inf, (double)err2 / NB_ITS / 64.0, (double)sysErrMax / NB_ITS, maxout, blockSumErrMax);
 330 #if 1 //Speed test
 331     /* speed test */
 332     for(i=0;i<64;i++)
 333         block1[i] = 0;
 334     switch(test){
 335     case 0:
 336         for(i=0;i<64;i++)
 337             block1[i] = (random() % 512) -256;
 338         if (is_idct){
 339             fdct(block1);
 340
 341             for(i=0;i<64;i++)
 342                 block1[i]>>=3;
 343         }
 344     break;
 345     case 1:{
 346     case 2:
 347         block1[0] = (random() % 512) -256;
 348         block1[1] = (random() % 512) -256;
 349         block1[2] = (random() % 512) -256;
 350         block1[3] = (random() % 512) -256;
 351     }break;
 352     }
 353
 354     if (form == MMX_PERM) {
 355         for(i=0;i<64;i++)
 356             block[idct_mmx_perm[i]] = block1[i];
 357     } else if(form == MMX_SIMPLE_PERM) {
 358         for(i=0;i<64;i++)
 359             block[idct_simple_mmx_perm[i]] = block1[i];
 360     } else {
 361         for(i=0; i<64; i++)
 362             block[i]= block1[i];
 363     }
 364
 365     ti = gettime();
 366     it1 = 0;
 367     do {
 368         for(it=0;it<NB_ITS_SPEED;it++) {
 369             for(i=0; i<64; i++)
 370                 block[i]= block1[i];
 371 //            memcpy(block, block1, sizeof(DCTELEM) * 64);
 372 // do not memcpy especially not fastmemcpy because it does movntq !!!
 373             fdct_func(block);
 374         }
 375         it1 += NB_ITS_SPEED;
 376         ti1 = gettime() - ti;
 377     } while (ti1 < 1000000);
 378     mmx_emms();
 379
 380     printf("%s %s: %0.1f kdct/s\n",
 381            is_idct ? "IDCT" : "DCT",
 382            name, (double)it1 * 1000.0 / (double)ti1);
 383 #endif
 384 }
 385
 386 static uint8_t img_dest[64] __attribute__ ((aligned (8)));
 387 static uint8_t img_dest1[64] __attribute__ ((aligned (8)));
 388
 389 void idct248_ref(uint8_t *dest, int linesize, int16_t *block)
 390 {
 391     static int init;
 392     static double c8[8][8];
 393     static double c4[4][4];
 394     double block1[64], block2[64], block3[64];
 395     double s, sum, v;
 396     int i, j, k;
 397
 398     if (!init) {
 399         init = 1;
 400
 401         for(i=0;i<8;i++) {
 402             sum = 0;
 403             for(j=0;j<8;j++) {
 404                 s = (i==0) ? sqrt(1.0/8.0) : sqrt(1.0/4.0);
 405                 c8[i][j] = s * cos(M_PI * i * (j + 0.5) / 8.0);
 406                 sum += c8[i][j] * c8[i][j];
 407             }
 408         }
 409
 410         for(i=0;i<4;i++) {
 411             sum = 0;
 412             for(j=0;j<4;j++) {
 413                 s = (i==0) ? sqrt(1.0/4.0) : sqrt(1.0/2.0);
 414                 c4[i][j] = s * cos(M_PI * i * (j + 0.5) / 4.0);
 415                 sum += c4[i][j] * c4[i][j];
 416             }
 417         }
 418     }
 419
 420     /* butterfly */
 421     s = 0.5 * sqrt(2.0);
 422     for(i=0;i<4;i++) {
 423         for(j=0;j<8;j++) {
 424             block1[8*(2*i)+j] = (block[8*(2*i)+j] + block[8*(2*i+1)+j]) * s;
 425             block1[8*(2*i+1)+j] = (block[8*(2*i)+j] - block[8*(2*i+1)+j]) * s;
 426         }
 427     }
 428
 429     /* idct8 on lines */
 430     for(i=0;i<8;i++) {
 431         for(j=0;j<8;j++) {
 432             sum = 0;
 433             for(k=0;k<8;k++)
 434                 sum += c8[k][j] * block1[8*i+k];
 435             block2[8*i+j] = sum;
 436         }
 437     }
 438
 439     /* idct4 */
 440     for(i=0;i<8;i++) {
 441         for(j=0;j<4;j++) {
 442             /* top */
 443             sum = 0;
 444             for(k=0;k<4;k++)
 445                 sum += c4[k][j] * block2[8*(2*k)+i];
 446             block3[8*(2*j)+i] = sum;
 447
 448             /* bottom */
 449             sum = 0;
 450             for(k=0;k<4;k++)
 451                 sum += c4[k][j] * block2[8*(2*k+1)+i];
 452             block3[8*(2*j+1)+i] = sum;
 453         }
 454     }
 455
 456     /* clamp and store the result */
 457     for(i=0;i<8;i++) {
 458         for(j=0;j<8;j++) {
 459             v = block3[8*i+j];
 460             if (v < 0)
 461                 v = 0;
 462             else if (v > 255)
 463                 v = 255;
 464             dest[i * linesize + j] = (int)rint(v);
 465         }
 466     }
 467 }
 468
 469 void idct248_error(const char *name,
 470                     void (*idct248_put)(uint8_t *dest, int line_size, int16_t *block))
 471 {
 472     int it, i, it1, ti, ti1, err_max, v;
 473
 474     srandom(0);
 475
 476     /* just one test to see if code is correct (precision is less
 477        important here) */
 478     err_max = 0;
 479     for(it=0;it<NB_ITS;it++) {
 480
 481         /* XXX: use forward transform to generate values */
 482         for(i=0;i<64;i++)
 483             block1[i] = (random() % 256) - 128;
 484         block1[0] += 1024;
 485
 486         for(i=0; i<64; i++)
 487             block[i]= block1[i];
 488         idct248_ref(img_dest1, 8, block);
 489
 490         for(i=0; i<64; i++)
 491             block[i]= block1[i];
 492         idct248_put(img_dest, 8, block);
 493
 494         for(i=0;i<64;i++) {
 495             v = abs((int)img_dest[i] - (int)img_dest1[i]);
 496             if (v == 255)
 497                 printf("%d %d\n", img_dest[i], img_dest1[i]);
 498             if (v > err_max)
 499                 err_max = v;
 500         }
 501 #if 0
 502         printf("ref=\n");
 503         for(i=0;i<8;i++) {
 504             int j;
 505             for(j=0;j<8;j++) {
 506                 printf(" %3d", img_dest1[i*8+j]);
 507             }
 508             printf("\n");
 509         }
 510
 511         printf("out=\n");
 512         for(i=0;i<8;i++) {
 513             int j;
 514             for(j=0;j<8;j++) {
 515                 printf(" %3d", img_dest[i*8+j]);
 516             }
 517             printf("\n");
 518         }
 519 #endif
 520     }
 521     printf("%s %s: err_inf=%d\n",
 522            1 ? "IDCT248" : "DCT248",
 523            name, err_max);
 524
 525     ti = gettime();
 526     it1 = 0;
 527     do {
 528         for(it=0;it<NB_ITS_SPEED;it++) {
 529             for(i=0; i<64; i++)
 530                 block[i]= block1[i];
 531 //            memcpy(block, block1, sizeof(DCTELEM) * 64);
 532 // do not memcpy especially not fastmemcpy because it does movntq !!!
 533             idct248_put(img_dest, 8, block);
 534         }
 535         it1 += NB_ITS_SPEED;
 536         ti1 = gettime() - ti;
 537     } while (ti1 < 1000000);
 538     mmx_emms();
 539
 540     printf("%s %s: %0.1f kdct/s\n",
 541            1 ? "IDCT248" : "DCT248",
 542            name, (double)it1 * 1000.0 / (double)ti1);
 543 }
 544
 545 void help(void)
 546 {
 547     printf("dct-test [-i] [<test-number>]\n"
 548            "test-number 0 -> test with random matrixes\n"
 549            "            1 -> test with random sparse matrixes\n"
 550            "            2 -> do 3. test from mpeg4 std\n"
 551            "-i          test IDCT implementations\n"
 552            "-4          test IDCT248 implementations\n");
 553 }
 554
 555 int main(int argc, char **argv)
 556 {
 557     int test_idct = 0, test_248_dct = 0;
 558     int c,i;
 559     int test=1;
 560     cpu_flags = mm_support();
 561
 562     init_fdct();
 563     idct_mmx_init();
 564
 565     for(i=0;i<256;i++) cropTbl[i + MAX_NEG_CROP] = i;
 566     for(i=0;i<MAX_NEG_CROP;i++) {
 567         cropTbl[i] = 0;
 568         cropTbl[i + MAX_NEG_CROP + 256] = 255;
 569     }
 570
 571     for(;;) {
 572         c = getopt(argc, argv, "ih4");
 573         if (c == -1)
 574             break;
 575         switch(c) {
 576         case 'i':
 577             test_idct = 1;
 578             break;
 579         case '4':
 580             test_248_dct = 1;
 581             break;
 582         default :
 583         case 'h':
 584             help();
 585             return 0;
 586         }
 587     }
 588
 589     if(optind <argc) test= atoi(argv[optind]);
 590
 591     printf("ffmpeg DCT/IDCT test\n");
 592
 593     if (test_248_dct) {
 594         idct248_error("SIMPLE-C", ff_simple_idct248_put);
 595     } else {
 596       for (i=0;algos[i].name;i++)
 597         if (algos[i].is_idct == test_idct && !(~cpu_flags & algos[i].mm_support)) {
 598           dct_error (algos[i].name, algos[i].is_idct, algos[i].func, algos[i].ref, algos[i].format, test);
 599         }
 600     }
 601     return 0;
 602 }