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

Windows CE

开发平台:

C/C++

  1. /**
  2.  * @file fdctref.c
  3.  * forward discrete cosine transform, double precision.
  4.  */
  5. /* Copyright (C) 1996, MPEG Software Simulation Group. All Rights Reserved. */
  6. /*
  7.  * Disclaimer of Warranty
  8.  *
  9.  * These software programs are available to the user without any license fee or
  10.  * royalty on an "as is" basis.  The MPEG Software Simulation Group disclaims
  11.  * any and all warranties, whether express, implied, or statuary, including any
  12.  * implied warranties or merchantability or of fitness for a particular
  13.  * purpose.  In no event shall the copyright-holder be liable for any
  14.  * incidental, punitive, or consequential damages of any kind whatsoever
  15.  * arising from the use of these programs.
  16.  *
  17.  * This disclaimer of warranty extends to the user of these programs and user's
  18.  * customers, employees, agents, transferees, successors, and assigns.
  19.  *
  20.  * The MPEG Software Simulation Group does not represent or warrant that the
  21.  * programs furnished hereunder are free of infringement of any third-party
  22.  * patents.
  23.  *
  24.  * Commercial implementations of MPEG-1 and MPEG-2 video, including shareware,
  25.  * are subject to royalty fees to patent holders.  Many of these patents are
  26.  * general enough such that they are unavoidable regardless of implementation
  27.  * design.
  28.  *
  29.  */
  30. #include <math.h>
  31. #ifndef PI
  32. # ifdef M_PI
  33. #  define PI M_PI
  34. # else
  35. #  define PI 3.14159265358979323846
  36. # endif
  37. #endif
  38. /* global declarations */
  39. void init_fdct (void);
  40. void fdct (short *block);
  41. /* private data */
  42. static double c[8][8]; /* transform coefficients */
  43. void init_fdct()
  44. {
  45.   int i, j;
  46.   double s;
  47.   for (i=0; i<8; i++)
  48.   {
  49.     s = (i==0) ? sqrt(0.125) : 0.5;
  50.     for (j=0; j<8; j++)
  51.       c[i][j] = s * cos((PI/8.0)*i*(j+0.5));
  52.   }
  53. }
  54. void fdct(block)
  55. short *block;
  56. {
  57. register int i, j;
  58. double s;
  59. double tmp[64];
  60. for(i = 0; i < 8; i++)
  61.      for(j = 0; j < 8; j++)
  62.      {
  63.      s = 0.0;
  64. /*
  65.  *      for(k = 0; k < 8; k++)
  66.  *          s += c[j][k] * block[8 * i + k];
  67.  */
  68.          s += c[j][0] * block[8 * i + 0];
  69.          s += c[j][1] * block[8 * i + 1];
  70.          s += c[j][2] * block[8 * i + 2];
  71.          s += c[j][3] * block[8 * i + 3];
  72.          s += c[j][4] * block[8 * i + 4];
  73.          s += c[j][5] * block[8 * i + 5];
  74.          s += c[j][6] * block[8 * i + 6];
  75.          s += c[j][7] * block[8 * i + 7];
  76.      tmp[8 * i + j] = s;
  77.      }
  78. for(j = 0; j < 8; j++)
  79.      for(i = 0; i < 8; i++)
  80.      {
  81.      s = 0.0;
  82. /*
  83.  *         for(k = 0; k < 8; k++)
  84.  *             s += c[i][k] * tmp[8 * k + j];
  85.  */
  86.          s += c[i][0] * tmp[8 * 0 + j];
  87.          s += c[i][1] * tmp[8 * 1 + j];
  88.          s += c[i][2] * tmp[8 * 2 + j];
  89.          s += c[i][3] * tmp[8 * 3 + j];
  90.          s += c[i][4] * tmp[8 * 4 + j];
  91.          s += c[i][5] * tmp[8 * 5 + j];
  92.          s += c[i][6] * tmp[8 * 6 + j];
  93.          s += c[i][7] * tmp[8 * 7 + j];
  94. s*=8.0;
  95.      block[8 * i + j] = (short)floor(s + 0.499999);
  96. /*
  97.  * reason for adding 0.499999 instead of 0.5:
  98.  * s is quite often x.5 (at least for i and/or j = 0 or 4)
  99.  * and setting the rounding threshold exactly to 0.5 leads to an
  100.  * extremely high arithmetic implementation dependency of the result;
  101.  * s being between x.5 and x.500001 (which is now incorrectly rounded
  102.  * downwards instead of upwards) is assumed to occur less often
  103.  * (if at all)
  104.  */
  105.       }
  106. }
  107. /* perform IDCT matrix multiply for 8x8 coefficient block */
  108. void idct(block)
  109. short *block;
  110. {
  111.   int i, j, k, v;
  112.   double partial_product;
  113.   double tmp[64];
  114.   for (i=0; i<8; i++)
  115.     for (j=0; j<8; j++)
  116.     {
  117.       partial_product = 0.0;
  118.       for (k=0; k<8; k++)
  119.         partial_product+= c[k][j]*block[8*i+k];
  120.       tmp[8*i+j] = partial_product;
  121.     }
  122.   /* Transpose operation is integrated into address mapping by switching 
  123.      loop order of i and j */
  124.   for (j=0; j<8; j++)
  125.     for (i=0; i<8; i++)
  126.     {
  127.       partial_product = 0.0;
  128.       for (k=0; k<8; k++)
  129.         partial_product+= c[k][i]*tmp[8*k+j];
  130.       v = (int) floor(partial_product+0.5);
  131.       block[8*i+j] = v;
  132.     }
  133. }