Thomas G. Lane | 2cbeb8a | 1991-10-07 00:00:00 +0000 | [diff] [blame^] | 1 | /* |
| 2 | * jfwddct.c |
| 3 | * |
| 4 | * Copyright (C) 1991, Thomas G. Lane. |
| 5 | * This file is part of the Independent JPEG Group's software. |
| 6 | * For conditions of distribution and use, see the accompanying README file. |
| 7 | * |
| 8 | * This file contains the basic DCT (Discrete Cosine Transform) |
| 9 | * transformation subroutine. |
| 10 | * |
| 11 | * This implementation is based on Appendix A.2 of the book |
| 12 | * "Discrete Cosine Transform---Algorithms, Advantages, Applications" |
| 13 | * by K.R. Rao and P. Yip (Academic Press, Inc, London, 1990). |
| 14 | * It uses scaled fixed-point arithmetic instead of floating point. |
| 15 | */ |
| 16 | |
| 17 | #include "jinclude.h" |
| 18 | |
| 19 | |
| 20 | /* The poop on this scaling stuff is as follows: |
| 21 | * |
| 22 | * Most of the numbers (after multiplication by the constants) are |
| 23 | * (logically) shifted left by LG2_DCT_SCALE. This is undone by UNFIXH |
| 24 | * before assignment to the output array. Note that we want an additional |
| 25 | * division by 2 on the output (required by the equations). |
| 26 | * |
| 27 | * If right shifts are unsigned, then there is a potential problem. |
| 28 | * However, shifting right by 16 and then assigning to a short |
| 29 | * (assuming short = 16 bits) will keep the sign right!! |
| 30 | * |
| 31 | * For other shifts, |
| 32 | * |
| 33 | * ((x + (1 << 30)) >> shft) - (1 << (30 - shft)) |
| 34 | * |
| 35 | * gives a nice right shift with sign (assuming no overflow). However, all the |
| 36 | * scaling is such that this isn't a problem. (Is this true?) |
| 37 | */ |
| 38 | |
| 39 | |
| 40 | #define ONE 1L /* remove L if long > 32 bits */ |
| 41 | |
| 42 | #ifdef RIGHT_SHIFT_IS_UNSIGNED |
| 43 | #define LG2_DCT_SCALE 15 |
| 44 | #define RIGHT_SHIFT(_x,_shft) ((((_x) + (ONE << 30)) >> (_shft)) - (ONE << (30 - (_shft)))) |
| 45 | #else |
| 46 | #define LG2_DCT_SCALE 16 |
| 47 | #define RIGHT_SHIFT(_x,_shft) ((_x) >> (_shft)) |
| 48 | #endif |
| 49 | |
| 50 | #define DCT_SCALE (ONE << LG2_DCT_SCALE) |
| 51 | |
| 52 | #define LG2_OVERSCALE 2 |
| 53 | #define OVERSCALE (ONE << LG2_OVERSCALE) |
| 54 | |
| 55 | #define FIX(x) ((INT32) ((x) * DCT_SCALE + 0.5)) |
| 56 | #define FIXO(x) ((INT32) ((x) * DCT_SCALE / OVERSCALE + 0.5)) |
| 57 | #define UNFIX(x) RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1)), LG2_DCT_SCALE) |
| 58 | #define UNFIXH(x) RIGHT_SHIFT((x) + (ONE << LG2_DCT_SCALE), LG2_DCT_SCALE+1) |
| 59 | #define UNFIXO(x) RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1-LG2_OVERSCALE)), LG2_DCT_SCALE-LG2_OVERSCALE) |
| 60 | #define OVERSH(x) ((x) << LG2_OVERSCALE) |
| 61 | |
| 62 | #define SIN_1_4 FIX(0.7071067811856476) |
| 63 | #define COS_1_4 SIN_1_4 |
| 64 | |
| 65 | #define SIN_1_8 FIX(0.3826834323650898) |
| 66 | #define COS_1_8 FIX(0.9238795325112870) |
| 67 | #define SIN_3_8 COS_1_8 |
| 68 | #define COS_3_8 SIN_1_8 |
| 69 | |
| 70 | #define SIN_1_16 FIX(0.1950903220161282) |
| 71 | #define COS_1_16 FIX(0.9807852804032300) |
| 72 | #define SIN_7_16 COS_1_16 |
| 73 | #define COS_7_16 SIN_1_16 |
| 74 | |
| 75 | #define SIN_3_16 FIX(0.5555702330196022) |
| 76 | #define COS_3_16 FIX(0.8314696123025450) |
| 77 | #define SIN_5_16 COS_3_16 |
| 78 | #define COS_5_16 SIN_3_16 |
| 79 | |
| 80 | #define OSIN_1_4 FIXO(0.707106781185647) |
| 81 | #define OCOS_1_4 OSIN_1_4 |
| 82 | |
| 83 | #define OSIN_1_8 FIXO(0.3826834323650898) |
| 84 | #define OCOS_1_8 FIXO(0.9238795325112870) |
| 85 | #define OSIN_3_8 OCOS_1_8 |
| 86 | #define OCOS_3_8 OSIN_1_8 |
| 87 | |
| 88 | #define OSIN_1_16 FIXO(0.1950903220161282) |
| 89 | #define OCOS_1_16 FIXO(0.9807852804032300) |
| 90 | #define OSIN_7_16 OCOS_1_16 |
| 91 | #define OCOS_7_16 OSIN_1_16 |
| 92 | |
| 93 | #define OSIN_3_16 FIXO(0.5555702330196022) |
| 94 | #define OCOS_3_16 FIXO(0.8314696123025450) |
| 95 | #define OSIN_5_16 OCOS_3_16 |
| 96 | #define OCOS_5_16 OSIN_3_16 |
| 97 | |
| 98 | |
| 99 | INLINE |
| 100 | LOCAL void |
| 101 | fast_dct_8 (DCTELEM *in, int stride) |
| 102 | { |
| 103 | /* tmp1x are new values of tmpx -- flashy register colourers |
| 104 | * should be able to do this lot very well |
| 105 | */ |
| 106 | INT32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7; |
| 107 | INT32 tmp10, tmp11, tmp12, tmp13, tmp14, tmp15, tmp16, tmp17; |
| 108 | INT32 tmp25, tmp26; |
| 109 | INT32 in0, in1, in2, in3, in4, in5, in6, in7; |
| 110 | |
| 111 | in0 = in[ 0]; |
| 112 | in1 = in[stride ]; |
| 113 | in2 = in[stride*2]; |
| 114 | in3 = in[stride*3]; |
| 115 | in4 = in[stride*4]; |
| 116 | in5 = in[stride*5]; |
| 117 | in6 = in[stride*6]; |
| 118 | in7 = in[stride*7]; |
| 119 | |
| 120 | tmp0 = in7 + in0; |
| 121 | tmp1 = in6 + in1; |
| 122 | tmp2 = in5 + in2; |
| 123 | tmp3 = in4 + in3; |
| 124 | tmp4 = in3 - in4; |
| 125 | tmp5 = in2 - in5; |
| 126 | tmp6 = in1 - in6; |
| 127 | tmp7 = in0 - in7; |
| 128 | |
| 129 | tmp10 = tmp3 + tmp0 ; |
| 130 | tmp11 = tmp2 + tmp1 ; |
| 131 | tmp12 = tmp1 - tmp2 ; |
| 132 | tmp13 = tmp0 - tmp3 ; |
| 133 | |
| 134 | /* Now using tmp10, tmp11, tmp12, tmp13 */ |
| 135 | |
| 136 | in[ 0] = UNFIXH((tmp10 + tmp11) * SIN_1_4); |
| 137 | in[stride*4] = UNFIXH((tmp10 - tmp11) * COS_1_4); |
| 138 | |
| 139 | in[stride*2] = UNFIXH(tmp13*COS_1_8 + tmp12*SIN_1_8); |
| 140 | in[stride*6] = UNFIXH(tmp13*SIN_1_8 - tmp12*COS_1_8); |
| 141 | |
| 142 | tmp16 = UNFIXO((tmp6 + tmp5) * SIN_1_4); |
| 143 | tmp15 = UNFIXO((tmp6 - tmp5) * COS_1_4); |
| 144 | |
| 145 | /* Now using tmp10, tmp11, tmp13, tmp14, tmp15, tmp16 */ |
| 146 | |
| 147 | tmp14 = OVERSH(tmp4) + tmp15; |
| 148 | tmp25 = OVERSH(tmp4) - tmp15; |
| 149 | tmp26 = OVERSH(tmp7) - tmp16; |
| 150 | tmp17 = OVERSH(tmp7) + tmp16; |
| 151 | |
| 152 | /* These are now overscaled by OVERSCALE */ |
| 153 | |
| 154 | /* tmp10, tmp11, tmp12, tmp13, tmp14, tmp25, tmp26, tmp17 */ |
| 155 | |
| 156 | in[stride ] = UNFIXH(tmp17*OCOS_1_16 + tmp14*OSIN_1_16); |
| 157 | in[stride*7] = UNFIXH(tmp17*OCOS_7_16 - tmp14*OSIN_7_16); |
| 158 | in[stride*5] = UNFIXH(tmp26*OCOS_5_16 + tmp25*OSIN_5_16); |
| 159 | in[stride*3] = UNFIXH(tmp26*OCOS_3_16 - tmp25*OSIN_3_16); |
| 160 | } |
| 161 | |
| 162 | |
| 163 | /* |
| 164 | * Perform the forward DCT on one block of samples. |
| 165 | * |
| 166 | * Note that this code is specialized to the case DCTSIZE = 8. |
| 167 | */ |
| 168 | |
| 169 | GLOBAL void |
| 170 | j_fwd_dct (DCTBLOCK data) |
| 171 | { |
| 172 | int i; |
| 173 | |
| 174 | for (i = 0; i < DCTSIZE; i++) |
| 175 | fast_dct_8(data+i*DCTSIZE, 1); |
| 176 | |
| 177 | for (i = 0; i < DCTSIZE; i++) |
| 178 | fast_dct_8(data+i, DCTSIZE); |
| 179 | } |