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