dct-test.c
上传用户:wstnjxml
上传日期:2014-04-03
资源大小:7248k
文件大小:14k
源码类别:

Windows CE

开发平台:

C/C++

  1. /**
  2.  * @file dct-test.c
  3.  * DCT test. (c) 2001 Fabrice Bellard. 
  4.  * Started from sample code by Juan J. Sierralta P.
  5.  */
  6. #include <stdlib.h>
  7. #include <stdio.h>
  8. #include <string.h>
  9. #include <sys/time.h>
  10. #include <unistd.h>
  11. #include "dsputil.h"
  12. #include "i386/mmx.h"
  13. #include "simple_idct.h"
  14. #include "faandct.h"
  15. #ifndef MAX
  16. #define MAX(a, b)  (((a) > (b)) ? (a) : (b))
  17. #endif
  18. #undef printf
  19. void *fast_memcpy(void *a, const void *b, size_t c){return memcpy(a,b,c);};
  20. /* reference fdct/idct */
  21. extern void fdct(DCTELEM *block);
  22. extern void idct(DCTELEM *block);
  23. extern void ff_idct_xvid_mmx(DCTELEM *block);
  24. extern void ff_idct_xvid_mmx2(DCTELEM *block);
  25. extern void init_fdct();
  26. extern void j_rev_dct(DCTELEM *data);
  27. extern void ff_mmx_idct(DCTELEM *data);
  28. extern void ff_mmxext_idct(DCTELEM *data);
  29. extern void odivx_idct_c (short *block);
  30. #define AANSCALE_BITS 12
  31. static const unsigned short aanscales[64] = {
  32.     /* precomputed values scaled up by 14 bits */
  33.     16384, 22725, 21407, 19266, 16384, 12873,  8867,  4520,
  34.     22725, 31521, 29692, 26722, 22725, 17855, 12299,  6270,
  35.     21407, 29692, 27969, 25172, 21407, 16819, 11585,  5906,
  36.     19266, 26722, 25172, 22654, 19266, 15137, 10426,  5315,
  37.     16384, 22725, 21407, 19266, 16384, 12873,  8867,  4520,
  38.     12873, 17855, 16819, 15137, 12873, 10114,  6967,  3552,
  39.     8867, 12299, 11585, 10426,  8867,  6967,  4799,  2446,
  40.     4520,  6270,  5906,  5315,  4520,  3552,  2446,  1247
  41. };
  42. uint8_t cropTbl[256 + 2 * MAX_NEG_CROP];
  43. int64_t gettime(void)
  44. {
  45.     struct timeval tv;
  46.     gettimeofday(&tv,NULL);
  47.     return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
  48. }
  49. #define NB_ITS 20000
  50. #define NB_ITS_SPEED 50000
  51. static short idct_mmx_perm[64];
  52. static short idct_simple_mmx_perm[64]={
  53. 0x00, 0x08, 0x04, 0x09, 0x01, 0x0C, 0x05, 0x0D, 
  54. 0x10, 0x18, 0x14, 0x19, 0x11, 0x1C, 0x15, 0x1D, 
  55. 0x20, 0x28, 0x24, 0x29, 0x21, 0x2C, 0x25, 0x2D, 
  56. 0x12, 0x1A, 0x16, 0x1B, 0x13, 0x1E, 0x17, 0x1F, 
  57. 0x02, 0x0A, 0x06, 0x0B, 0x03, 0x0E, 0x07, 0x0F, 
  58. 0x30, 0x38, 0x34, 0x39, 0x31, 0x3C, 0x35, 0x3D, 
  59. 0x22, 0x2A, 0x26, 0x2B, 0x23, 0x2E, 0x27, 0x2F, 
  60. 0x32, 0x3A, 0x36, 0x3B, 0x33, 0x3E, 0x37, 0x3F,
  61. };
  62. void idct_mmx_init(void)
  63. {
  64.     int i;
  65.     /* the mmx/mmxext idct uses a reordered input, so we patch scan tables */
  66.     for (i = 0; i < 64; i++) {
  67. idct_mmx_perm[i] = (i & 0x38) | ((i & 6) >> 1) | ((i & 1) << 2);
  68. // idct_simple_mmx_perm[i] = simple_block_permute_op(i);
  69.     }
  70. }
  71. static DCTELEM block[64] __attribute__ ((aligned (8)));
  72. static DCTELEM block1[64] __attribute__ ((aligned (8)));
  73. static DCTELEM block_org[64] __attribute__ ((aligned (8)));
  74. void dct_error(const char *name, int is_idct,
  75.                void (*fdct_func)(DCTELEM *block),
  76.                void (*fdct_ref)(DCTELEM *block), int test)
  77. {
  78.     int it, i, scale;
  79.     int err_inf, v;
  80.     int64_t err2, ti, ti1, it1;
  81.     int64_t sysErr[64], sysErrMax=0;
  82.     int maxout=0;
  83.     int blockSumErrMax=0, blockSumErr;
  84.     srandom(0);
  85.     err_inf = 0;
  86.     err2 = 0;
  87.     for(i=0; i<64; i++) sysErr[i]=0;
  88.     for(it=0;it<NB_ITS;it++) {
  89.         for(i=0;i<64;i++)
  90.             block1[i] = 0;
  91.         switch(test){
  92.         case 0: 
  93.             for(i=0;i<64;i++)
  94.                 block1[i] = (random() % 512) -256;
  95.             if (is_idct){
  96.                 fdct(block1);
  97.                 for(i=0;i<64;i++)
  98.                     block1[i]>>=3;
  99.             }
  100.         break;
  101.         case 1:{
  102.             int num= (random()%10)+1;
  103.             for(i=0;i<num;i++)
  104.                 block1[random()%64] = (random() % 512) -256;
  105.         }break;
  106.         case 2:
  107.             block1[0]= (random()%4096)-2048;
  108.             block1[63]= (block1[0]&1)^1;
  109.         break;
  110.         }
  111. #if 0 // simulate mismatch control
  112. { int sum=0;
  113.         for(i=0;i<64;i++)
  114.            sum+=block1[i];
  115.         if((sum&1)==0) block1[63]^=1; 
  116. }
  117. #endif
  118.         for(i=0; i<64; i++)
  119.             block_org[i]= block1[i];
  120.         if (fdct_func == ff_mmx_idct ||
  121.             fdct_func == j_rev_dct || fdct_func == ff_mmxext_idct) {
  122.             for(i=0;i<64;i++)
  123.                 block[idct_mmx_perm[i]] = block1[i];
  124.         } else if(fdct_func == ff_simple_idct_mmx ) {
  125.             for(i=0;i<64;i++)
  126.                 block[idct_simple_mmx_perm[i]] = block1[i];
  127. } else {
  128.             for(i=0; i<64; i++)
  129.                 block[i]= block1[i];
  130.         }
  131. #if 0 // simulate mismatch control for tested IDCT but not the ref
  132. { int sum=0;
  133.         for(i=0;i<64;i++)
  134.            sum+=block[i];
  135.         if((sum&1)==0) block[63]^=1; 
  136. }
  137. #endif
  138.         fdct_func(block);
  139.         emms(); /* for ff_mmx_idct */
  140.         if (fdct_func == fdct_ifast 
  141. #ifndef FAAN_POSTSCALE        
  142.             || fdct_func == ff_faandct
  143. #endif
  144.             ) {
  145.             for(i=0; i<64; i++) {
  146.                 scale = 8*(1 << (AANSCALE_BITS + 11)) / aanscales[i];
  147.                 block[i] = (block[i] * scale /*+ (1<<(AANSCALE_BITS-1))*/) >> AANSCALE_BITS;
  148.             }
  149.         }
  150.         fdct_ref(block1);
  151.         blockSumErr=0;
  152.         for(i=0;i<64;i++) {
  153.             v = abs(block[i] - block1[i]);
  154.             if (v > err_inf)
  155.                 err_inf = v;
  156.             err2 += v * v;
  157.     sysErr[i] += block[i] - block1[i];
  158.     blockSumErr += v;
  159.     if( abs(block[i])>maxout) maxout=abs(block[i]);
  160.         }
  161.         if(blockSumErrMax < blockSumErr) blockSumErrMax= blockSumErr;
  162. #if 0 // print different matrix pairs
  163.         if(blockSumErr){
  164.             printf("n");
  165.             for(i=0; i<64; i++){
  166.                 if((i&7)==0) printf("n");
  167.                 printf("%4d ", block_org[i]);
  168.             }
  169.             for(i=0; i<64; i++){
  170.                 if((i&7)==0) printf("n");
  171.                 printf("%4d ", block[i] - block1[i]);
  172.             }
  173.         }
  174. #endif
  175.     }
  176.     for(i=0; i<64; i++) sysErrMax= MAX(sysErrMax, ABS(sysErr[i]));
  177.     
  178. #if 1 // dump systematic errors
  179.     for(i=0; i<64; i++){
  180. if(i%8==0) printf("n");
  181.         printf("%5d ", (int)sysErr[i]);
  182.     }
  183.     printf("n");
  184. #endif
  185.     
  186.     printf("%s %s: err_inf=%d err2=%0.8f syserr=%0.8f maxout=%d blockSumErr=%dn",
  187.            is_idct ? "IDCT" : "DCT",
  188.            name, err_inf, (double)err2 / NB_ITS / 64.0, (double)sysErrMax / NB_ITS, maxout, blockSumErrMax);
  189. #if 1 //Speed test
  190.     /* speed test */
  191.     for(i=0;i<64;i++)
  192.         block1[i] = 0;
  193.     switch(test){
  194.     case 0: 
  195.         for(i=0;i<64;i++)
  196.             block1[i] = (random() % 512) -256;
  197.         if (is_idct){
  198.             fdct(block1);
  199.             for(i=0;i<64;i++)
  200.                 block1[i]>>=3;
  201.         }
  202.     break;
  203.     case 1:{
  204.     case 2:
  205.         block1[0] = (random() % 512) -256;
  206.         block1[1] = (random() % 512) -256;
  207.         block1[2] = (random() % 512) -256;
  208.         block1[3] = (random() % 512) -256;
  209.     }break;
  210.     }
  211.     if (fdct_func == ff_mmx_idct ||
  212.         fdct_func == j_rev_dct || fdct_func == ff_mmxext_idct) {
  213.         for(i=0;i<64;i++)
  214.             block[idct_mmx_perm[i]] = block1[i];
  215.     } else if(fdct_func == ff_simple_idct_mmx ) {
  216.         for(i=0;i<64;i++)
  217.             block[idct_simple_mmx_perm[i]] = block1[i];
  218.     } else {
  219.         for(i=0; i<64; i++)
  220.             block[i]= block1[i];
  221.     }
  222.     ti = gettime();
  223.     it1 = 0;
  224.     do {
  225.         for(it=0;it<NB_ITS_SPEED;it++) {
  226.             for(i=0; i<64; i++)
  227.                 block[i]= block1[i];
  228. //            memcpy(block, block1, sizeof(DCTELEM) * 64);
  229. // dont memcpy especially not fastmemcpy because it does movntq !!!
  230.             fdct_func(block);
  231.         }
  232.         it1 += NB_ITS_SPEED;
  233.         ti1 = gettime() - ti;
  234.     } while (ti1 < 1000000);
  235.     emms();
  236.     printf("%s %s: %0.1f kdct/sn",
  237.            is_idct ? "IDCT" : "DCT",
  238.            name, (double)it1 * 1000.0 / (double)ti1);
  239. #endif
  240. }
  241. static uint8_t img_dest[64] __attribute__ ((aligned (8)));
  242. static uint8_t img_dest1[64] __attribute__ ((aligned (8)));
  243. void idct248_ref(uint8_t *dest, int linesize, int16_t *block)
  244. {
  245.     static int init;
  246.     static double c8[8][8];
  247.     static double c4[4][4];
  248.     double block1[64], block2[64], block3[64];
  249.     double s, sum, v;
  250.     int i, j, k;
  251.     if (!init) {
  252.         init = 1;
  253.         for(i=0;i<8;i++) {
  254.             sum = 0;
  255.             for(j=0;j<8;j++) {
  256.                 s = (i==0) ? sqrt(1.0/8.0) : sqrt(1.0/4.0);
  257.                 c8[i][j] = s * cos(M_PI * i * (j + 0.5) / 8.0);
  258.                 sum += c8[i][j] * c8[i][j];
  259.             }
  260.         }
  261.         
  262.         for(i=0;i<4;i++) {
  263.             sum = 0;
  264.             for(j=0;j<4;j++) {
  265.                 s = (i==0) ? sqrt(1.0/4.0) : sqrt(1.0/2.0);
  266.                 c4[i][j] = s * cos(M_PI * i * (j + 0.5) / 4.0);
  267.                 sum += c4[i][j] * c4[i][j];
  268.             }
  269.         }
  270.     }
  271.     /* butterfly */
  272.     s = 0.5 * sqrt(2.0);
  273.     for(i=0;i<4;i++) {
  274.         for(j=0;j<8;j++) {
  275.             block1[8*(2*i)+j] = (block[8*(2*i)+j] + block[8*(2*i+1)+j]) * s;
  276.             block1[8*(2*i+1)+j] = (block[8*(2*i)+j] - block[8*(2*i+1)+j]) * s;
  277.         }
  278.     }
  279.     /* idct8 on lines */
  280.     for(i=0;i<8;i++) {
  281.         for(j=0;j<8;j++) {
  282.             sum = 0;
  283.             for(k=0;k<8;k++)
  284.                 sum += c8[k][j] * block1[8*i+k];
  285.             block2[8*i+j] = sum;
  286.         }
  287.     }
  288.     /* idct4 */
  289.     for(i=0;i<8;i++) {
  290.         for(j=0;j<4;j++) {
  291.             /* top */
  292.             sum = 0;
  293.             for(k=0;k<4;k++)
  294.                 sum += c4[k][j] * block2[8*(2*k)+i];
  295.             block3[8*(2*j)+i] = sum;
  296.             /* bottom */
  297.             sum = 0;
  298.             for(k=0;k<4;k++)
  299.                 sum += c4[k][j] * block2[8*(2*k+1)+i];
  300.             block3[8*(2*j+1)+i] = sum;
  301.         }
  302.     }
  303.     /* clamp and store the result */
  304.     for(i=0;i<8;i++) {
  305.         for(j=0;j<8;j++) {
  306.             v = block3[8*i+j];
  307.             if (v < 0)
  308.                 v = 0;
  309.             else if (v > 255)
  310.                 v = 255;
  311.             dest[i * linesize + j] = (int)rint(v);
  312.         }
  313.     }
  314. }
  315. void idct248_error(const char *name, 
  316.                     void (*idct248_put)(uint8_t *dest, int line_size, int16_t *block))
  317. {
  318.     int it, i, it1, ti, ti1, err_max, v;
  319.     srandom(0);
  320.     
  321.     /* just one test to see if code is correct (precision is less
  322.        important here) */
  323.     err_max = 0;
  324.     for(it=0;it<NB_ITS;it++) {
  325.         
  326.         /* XXX: use forward transform to generate values */
  327.         for(i=0;i<64;i++)
  328.             block1[i] = (random() % 256) - 128;
  329.         block1[0] += 1024;
  330.         for(i=0; i<64; i++)
  331.             block[i]= block1[i];
  332.         idct248_ref(img_dest1, 8, block);
  333.         
  334.         for(i=0; i<64; i++)
  335.             block[i]= block1[i];
  336.         idct248_put(img_dest, 8, block);
  337.         
  338.         for(i=0;i<64;i++) {
  339.             v = abs((int)img_dest[i] - (int)img_dest1[i]);
  340.             if (v == 255)
  341.                 printf("%d %dn", img_dest[i], img_dest1[i]);
  342.             if (v > err_max)
  343.                 err_max = v;
  344.         }
  345. #if 0
  346.         printf("ref=n");
  347.         for(i=0;i<8;i++) {
  348.             int j;
  349.             for(j=0;j<8;j++) {
  350.                 printf(" %3d", img_dest1[i*8+j]);
  351.             }
  352.             printf("n");
  353.         }
  354.         
  355.         printf("out=n");
  356.         for(i=0;i<8;i++) {
  357.             int j;
  358.             for(j=0;j<8;j++) {
  359.                 printf(" %3d", img_dest[i*8+j]);
  360.             }
  361.             printf("n");
  362.         }
  363. #endif
  364.     }
  365.     printf("%s %s: err_inf=%dn",
  366.            1 ? "IDCT248" : "DCT248",
  367.            name, err_max);
  368.     ti = gettime();
  369.     it1 = 0;
  370.     do {
  371.         for(it=0;it<NB_ITS_SPEED;it++) {
  372.             for(i=0; i<64; i++)
  373.                 block[i]= block1[i];
  374. //            memcpy(block, block1, sizeof(DCTELEM) * 64);
  375. // dont memcpy especially not fastmemcpy because it does movntq !!!
  376.             idct248_put(img_dest, 8, block);
  377.         }
  378.         it1 += NB_ITS_SPEED;
  379.         ti1 = gettime() - ti;
  380.     } while (ti1 < 1000000);
  381.     emms();
  382.     printf("%s %s: %0.1f kdct/sn",
  383.            1 ? "IDCT248" : "DCT248",
  384.            name, (double)it1 * 1000.0 / (double)ti1);
  385. }
  386. void help(void)
  387. {
  388.     printf("dct-test [-i] [<test-number>]n"
  389.            "test-number 0 -> test with random matrixesn"
  390.            "            1 -> test with random sparse matrixesn"
  391.            "            2 -> do 3. test from mpeg4 stdn"
  392.            "-i          test IDCT implementationsn"
  393.            "-4          test IDCT248 implementationsn");
  394.     exit(1);
  395. }
  396. int main(int argc, char **argv)
  397. {
  398.     int test_idct = 0, test_248_dct = 0;
  399.     int c,i;
  400.     int test=1;
  401.     init_fdct();
  402.     idct_mmx_init();
  403.     for(i=0;i<256;i++) cropTbl[i + MAX_NEG_CROP] = i;
  404.     for(i=0;i<MAX_NEG_CROP;i++) {
  405.         cropTbl[i] = 0;
  406.         cropTbl[i + MAX_NEG_CROP + 256] = 255;
  407.     }
  408.     
  409.     for(;;) {
  410.         c = getopt(argc, argv, "ih4");
  411.         if (c == -1)
  412.             break;
  413.         switch(c) {
  414.         case 'i':
  415.             test_idct = 1;
  416.             break;
  417.         case '4':
  418.             test_248_dct = 1;
  419.             break;
  420.         default :
  421.         case 'h':
  422.             help();
  423.             break;
  424.         }
  425.     }
  426.     
  427.     if(optind <argc) test= atoi(argv[optind]);
  428.                
  429.     printf("ffmpeg DCT/IDCT testn");
  430.     if (test_248_dct) {
  431.         idct248_error("SIMPLE-C", simple_idct248_put);
  432.     } else {
  433.         if (!test_idct) {
  434.             dct_error("REF-DBL", 0, fdct, fdct, test); /* only to verify code ! */
  435.             dct_error("IJG-AAN-INT", 0, fdct_ifast, fdct, test);
  436.             dct_error("IJG-LLM-INT", 0, ff_jpeg_fdct_islow, fdct, test);
  437.             dct_error("MMX", 0, ff_fdct_mmx, fdct, test);
  438.             dct_error("MMX2", 0, ff_fdct_mmx2, fdct, test);
  439.             dct_error("FAAN", 0, ff_faandct, fdct, test);
  440.         } else {
  441.             dct_error("REF-DBL", 1, idct, idct, test);
  442.             dct_error("INT", 1, j_rev_dct, idct, test);
  443.             dct_error("LIBMPEG2-MMX", 1, ff_mmx_idct, idct, test);
  444.             dct_error("LIBMPEG2-MMXEXT", 1, ff_mmxext_idct, idct, test);
  445.             dct_error("SIMPLE-C", 1, simple_idct, idct, test);
  446.             dct_error("SIMPLE-MMX", 1, ff_simple_idct_mmx, idct, test);
  447.             dct_error("XVID-MMX", 1, ff_idct_xvid_mmx, idct, test);
  448.             dct_error("XVID-MMX2", 1, ff_idct_xvid_mmx2, idct, test);
  449.             //        dct_error("ODIVX-C", 1, odivx_idct_c, idct);
  450.             //printf(" test against odivx idctn");
  451.             // dct_error("REF", 1, idct, odivx_idct_c);
  452.             //        dct_error("INT", 1, j_rev_dct, odivx_idct_c);
  453.             //        dct_error("MMX", 1, ff_mmx_idct, odivx_idct_c);
  454.             //        dct_error("MMXEXT", 1, ff_mmxext_idct, odivx_idct_c);
  455.             //        dct_error("SIMPLE-C", 1, simple_idct, odivx_idct_c);
  456.             //        dct_error("SIMPLE-MMX", 1, ff_simple_idct_mmx, odivx_idct_c);
  457.             //        dct_error("ODIVX-C", 1, odivx_idct_c, odivx_idct_c);
  458.         }
  459.     }
  460.     return 0;
  461. }