dct.c
1.81 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#include <math.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#include "dct.h"
/* prototypes */
void dct_init(void);
void ref_fdct(DCTELEM block[8][8]);
void ref_idct(DCTELEM block[8][8]);
/* Reference double-precision FDCT and IDCT */
/* The cosine lookup table */
/* coslu[a][b] = C(b)/2 * cos[(2a+1)b*pi/16] */
double coslu[8][8];
/* Routine to initialise the cosine lookup table */
void dct_init(void)
{
int a,b;
double tmp;
for(a=0;a<8;a++)
for(b=0;b<8;b++) {
tmp = cos((double)((a+a+1)*b) * (M_PI / 16.0));
if(b==0)
tmp /= sqrt(2.0);
coslu[a][b] = tmp * 0.5;
}
}
void ref_fdct (DCTELEM block[8][8])
{
int x,y,u,v;
double tmp, tmp2;
double res[8][8];
if( coslu[0][0] == 0.0 )
dct_init();
for (v=0; v<8; v++) {
for (u=0; u<8; u++) {
tmp = 0.0;
for (y=0; y<8; y++) {
tmp2 = 0.0;
for (x=0; x<8; x++) {
tmp2 += (double) block[y][x] * coslu[x][u];
}
tmp += coslu[y][v] * tmp2;
}
res[v][u] = tmp;
}
}
for (v=0; v<8; v++) {
for (u=0; u<8; u++) {
tmp = res[v][u];
if (tmp < 0.0) {
x = - ((int) (0.5 - tmp));
} else {
x = (int) (tmp + 0.5);
}
block[v][u] = (DCTELEM) x;
}
}
}
void ref_idct (DCTELEM block[8][8])
{
int x,y,u,v;
double tmp, tmp2;
double res[8][8];
if( coslu[0][0] == 0.0 )
dct_init();
for (y=0; y<8; y++) {
for (x=0; x<8; x++) {
tmp = 0.0;
for (v=0; v<8; v++) {
tmp2 = 0.0;
for (u=0; u<8; u++) {
tmp2 += (double) block[v][u] * coslu[x][u];
}
tmp += coslu[y][v] * tmp2;
}
res[y][x] = tmp;
}
}
for (v=0; v<8; v++) {
for (u=0; u<8; u++) {
tmp = res[v][u];
if (tmp < 0.0) {
x = - ((int) (0.5 - tmp));
} else {
x = (int) (tmp + 0.5);
}
block[v][u] = (DCTELEM) x;
}
}
}