4 * Copyright (C) 1994-1996, Thomas G. Lane.
5 * Modified 2003-2009 by Guido Vollbeding.
6 * This file is part of the Independent JPEG Group's software.
7 * For conditions of distribution and use, see the accompanying README file.
9 * This file contains a floating-point implementation of the
10 * forward DCT (Discrete Cosine Transform).
12 * This implementation should be more accurate than either of the integer
13 * DCT implementations. However, it may not give the same results on all
14 * machines because of differences in roundoff behavior. Speed will depend
15 * on the hardware's floating point capacity.
17 * A 2-D DCT can be done by 1-D DCT on each row followed by 1-D DCT
18 * on each column. Direct algorithms are also available, but they are
19 * much more complex and seem not to be any faster when reduced to code.
21 * This implementation is based on Arai, Agui, and Nakajima's algorithm for
22 * scaled DCT. Their original paper (Trans. IEICE E-71(11):1095) is in
23 * Japanese, but the algorithm is described in the Pennebaker & Mitchell
24 * JPEG textbook (see REFERENCES section in file README). The following code
25 * is based directly on figure 4-8 in P&M.
26 * While an 8-point DCT cannot be done in less than 11 multiplies, it is
27 * possible to arrange the computation so that many of the multiplies are
28 * simple scalings of the final outputs. These multiplies can then be
29 * folded into the multiplications or divisions by the JPEG quantization
30 * table entries. The AA&N method leaves only 5 multiplies and 29 adds
31 * to be done in the DCT itself.
32 * The primary disadvantage of this method is that with a fixed-point
33 * implementation, accuracy is lost due to imprecise representation of the
34 * scaled quantization values. However, that problem does not arise if
35 * we use floating point arithmetic.
38 #define JPEG_INTERNALS
41 #include "jdct.h" /* Private declarations for DCT subsystem */
43 #ifdef DCT_FLOAT_SUPPORTED
47 * This module is specialized to the case DCTSIZE = 8.
51 Sorry
, this code only copes with
8x8 DCTs
. /* deliberate syntax err */
56 * Perform the forward DCT on one block of samples.
60 jpeg_fdct_float (FAST_FLOAT
* data
, JSAMPARRAY sample_data
, JDIMENSION start_col
)
62 FAST_FLOAT tmp0
, tmp1
, tmp2
, tmp3
, tmp4
, tmp5
, tmp6
, tmp7
;
63 FAST_FLOAT tmp10
, tmp11
, tmp12
, tmp13
;
64 FAST_FLOAT z1
, z2
, z3
, z4
, z5
, z11
, z13
;
69 /* Pass 1: process rows. */
72 for (ctr
= 0; ctr
< DCTSIZE
; ctr
++) {
73 elemptr
= sample_data
[ctr
] + start_col
;
75 /* Load data into workspace */
76 tmp0
= (FAST_FLOAT
) (GETJSAMPLE(elemptr
[0]) + GETJSAMPLE(elemptr
[7]));
77 tmp7
= (FAST_FLOAT
) (GETJSAMPLE(elemptr
[0]) - GETJSAMPLE(elemptr
[7]));
78 tmp1
= (FAST_FLOAT
) (GETJSAMPLE(elemptr
[1]) + GETJSAMPLE(elemptr
[6]));
79 tmp6
= (FAST_FLOAT
) (GETJSAMPLE(elemptr
[1]) - GETJSAMPLE(elemptr
[6]));
80 tmp2
= (FAST_FLOAT
) (GETJSAMPLE(elemptr
[2]) + GETJSAMPLE(elemptr
[5]));
81 tmp5
= (FAST_FLOAT
) (GETJSAMPLE(elemptr
[2]) - GETJSAMPLE(elemptr
[5]));
82 tmp3
= (FAST_FLOAT
) (GETJSAMPLE(elemptr
[3]) + GETJSAMPLE(elemptr
[4]));
83 tmp4
= (FAST_FLOAT
) (GETJSAMPLE(elemptr
[3]) - GETJSAMPLE(elemptr
[4]));
87 tmp10
= tmp0
+ tmp3
; /* phase 2 */
92 /* Apply unsigned->signed conversion */
93 dataptr
[0] = tmp10
+ tmp11
- 8 * CENTERJSAMPLE
; /* phase 3 */
94 dataptr
[4] = tmp10
- tmp11
;
96 z1
= (tmp12
+ tmp13
) * ((FAST_FLOAT
) 0.707106781); /* c4 */
97 dataptr
[2] = tmp13
+ z1
; /* phase 5 */
98 dataptr
[6] = tmp13
- z1
;
102 tmp10
= tmp4
+ tmp5
; /* phase 2 */
106 /* The rotator is modified from fig 4-8 to avoid extra negations. */
107 z5
= (tmp10
- tmp12
) * ((FAST_FLOAT
) 0.382683433); /* c6 */
108 z2
= ((FAST_FLOAT
) 0.541196100) * tmp10
+ z5
; /* c2-c6 */
109 z4
= ((FAST_FLOAT
) 1.306562965) * tmp12
+ z5
; /* c2+c6 */
110 z3
= tmp11
* ((FAST_FLOAT
) 0.707106781); /* c4 */
112 z11
= tmp7
+ z3
; /* phase 5 */
115 dataptr
[5] = z13
+ z2
; /* phase 6 */
116 dataptr
[3] = z13
- z2
;
117 dataptr
[1] = z11
+ z4
;
118 dataptr
[7] = z11
- z4
;
120 dataptr
+= DCTSIZE
; /* advance pointer to next row */
123 /* Pass 2: process columns. */
126 for (ctr
= DCTSIZE
-1; ctr
>= 0; ctr
--) {
127 tmp0
= dataptr
[DCTSIZE
*0] + dataptr
[DCTSIZE
*7];
128 tmp7
= dataptr
[DCTSIZE
*0] - dataptr
[DCTSIZE
*7];
129 tmp1
= dataptr
[DCTSIZE
*1] + dataptr
[DCTSIZE
*6];
130 tmp6
= dataptr
[DCTSIZE
*1] - dataptr
[DCTSIZE
*6];
131 tmp2
= dataptr
[DCTSIZE
*2] + dataptr
[DCTSIZE
*5];
132 tmp5
= dataptr
[DCTSIZE
*2] - dataptr
[DCTSIZE
*5];
133 tmp3
= dataptr
[DCTSIZE
*3] + dataptr
[DCTSIZE
*4];
134 tmp4
= dataptr
[DCTSIZE
*3] - dataptr
[DCTSIZE
*4];
138 tmp10
= tmp0
+ tmp3
; /* phase 2 */
143 dataptr
[DCTSIZE
*0] = tmp10
+ tmp11
; /* phase 3 */
144 dataptr
[DCTSIZE
*4] = tmp10
- tmp11
;
146 z1
= (tmp12
+ tmp13
) * ((FAST_FLOAT
) 0.707106781); /* c4 */
147 dataptr
[DCTSIZE
*2] = tmp13
+ z1
; /* phase 5 */
148 dataptr
[DCTSIZE
*6] = tmp13
- z1
;
152 tmp10
= tmp4
+ tmp5
; /* phase 2 */
156 /* The rotator is modified from fig 4-8 to avoid extra negations. */
157 z5
= (tmp10
- tmp12
) * ((FAST_FLOAT
) 0.382683433); /* c6 */
158 z2
= ((FAST_FLOAT
) 0.541196100) * tmp10
+ z5
; /* c2-c6 */
159 z4
= ((FAST_FLOAT
) 1.306562965) * tmp12
+ z5
; /* c2+c6 */
160 z3
= tmp11
* ((FAST_FLOAT
) 0.707106781); /* c4 */
162 z11
= tmp7
+ z3
; /* phase 5 */
165 dataptr
[DCTSIZE
*5] = z13
+ z2
; /* phase 6 */
166 dataptr
[DCTSIZE
*3] = z13
- z2
;
167 dataptr
[DCTSIZE
*1] = z11
+ z4
;
168 dataptr
[DCTSIZE
*7] = z11
- z4
;
170 dataptr
++; /* advance pointer to next column */
174 #endif /* DCT_FLOAT_SUPPORTED */