1 /*
2 Short Discrete Cosine Transform
3 data length :8x8, 16x16
4 method :row-column, radix 4 FFT
5 functions
6 ddct8x8s : 8x8 DCT
7 ddct16x16s: 16x16 DCT
8 function prototypes
9 void ddct8x8s(int isgn, double **a);
10 void ddct16x16s(int isgn, double **a);
11 */
12
13
14 /*
15 -------- 8x8 DCT (Discrete Cosine Transform) / Inverse of DCT --------
16 [definition]
17 <case1> Normalized 8x8 IDCT
18 C[k1][k2] = (1/4) * sum_j1=0^7 sum_j2=0^7
19 a[j1][j2] * s[j1] * s[j2] *
20 cos(pi*j1*(k1+1/2)/8) *
21 cos(pi*j2*(k2+1/2)/8), 0<=k1<8, 0<=k2<8
22 (s[0] = 1/sqrt(2), s[j] = 1, j > 0)
23 <case2> Normalized 8x8 DCT
24 C[k1][k2] = (1/4) * s[k1] * s[k2] * sum_j1=0^7 sum_j2=0^7
25 a[j1][j2] *
26 cos(pi*(j1+1/2)*k1/8) *
27 cos(pi*(j2+1/2)*k2/8), 0<=k1<8, 0<=k2<8
28 (s[0] = 1/sqrt(2), s[j] = 1, j > 0)
29 [usage]
30 <case1>
31 ddct8x8s(1, a);
32 <case2>
33 ddct8x8s(-1, a);
34 [parameters]
35 a[0...7][0...7] :input/output data (double **)
36 output data
37 a[k1][k2] = C[k1][k2], 0<=k1<8, 0<=k2<8
38 */
39
40
41 /* Cn_kR = sqrt(2.0/n) * cos(pi/2*k/n) */
42 /* Cn_kI = sqrt(2.0/n) * sin(pi/2*k/n) */
43 /* Wn_kR = cos(pi/2*k/n) */
44 /* Wn_kI = sin(pi/2*k/n) */
45 #define C8_1R 0.49039264020161522456
46 #define C8_1I 0.09754516100806413392
47 #define C8_2R 0.46193976625564337806
48 #define C8_2I 0.19134171618254488586
49 #define C8_3R 0.41573480615127261854
50 #define C8_3I 0.27778511650980111237
51 #define C8_4R 0.35355339059327376220
52 #define W8_4R 0.70710678118654752440
53
54
ddct8x8s(int isgn,double ** a)55 void ddct8x8s(int isgn, double **a)
56 {
57 int j;
58 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
59 double xr, xi;
60
61 if (isgn < 0) {
62 for (j = 0; j <= 7; j++) {
63 x0r = a[0][j] + a[7][j];
64 x1r = a[0][j] - a[7][j];
65 x0i = a[2][j] + a[5][j];
66 x1i = a[2][j] - a[5][j];
67 x2r = a[4][j] + a[3][j];
68 x3r = a[4][j] - a[3][j];
69 x2i = a[6][j] + a[1][j];
70 x3i = a[6][j] - a[1][j];
71 xr = x0r + x2r;
72 xi = x0i + x2i;
73 a[0][j] = C8_4R * (xr + xi);
74 a[4][j] = C8_4R * (xr - xi);
75 xr = x0r - x2r;
76 xi = x0i - x2i;
77 a[2][j] = C8_2R * xr - C8_2I * xi;
78 a[6][j] = C8_2R * xi + C8_2I * xr;
79 xr = W8_4R * (x1i - x3i);
80 x1i = W8_4R * (x1i + x3i);
81 x3i = x1i - x3r;
82 x1i += x3r;
83 x3r = x1r - xr;
84 x1r += xr;
85 a[1][j] = C8_1R * x1r - C8_1I * x1i;
86 a[7][j] = C8_1R * x1i + C8_1I * x1r;
87 a[3][j] = C8_3R * x3r - C8_3I * x3i;
88 a[5][j] = C8_3R * x3i + C8_3I * x3r;
89 }
90 for (j = 0; j <= 7; j++) {
91 x0r = a[j][0] + a[j][7];
92 x1r = a[j][0] - a[j][7];
93 x0i = a[j][2] + a[j][5];
94 x1i = a[j][2] - a[j][5];
95 x2r = a[j][4] + a[j][3];
96 x3r = a[j][4] - a[j][3];
97 x2i = a[j][6] + a[j][1];
98 x3i = a[j][6] - a[j][1];
99 xr = x0r + x2r;
100 xi = x0i + x2i;
101 a[j][0] = C8_4R * (xr + xi);
102 a[j][4] = C8_4R * (xr - xi);
103 xr = x0r - x2r;
104 xi = x0i - x2i;
105 a[j][2] = C8_2R * xr - C8_2I * xi;
106 a[j][6] = C8_2R * xi + C8_2I * xr;
107 xr = W8_4R * (x1i - x3i);
108 x1i = W8_4R * (x1i + x3i);
109 x3i = x1i - x3r;
110 x1i += x3r;
111 x3r = x1r - xr;
112 x1r += xr;
113 a[j][1] = C8_1R * x1r - C8_1I * x1i;
114 a[j][7] = C8_1R * x1i + C8_1I * x1r;
115 a[j][3] = C8_3R * x3r - C8_3I * x3i;
116 a[j][5] = C8_3R * x3i + C8_3I * x3r;
117 }
118 } else {
119 for (j = 0; j <= 7; j++) {
120 x1r = C8_1R * a[1][j] + C8_1I * a[7][j];
121 x1i = C8_1R * a[7][j] - C8_1I * a[1][j];
122 x3r = C8_3R * a[3][j] + C8_3I * a[5][j];
123 x3i = C8_3R * a[5][j] - C8_3I * a[3][j];
124 xr = x1r - x3r;
125 xi = x1i + x3i;
126 x1r += x3r;
127 x3i -= x1i;
128 x1i = W8_4R * (xr + xi);
129 x3r = W8_4R * (xr - xi);
130 xr = C8_2R * a[2][j] + C8_2I * a[6][j];
131 xi = C8_2R * a[6][j] - C8_2I * a[2][j];
132 x0r = C8_4R * (a[0][j] + a[4][j]);
133 x0i = C8_4R * (a[0][j] - a[4][j]);
134 x2r = x0r - xr;
135 x2i = x0i - xi;
136 x0r += xr;
137 x0i += xi;
138 a[0][j] = x0r + x1r;
139 a[7][j] = x0r - x1r;
140 a[2][j] = x0i + x1i;
141 a[5][j] = x0i - x1i;
142 a[4][j] = x2r - x3i;
143 a[3][j] = x2r + x3i;
144 a[6][j] = x2i - x3r;
145 a[1][j] = x2i + x3r;
146 }
147 for (j = 0; j <= 7; j++) {
148 x1r = C8_1R * a[j][1] + C8_1I * a[j][7];
149 x1i = C8_1R * a[j][7] - C8_1I * a[j][1];
150 x3r = C8_3R * a[j][3] + C8_3I * a[j][5];
151 x3i = C8_3R * a[j][5] - C8_3I * a[j][3];
152 xr = x1r - x3r;
153 xi = x1i + x3i;
154 x1r += x3r;
155 x3i -= x1i;
156 x1i = W8_4R * (xr + xi);
157 x3r = W8_4R * (xr - xi);
158 xr = C8_2R * a[j][2] + C8_2I * a[j][6];
159 xi = C8_2R * a[j][6] - C8_2I * a[j][2];
160 x0r = C8_4R * (a[j][0] + a[j][4]);
161 x0i = C8_4R * (a[j][0] - a[j][4]);
162 x2r = x0r - xr;
163 x2i = x0i - xi;
164 x0r += xr;
165 x0i += xi;
166 a[j][0] = x0r + x1r;
167 a[j][7] = x0r - x1r;
168 a[j][2] = x0i + x1i;
169 a[j][5] = x0i - x1i;
170 a[j][4] = x2r - x3i;
171 a[j][3] = x2r + x3i;
172 a[j][6] = x2i - x3r;
173 a[j][1] = x2i + x3r;
174 }
175 }
176 }
177
178
179
180 /*
181 -------- 16x16 DCT (Discrete Cosine Transform) / Inverse of DCT --------
182 [definition]
183 <case1> Normalized 16x16 IDCT
184 C[k1][k2] = (1/8) * sum_j1=0^15 sum_j2=0^15
185 a[j1][j2] * s[j1] * s[j2] *
186 cos(pi*j1*(k1+1/2)/16) *
187 cos(pi*j2*(k2+1/2)/16), 0<=k1<16, 0<=k2<16
188 (s[0] = 1/sqrt(2), s[j] = 1, j > 0)
189 <case2> Normalized 16x16 DCT
190 C[k1][k2] = (1/8) * s[k1] * s[k2] * sum_j1=0^15 sum_j2=0^15
191 a[j1][j2] *
192 cos(pi*(j1+1/2)*k1/16) *
193 cos(pi*(j2+1/2)*k2/16), 0<=k1<16, 0<=k2<16
194 (s[0] = 1/sqrt(2), s[j] = 1, j > 0)
195 [usage]
196 <case1>
197 ddct16x16s(1, a);
198 <case2>
199 ddct16x16s(-1, a);
200 [parameters]
201 a[0...15][0...15] :input/output data (double **)
202 output data
203 a[k1][k2] = C[k1][k2], 0<=k1<16, 0<=k2<16
204 */
205
206
207 /* Cn_kR = sqrt(2.0/n) * cos(pi/2*k/n) */
208 /* Cn_kI = sqrt(2.0/n) * sin(pi/2*k/n) */
209 /* Wn_kR = cos(pi/2*k/n) */
210 /* Wn_kI = sin(pi/2*k/n) */
211 #define C16_1R 0.35185093438159561476
212 #define C16_1I 0.03465429229977286565
213 #define C16_2R 0.34675996133053686546
214 #define C16_2I 0.06897484482073575308
215 #define C16_3R 0.33832950029358816957
216 #define C16_3I 0.10263113188058934529
217 #define C16_4R 0.32664074121909413196
218 #define C16_4I 0.13529902503654924610
219 #define C16_5R 0.31180625324666780814
220 #define C16_5I 0.16666391461943662432
221 #define C16_6R 0.29396890060483967924
222 #define C16_6I 0.19642373959677554532
223 #define C16_7R 0.27330046675043937206
224 #define C16_7I 0.22429189658565907106
225 #define C16_8R 0.25
226 #define W16_4R 0.92387953251128675613
227 #define W16_4I 0.38268343236508977173
228 #define W16_8R 0.70710678118654752440
229
230
ddct16x16s(int isgn,double ** a)231 void ddct16x16s(int isgn, double **a)
232 {
233 int j;
234 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
235 double x4r, x4i, x5r, x5i, x6r, x6i, x7r, x7i;
236 double xr, xi;
237
238 if (isgn < 0) {
239 for (j = 0; j <= 15; j++) {
240 x4r = a[0][j] - a[15][j];
241 xr = a[0][j] + a[15][j];
242 x4i = a[8][j] - a[7][j];
243 xi = a[8][j] + a[7][j];
244 x0r = xr + xi;
245 x0i = xr - xi;
246 x5r = a[2][j] - a[13][j];
247 xr = a[2][j] + a[13][j];
248 x5i = a[10][j] - a[5][j];
249 xi = a[10][j] + a[5][j];
250 x1r = xr + xi;
251 x1i = xr - xi;
252 x6r = a[4][j] - a[11][j];
253 xr = a[4][j] + a[11][j];
254 x6i = a[12][j] - a[3][j];
255 xi = a[12][j] + a[3][j];
256 x2r = xr + xi;
257 x2i = xr - xi;
258 x7r = a[6][j] - a[9][j];
259 xr = a[6][j] + a[9][j];
260 x7i = a[14][j] - a[1][j];
261 xi = a[14][j] + a[1][j];
262 x3r = xr + xi;
263 x3i = xr - xi;
264 xr = x0r + x2r;
265 xi = x1r + x3r;
266 a[0][j] = C16_8R * (xr + xi);
267 a[8][j] = C16_8R * (xr - xi);
268 xr = x0r - x2r;
269 xi = x1r - x3r;
270 a[4][j] = C16_4R * xr - C16_4I * xi;
271 a[12][j] = C16_4R * xi + C16_4I * xr;
272 x0r = W16_8R * (x1i - x3i);
273 x2r = W16_8R * (x1i + x3i);
274 xr = x0i + x0r;
275 xi = x2r + x2i;
276 a[2][j] = C16_2R * xr - C16_2I * xi;
277 a[14][j] = C16_2R * xi + C16_2I * xr;
278 xr = x0i - x0r;
279 xi = x2r - x2i;
280 a[6][j] = C16_6R * xr - C16_6I * xi;
281 a[10][j] = C16_6R * xi + C16_6I * xr;
282 xr = W16_8R * (x6r - x6i);
283 xi = W16_8R * (x6i + x6r);
284 x6r = x4r - xr;
285 x6i = x4i - xi;
286 x4r += xr;
287 x4i += xi;
288 xr = W16_4I * x7r - W16_4R * x7i;
289 xi = W16_4I * x7i + W16_4R * x7r;
290 x7r = W16_4R * x5r - W16_4I * x5i;
291 x7i = W16_4R * x5i + W16_4I * x5r;
292 x5r = x7r + xr;
293 x5i = x7i + xi;
294 x7r -= xr;
295 x7i -= xi;
296 xr = x4r + x5r;
297 xi = x5i + x4i;
298 a[1][j] = C16_1R * xr - C16_1I * xi;
299 a[15][j] = C16_1R * xi + C16_1I * xr;
300 xr = x4r - x5r;
301 xi = x5i - x4i;
302 a[7][j] = C16_7R * xr - C16_7I * xi;
303 a[9][j] = C16_7R * xi + C16_7I * xr;
304 xr = x6r - x7i;
305 xi = x7r + x6i;
306 a[5][j] = C16_5R * xr - C16_5I * xi;
307 a[11][j] = C16_5R * xi + C16_5I * xr;
308 xr = x6r + x7i;
309 xi = x7r - x6i;
310 a[3][j] = C16_3R * xr - C16_3I * xi;
311 a[13][j] = C16_3R * xi + C16_3I * xr;
312 }
313 for (j = 0; j <= 15; j++) {
314 x4r = a[j][0] - a[j][15];
315 xr = a[j][0] + a[j][15];
316 x4i = a[j][8] - a[j][7];
317 xi = a[j][8] + a[j][7];
318 x0r = xr + xi;
319 x0i = xr - xi;
320 x5r = a[j][2] - a[j][13];
321 xr = a[j][2] + a[j][13];
322 x5i = a[j][10] - a[j][5];
323 xi = a[j][10] + a[j][5];
324 x1r = xr + xi;
325 x1i = xr - xi;
326 x6r = a[j][4] - a[j][11];
327 xr = a[j][4] + a[j][11];
328 x6i = a[j][12] - a[j][3];
329 xi = a[j][12] + a[j][3];
330 x2r = xr + xi;
331 x2i = xr - xi;
332 x7r = a[j][6] - a[j][9];
333 xr = a[j][6] + a[j][9];
334 x7i = a[j][14] - a[j][1];
335 xi = a[j][14] + a[j][1];
336 x3r = xr + xi;
337 x3i = xr - xi;
338 xr = x0r + x2r;
339 xi = x1r + x3r;
340 a[j][0] = C16_8R * (xr + xi);
341 a[j][8] = C16_8R * (xr - xi);
342 xr = x0r - x2r;
343 xi = x1r - x3r;
344 a[j][4] = C16_4R * xr - C16_4I * xi;
345 a[j][12] = C16_4R * xi + C16_4I * xr;
346 x0r = W16_8R * (x1i - x3i);
347 x2r = W16_8R * (x1i + x3i);
348 xr = x0i + x0r;
349 xi = x2r + x2i;
350 a[j][2] = C16_2R * xr - C16_2I * xi;
351 a[j][14] = C16_2R * xi + C16_2I * xr;
352 xr = x0i - x0r;
353 xi = x2r - x2i;
354 a[j][6] = C16_6R * xr - C16_6I * xi;
355 a[j][10] = C16_6R * xi + C16_6I * xr;
356 xr = W16_8R * (x6r - x6i);
357 xi = W16_8R * (x6i + x6r);
358 x6r = x4r - xr;
359 x6i = x4i - xi;
360 x4r += xr;
361 x4i += xi;
362 xr = W16_4I * x7r - W16_4R * x7i;
363 xi = W16_4I * x7i + W16_4R * x7r;
364 x7r = W16_4R * x5r - W16_4I * x5i;
365 x7i = W16_4R * x5i + W16_4I * x5r;
366 x5r = x7r + xr;
367 x5i = x7i + xi;
368 x7r -= xr;
369 x7i -= xi;
370 xr = x4r + x5r;
371 xi = x5i + x4i;
372 a[j][1] = C16_1R * xr - C16_1I * xi;
373 a[j][15] = C16_1R * xi + C16_1I * xr;
374 xr = x4r - x5r;
375 xi = x5i - x4i;
376 a[j][7] = C16_7R * xr - C16_7I * xi;
377 a[j][9] = C16_7R * xi + C16_7I * xr;
378 xr = x6r - x7i;
379 xi = x7r + x6i;
380 a[j][5] = C16_5R * xr - C16_5I * xi;
381 a[j][11] = C16_5R * xi + C16_5I * xr;
382 xr = x6r + x7i;
383 xi = x7r - x6i;
384 a[j][3] = C16_3R * xr - C16_3I * xi;
385 a[j][13] = C16_3R * xi + C16_3I * xr;
386 }
387 } else {
388 for (j = 0; j <= 15; j++) {
389 x5r = C16_1R * a[1][j] + C16_1I * a[15][j];
390 x5i = C16_1R * a[15][j] - C16_1I * a[1][j];
391 xr = C16_7R * a[7][j] + C16_7I * a[9][j];
392 xi = C16_7R * a[9][j] - C16_7I * a[7][j];
393 x4r = x5r + xr;
394 x4i = x5i - xi;
395 x5r -= xr;
396 x5i += xi;
397 x7r = C16_5R * a[5][j] + C16_5I * a[11][j];
398 x7i = C16_5R * a[11][j] - C16_5I * a[5][j];
399 xr = C16_3R * a[3][j] + C16_3I * a[13][j];
400 xi = C16_3R * a[13][j] - C16_3I * a[3][j];
401 x6r = x7r + xr;
402 x6i = x7i - xi;
403 x7r -= xr;
404 x7i += xi;
405 xr = x4r - x6r;
406 xi = x4i - x6i;
407 x4r += x6r;
408 x4i += x6i;
409 x6r = W16_8R * (xi + xr);
410 x6i = W16_8R * (xi - xr);
411 xr = x5r + x7i;
412 xi = x5i - x7r;
413 x5r -= x7i;
414 x5i += x7r;
415 x7r = W16_4I * x5r + W16_4R * x5i;
416 x7i = W16_4I * x5i - W16_4R * x5r;
417 x5r = W16_4R * xr + W16_4I * xi;
418 x5i = W16_4R * xi - W16_4I * xr;
419 xr = C16_4R * a[4][j] + C16_4I * a[12][j];
420 xi = C16_4R * a[12][j] - C16_4I * a[4][j];
421 x2r = C16_8R * (a[0][j] + a[8][j]);
422 x3r = C16_8R * (a[0][j] - a[8][j]);
423 x0r = x2r + xr;
424 x1r = x3r + xi;
425 x2r -= xr;
426 x3r -= xi;
427 x0i = C16_2R * a[2][j] + C16_2I * a[14][j];
428 x2i = C16_2R * a[14][j] - C16_2I * a[2][j];
429 x1i = C16_6R * a[6][j] + C16_6I * a[10][j];
430 x3i = C16_6R * a[10][j] - C16_6I * a[6][j];
431 xr = x0i - x1i;
432 xi = x2i + x3i;
433 x0i += x1i;
434 x2i -= x3i;
435 x1i = W16_8R * (xi + xr);
436 x3i = W16_8R * (xi - xr);
437 xr = x0r + x0i;
438 xi = x0r - x0i;
439 a[0][j] = xr + x4r;
440 a[15][j] = xr - x4r;
441 a[8][j] = xi + x4i;
442 a[7][j] = xi - x4i;
443 xr = x1r + x1i;
444 xi = x1r - x1i;
445 a[2][j] = xr + x5r;
446 a[13][j] = xr - x5r;
447 a[10][j] = xi + x5i;
448 a[5][j] = xi - x5i;
449 xr = x2r + x2i;
450 xi = x2r - x2i;
451 a[4][j] = xr + x6r;
452 a[11][j] = xr - x6r;
453 a[12][j] = xi + x6i;
454 a[3][j] = xi - x6i;
455 xr = x3r + x3i;
456 xi = x3r - x3i;
457 a[6][j] = xr + x7r;
458 a[9][j] = xr - x7r;
459 a[14][j] = xi + x7i;
460 a[1][j] = xi - x7i;
461 }
462 for (j = 0; j <= 15; j++) {
463 x5r = C16_1R * a[j][1] + C16_1I * a[j][15];
464 x5i = C16_1R * a[j][15] - C16_1I * a[j][1];
465 xr = C16_7R * a[j][7] + C16_7I * a[j][9];
466 xi = C16_7R * a[j][9] - C16_7I * a[j][7];
467 x4r = x5r + xr;
468 x4i = x5i - xi;
469 x5r -= xr;
470 x5i += xi;
471 x7r = C16_5R * a[j][5] + C16_5I * a[j][11];
472 x7i = C16_5R * a[j][11] - C16_5I * a[j][5];
473 xr = C16_3R * a[j][3] + C16_3I * a[j][13];
474 xi = C16_3R * a[j][13] - C16_3I * a[j][3];
475 x6r = x7r + xr;
476 x6i = x7i - xi;
477 x7r -= xr;
478 x7i += xi;
479 xr = x4r - x6r;
480 xi = x4i - x6i;
481 x4r += x6r;
482 x4i += x6i;
483 x6r = W16_8R * (xi + xr);
484 x6i = W16_8R * (xi - xr);
485 xr = x5r + x7i;
486 xi = x5i - x7r;
487 x5r -= x7i;
488 x5i += x7r;
489 x7r = W16_4I * x5r + W16_4R * x5i;
490 x7i = W16_4I * x5i - W16_4R * x5r;
491 x5r = W16_4R * xr + W16_4I * xi;
492 x5i = W16_4R * xi - W16_4I * xr;
493 xr = C16_4R * a[j][4] + C16_4I * a[j][12];
494 xi = C16_4R * a[j][12] - C16_4I * a[j][4];
495 x2r = C16_8R * (a[j][0] + a[j][8]);
496 x3r = C16_8R * (a[j][0] - a[j][8]);
497 x0r = x2r + xr;
498 x1r = x3r + xi;
499 x2r -= xr;
500 x3r -= xi;
501 x0i = C16_2R * a[j][2] + C16_2I * a[j][14];
502 x2i = C16_2R * a[j][14] - C16_2I * a[j][2];
503 x1i = C16_6R * a[j][6] + C16_6I * a[j][10];
504 x3i = C16_6R * a[j][10] - C16_6I * a[j][6];
505 xr = x0i - x1i;
506 xi = x2i + x3i;
507 x0i += x1i;
508 x2i -= x3i;
509 x1i = W16_8R * (xi + xr);
510 x3i = W16_8R * (xi - xr);
511 xr = x0r + x0i;
512 xi = x0r - x0i;
513 a[j][0] = xr + x4r;
514 a[j][15] = xr - x4r;
515 a[j][8] = xi + x4i;
516 a[j][7] = xi - x4i;
517 xr = x1r + x1i;
518 xi = x1r - x1i;
519 a[j][2] = xr + x5r;
520 a[j][13] = xr - x5r;
521 a[j][10] = xi + x5i;
522 a[j][5] = xi - x5i;
523 xr = x2r + x2i;
524 xi = x2r - x2i;
525 a[j][4] = xr + x6r;
526 a[j][11] = xr - x6r;
527 a[j][12] = xi + x6i;
528 a[j][3] = xi - x6i;
529 xr = x3r + x3i;
530 xi = x3r - x3i;
531 a[j][6] = xr + x7r;
532 a[j][9] = xr - x7r;
533 a[j][14] = xi + x7i;
534 a[j][1] = xi - x7i;
535 }
536 }
537 }
538
539