Jamoma API  0.6.0.a19
shrtdct.cpp
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 
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 
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