ffmpeg / libavcodec / fdctref.c @ 186447f8
History  View  Annotate  Download (3.89 KB)
1 
/**


2 
* @file fdctref.c

3 
* forward discrete cosine transform, double precision.

4 
*/

5  
6 
/* Copyright (C) 1996, MPEG Software Simulation Group. All Rights Reserved. */

7  
8 
/*

9 
* Disclaimer of Warranty

10 
*

11 
* These software programs are available to the user without any license fee or

12 
* royalty on an "as is" basis. The MPEG Software Simulation Group disclaims

13 
* any and all warranties, whether express, implied, or statuary, including any

14 
* implied warranties or merchantability or of fitness for a particular

15 
* purpose. In no event shall the copyrightholder be liable for any

16 
* incidental, punitive, or consequential damages of any kind whatsoever

17 
* arising from the use of these programs.

18 
*

19 
* This disclaimer of warranty extends to the user of these programs and user's

20 
* customers, employees, agents, transferees, successors, and assigns.

21 
*

22 
* The MPEG Software Simulation Group does not represent or warrant that the

23 
* programs furnished hereunder are free of infringement of any thirdparty

24 
* patents.

25 
*

26 
* Commercial implementations of MPEG1 and MPEG2 video, including shareware,

27 
* are subject to royalty fees to patent holders. Many of these patents are

28 
* general enough such that they are unavoidable regardless of implementation

29 
* design.

30 
*

31 
*/

32  
33 
#include <math.h> 
34  
35 
#ifndef PI

36 
# ifdef M_PI

37 
# define PI M_PI

38 
# else

39 
# define PI 3.14159265358979323846 
40 
# endif

41 
#endif

42  
43 
/* global declarations */

44 
void init_fdct (void); 
45 
void fdct (short *block); 
46  
47 
/* private data */

48 
static double c[8][8]; /* transform coefficients */ 
49  
50 
void init_fdct()

51 
{ 
52 
int i, j;

53 
double s;

54  
55 
for (i=0; i<8; i++) 
56 
{ 
57 
s = (i==0) ? sqrt(0.125) : 0.5; 
58  
59 
for (j=0; j<8; j++) 
60 
c[i][j] = s * cos((PI/8.0)*i*(j+0.5)); 
61 
} 
62 
} 
63  
64 
void fdct(block)

65 
short *block;

66 
{ 
67 
register int i, j; 
68 
double s;

69 
double tmp[64]; 
70  
71 
for(i = 0; i < 8; i++) 
72 
for(j = 0; j < 8; j++) 
73 
{ 
74 
s = 0.0; 
75  
76 
/*

77 
* for(k = 0; k < 8; k++)

78 
* s += c[j][k] * block[8 * i + k];

79 
*/

80 
s += c[j][0] * block[8 * i + 0]; 
81 
s += c[j][1] * block[8 * i + 1]; 
82 
s += c[j][2] * block[8 * i + 2]; 
83 
s += c[j][3] * block[8 * i + 3]; 
84 
s += c[j][4] * block[8 * i + 4]; 
85 
s += c[j][5] * block[8 * i + 5]; 
86 
s += c[j][6] * block[8 * i + 6]; 
87 
s += c[j][7] * block[8 * i + 7]; 
88  
89 
tmp[8 * i + j] = s;

90 
} 
91  
92 
for(j = 0; j < 8; j++) 
93 
for(i = 0; i < 8; i++) 
94 
{ 
95 
s = 0.0; 
96  
97 
/*

98 
* for(k = 0; k < 8; k++)

99 
* s += c[i][k] * tmp[8 * k + j];

100 
*/

101 
s += c[i][0] * tmp[8 * 0 + j]; 
102 
s += c[i][1] * tmp[8 * 1 + j]; 
103 
s += c[i][2] * tmp[8 * 2 + j]; 
104 
s += c[i][3] * tmp[8 * 3 + j]; 
105 
s += c[i][4] * tmp[8 * 4 + j]; 
106 
s += c[i][5] * tmp[8 * 5 + j]; 
107 
s += c[i][6] * tmp[8 * 6 + j]; 
108 
s += c[i][7] * tmp[8 * 7 + j]; 
109 
s*=8.0; 
110  
111 
block[8 * i + j] = (short)floor(s + 0.499999); 
112 
/*

113 
* reason for adding 0.499999 instead of 0.5:

114 
* s is quite often x.5 (at least for i and/or j = 0 or 4)

115 
* and setting the rounding threshold exactly to 0.5 leads to an

116 
* extremely high arithmetic implementation dependency of the result;

117 
* s being between x.5 and x.500001 (which is now incorrectly rounded

118 
* downwards instead of upwards) is assumed to occur less often

119 
* (if at all)

120 
*/

121 
} 
122 
} 
123  
124 
/* perform IDCT matrix multiply for 8x8 coefficient block */

125  
126 
void idct(block)

127 
short *block;

128 
{ 
129 
int i, j, k, v;

130 
double partial_product;

131 
double tmp[64]; 
132  
133 
for (i=0; i<8; i++) 
134 
for (j=0; j<8; j++) 
135 
{ 
136 
partial_product = 0.0; 
137  
138 
for (k=0; k<8; k++) 
139 
partial_product+= c[k][j]*block[8*i+k];

140  
141 
tmp[8*i+j] = partial_product;

142 
} 
143  
144 
/* Transpose operation is integrated into address mapping by switching

145 
loop order of i and j */

146  
147 
for (j=0; j<8; j++) 
148 
for (i=0; i<8; i++) 
149 
{ 
150 
partial_product = 0.0; 
151  
152 
for (k=0; k<8; k++) 
153 
partial_product+= c[k][i]*tmp[8*k+j];

154  
155 
v = (int) floor(partial_product+0.5); 
156 
block[8*i+j] = v;

157 
} 
158 
} 