home *** CD-ROM | disk | FTP | other *** search
/ MPEG Toolkit / MPEG Toolkit.iso / os2 / mpegenc / src / mfwddct.c < prev    next >
Encoding:
C/C++ Source or Header  |  1997-01-01  |  8.1 KB  |  260 lines
  1.  
  2. /*
  3.  * mfwddct.c (derived from jfwddct.c, which carries the following info)
  4.  *
  5.  * Copyright (C) 1991, 1992, Thomas G. Lane. This file is part of the
  6.  * Independent JPEG Group's software. For conditions of distribution and use,
  7.  * see the accompanying README file.
  8.  *
  9.  * This file contains the basic DCT (Discrete Cosine Transform) transformation
  10.  * subroutine.
  11.  *
  12.  * This implementation is based on Appendix A.2 of the book "Discrete Cosine
  13.  * Transform---Algorithms, Advantages, Applications" by K.R. Rao and P. Yip
  14.  * (Academic Press, Inc, London, 1990). It uses scaled fixed-point arithmetic
  15.  * instead of floating point.
  16.  */
  17.  
  18. #include "all.h"
  19.  
  20. #include "dct.h"
  21. #include "mtypes.h"
  22.  
  23. /*
  24.  * The poop on this scaling stuff is as follows:
  25.  *
  26.  * We have to do addition and subtraction of the integer inputs, which is no
  27.  * problem, and multiplication by fractional constants, which is a problem to
  28.  * do in integer arithmetic.  We multiply all the constants by DCT_SCALE and
  29.  * convert them to integer constants (thus retaining LG2_DCT_SCALE bits of
  30.  * precision in the constants).  After doing a multiplication we have to
  31.  * divide the product by DCT_SCALE, with proper rounding, to produce the
  32.  * correct output.  The division can be implemented cheaply as a right shift
  33.  * of LG2_DCT_SCALE bits.  The DCT equations also specify an additional
  34.  * division by 2 on the final outputs; this can be folded into the
  35.  * right-shift by shifting one more bit (see UNFIXH).
  36.  *
  37.  * If you are planning to recode this in assembler, you might want to set
  38.  * LG2_DCT_SCALE to 15.  This loses a bit of precision, but then all the
  39.  * multiplications are between 16-bit quantities (given 8-bit JSAMPLEs!) so
  40.  * you could use a signed 16x16=>32 bit multiply instruction instead of full
  41.  * 32x32 multiply.  Unfortunately there's no way to describe such a multiply
  42.  * portably in C, so we've gone for the extra bit of accuracy here.
  43.  */
  44.  
  45. #define EIGHT_BIT_SAMPLES
  46. #ifdef EIGHT_BIT_SAMPLES
  47. #define LG2_DCT_SCALE 16
  48. #else
  49. #define LG2_DCT_SCALE 15    /* lose a little precision to avoid overflow */
  50. #endif
  51.  
  52. #define ONE    ((int32) 1)
  53.  
  54. #define DCT_SCALE (ONE << LG2_DCT_SCALE)
  55.  
  56. /* In some places we shift the inputs left by a couple more bits, */
  57. /* so that they can be added to fractional results without too much */
  58. /* loss of precision. */
  59. #define LG2_OVERSCALE 2
  60. #define OVERSCALE  (ONE << LG2_OVERSCALE)
  61. #define OVERSHIFT(x)  ((x) <<= LG2_OVERSCALE)
  62.  
  63. /* Scale a fractional constant by DCT_SCALE */
  64. #define FIX(x)    ((int32) ((x) * DCT_SCALE + 0.5))
  65.  
  66. /* Scale a fractional constant by DCT_SCALE/OVERSCALE */
  67. /* Such a constant can be multiplied with an overscaled input */
  68. /* to produce something that's scaled by DCT_SCALE */
  69. #define FIXO(x)  ((int32) ((x) * DCT_SCALE / OVERSCALE + 0.5))
  70.  
  71. /* Descale and correctly round a value that's scaled by DCT_SCALE */
  72. #define UNFIX(x)   RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1)), LG2_DCT_SCALE)
  73.  
  74. /* Same with an additional division by 2, ie, correctly rounded UNFIX(x/2) */
  75. #define UNFIXH(x)  RIGHT_SHIFT((x) + (ONE << LG2_DCT_SCALE), LG2_DCT_SCALE+1)
  76.  
  77. /* Take a value scaled by DCT_SCALE and round to integer scaled by OVERSCALE */
  78. #define UNFIXO(x)  RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1-LG2_OVERSCALE)),\
  79.                    LG2_DCT_SCALE-LG2_OVERSCALE)
  80.  
  81. /* Here are the constants we need */
  82. /* SIN_i_j is sine of i*pi/j, scaled by DCT_SCALE */
  83. /* COS_i_j is cosine of i*pi/j, scaled by DCT_SCALE */
  84.  
  85. #define SIN_1_4 FIX(0.707106781)
  86. #define COS_1_4 SIN_1_4
  87.  
  88. #define SIN_1_8 FIX(0.382683432)
  89. #define COS_1_8 FIX(0.923879533)
  90. #define SIN_3_8 COS_1_8
  91. #define COS_3_8 SIN_1_8
  92.  
  93. #define SIN_1_16 FIX(0.195090322)
  94. #define COS_1_16 FIX(0.980785280)
  95. #define SIN_7_16 COS_1_16
  96. #define COS_7_16 SIN_1_16
  97.  
  98. #define SIN_3_16 FIX(0.555570233)
  99. #define COS_3_16 FIX(0.831469612)
  100. #define SIN_5_16 COS_3_16
  101. #define COS_5_16 SIN_3_16
  102.  
  103. /* OSIN_i_j is sine of i*pi/j, scaled by DCT_SCALE/OVERSCALE */
  104. /* OCOS_i_j is cosine of i*pi/j, scaled by DCT_SCALE/OVERSCALE */
  105.  
  106. #define OSIN_1_4 FIXO(0.707106781)
  107. #define OCOS_1_4 OSIN_1_4
  108.  
  109. #define OSIN_1_8 FIXO(0.382683432)
  110. #define OCOS_1_8 FIXO(0.923879533)
  111. #define OSIN_3_8 OCOS_1_8
  112. #define OCOS_3_8 OSIN_1_8
  113.  
  114. #define OSIN_1_16 FIXO(0.195090322)
  115. #define OCOS_1_16 FIXO(0.980785280)
  116. #define OSIN_7_16 OCOS_1_16
  117. #define OCOS_7_16 OSIN_1_16
  118.  
  119. #define OSIN_3_16 FIXO(0.555570233)
  120. #define OCOS_3_16 FIXO(0.831469612)
  121. #define OSIN_5_16 OCOS_3_16
  122. #define OCOS_5_16 OSIN_3_16
  123.  
  124.  
  125.  
  126.  
  127. /*
  128.  * --------------------------------------------------------------
  129.  *
  130.  * mp_fwd_dct_block --
  131.  *
  132.  * Perform the forward DCT on one block of samples.
  133.  *
  134.  * A 2-D DCT can be done by 1-D DCT on each row followed by 1-D DCT on each
  135.  * column.
  136.  *
  137.  * Results: None
  138.  *
  139.  * Side effects: Overwrites the input data
  140.  *
  141.  * --------------------------------------------------------------
  142.  */
  143.  
  144. void
  145. mp_fwd_dct_block(data2d)
  146.     Block data2d;
  147. {
  148.     int16 *data = (int16 *) data2d;    /* this algorithm wants
  149.                      * a 1-d array */
  150.     int pass, rowctr;
  151.     register int16 *inptr, *outptr;
  152.     int16 workspace[DCTSIZE_SQ];
  153.     SHIFT_TEMPS
  154.  
  155. #ifdef ndef
  156.     {
  157.     int y;
  158.  
  159.     printf("fwd_dct (beforehand):\n");
  160.     for (y = 0; y < 8; y++)
  161.         printf("%4d %4d %4d %4d %4d %4d %4d %4d\n",
  162.            data2d[y][0], data2d[y][1],
  163.            data2d[y][2], data2d[y][3],
  164.            data2d[y][4], data2d[y][5],
  165.            data2d[y][6], data2d[y][7]);
  166.     }
  167. #endif
  168.  
  169.     /*
  170.      * Each iteration of the inner loop performs one 8-point 1-D DCT. It
  171.      * reads from a *row* of the input matrix and stores into a *column*
  172.      * of the output matrix.  In the first pass, we read from the data[]
  173.      * array and store into the local workspace[].  In the second pass,
  174.      * we read from the workspace[] array and store into data[], thus
  175.      * performing the equivalent of a columnar DCT pass with no variable
  176.      * array indexing.
  177.      */
  178.  
  179.     inptr = data;        /* initialize pointers for first pass */
  180.     outptr = workspace;
  181.     for (pass = 1; pass >= 0; pass--) {
  182.     for (rowctr = DCTSIZE - 1; rowctr >= 0; rowctr--) {
  183.         /*
  184.          * many tmps have nonoverlapping lifetime -- flashy
  185.          * register colourers should be able to do this lot
  186.          * very well
  187.          */
  188.         int32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
  189.         int32 tmp10, tmp11, tmp12, tmp13;
  190.         int32 tmp14, tmp15, tmp16, tmp17;
  191.         int32 tmp25, tmp26;
  192.         /* SHIFT_TEMPS */
  193.  
  194.         /* temp0 through tmp7:  -512 to +512 */
  195.         /* if I-block, then -256 to +256 */
  196.         tmp0 = inptr[7] + inptr[0];
  197.         tmp1 = inptr[6] + inptr[1];
  198.         tmp2 = inptr[5] + inptr[2];
  199.         tmp3 = inptr[4] + inptr[3];
  200.         tmp4 = inptr[3] - inptr[4];
  201.         tmp5 = inptr[2] - inptr[5];
  202.         tmp6 = inptr[1] - inptr[6];
  203.         tmp7 = inptr[0] - inptr[7];
  204.  
  205.         /* tmp10 through tmp13:  -1024 to +1024 */
  206.         /* if I-block, then -512 to +512 */
  207.         tmp10 = tmp3 + tmp0;
  208.         tmp11 = tmp2 + tmp1;
  209.         tmp12 = tmp1 - tmp2;
  210.         tmp13 = tmp0 - tmp3;
  211.  
  212.         outptr[0] = (int16) UNFIXH((tmp10 + tmp11) * SIN_1_4);
  213.         outptr[DCTSIZE * 4] = (int16) UNFIXH((tmp10 - tmp11) * COS_1_4);
  214.  
  215.         outptr[DCTSIZE * 2] = (int16) UNFIXH(tmp13 * COS_1_8 + tmp12 * SIN_1_8);
  216.         outptr[DCTSIZE * 6] = (int16) UNFIXH(tmp13 * SIN_1_8 - tmp12 * COS_1_8);
  217.  
  218.         tmp16 = UNFIXO((tmp6 + tmp5) * SIN_1_4);
  219.         tmp15 = UNFIXO((tmp6 - tmp5) * COS_1_4);
  220.  
  221.         OVERSHIFT(tmp4);
  222.         OVERSHIFT(tmp7);
  223.  
  224.         /*
  225.          * tmp4, tmp7, tmp15, tmp16 are overscaled by
  226.          * OVERSCALE
  227.          */
  228.  
  229.         tmp14 = tmp4 + tmp15;
  230.         tmp25 = tmp4 - tmp15;
  231.         tmp26 = tmp7 - tmp16;
  232.         tmp17 = tmp7 + tmp16;
  233.  
  234.         outptr[DCTSIZE] = (int16) UNFIXH(tmp17 * OCOS_1_16 + tmp14 * OSIN_1_16);
  235.         outptr[DCTSIZE * 7] = (int16) UNFIXH(tmp17 * OCOS_7_16 - tmp14 * OSIN_7_16);
  236.         outptr[DCTSIZE * 5] = (int16) UNFIXH(tmp26 * OCOS_5_16 + tmp25 * OSIN_5_16);
  237.         outptr[DCTSIZE * 3] = (int16) UNFIXH(tmp26 * OCOS_3_16 - tmp25 * OSIN_3_16);
  238.  
  239.         inptr += DCTSIZE;    /* advance inptr to next row */
  240.         outptr++;        /* advance outptr to next column */
  241.     }
  242.     /* end of pass; in case it was pass 1, set up for pass 2 */
  243.     inptr = workspace;
  244.     outptr = data;
  245.     }
  246. #ifdef ndef
  247.     {
  248.     int y;
  249.  
  250.     printf("fwd_dct (afterward):\n");
  251.     for (y = 0; y < 8; y++)
  252.         printf("%4d %4d %4d %4d %4d %4d %4d %4d\n",
  253.            data2d[y][0], data2d[y][1],
  254.            data2d[y][2], data2d[y][3],
  255.            data2d[y][4], data2d[y][5],
  256.            data2d[y][6], data2d[y][7]);
  257.     }
  258. #endif
  259. }
  260.