ffmpeg / libavcodec / fdctref.c @ bb7d4939
History  View  Annotate  Download (3.87 KB)
1 
/* fdctref.c, forward discrete cosine transform, double precision */


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

4  
5 
/*

6 
* Disclaimer of Warranty

7 
*

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

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

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

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

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

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

14 
* arising from the use of these programs.

15 
*

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

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

18 
*

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

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

21 
* patents.

22 
*

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

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

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

26 
* design.

27 
*

28 
*/

29  
30 
#include <math.h> 
31  
32 
#ifndef PI

33 
# ifdef M_PI

34 
# define PI M_PI

35 
# else

36 
# define PI 3.14159265358979323846 
37 
# endif

38 
#endif

39  
40 
/* global declarations */

41 
void init_fdct (void); 
42 
void fdct (short *block); 
43  
44 
/* private data */

45 
static double c[8][8]; /* transform coefficients */ 
46  
47 
void init_fdct()

48 
{ 
49 
int i, j;

50 
double s;

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

62 
short *block;

63 
{ 
64 
register int i, j; 
65 
double s;

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

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

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

76 
*/

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

87 
} 
88  
89 
for(j = 0; j < 8; j++) 
90 
for(i = 0; i < 8; i++) 
91 
{ 
92 
s = 0.0; 
93  
94 
/*

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

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

97 
*/

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

109 
* reason for adding 0.499999 instead of 0.5:

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

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

112 
* extremely high arithmetic implementation dependency of the result;

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

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

115 
* (if at all)

116 
*/

117 
} 
118 
} 
119  
120 
/* perform IDCT matrix multiply for 8x8 coefficient block */

121  
122 
void idct(block)

123 
short *block;

124 
{ 
125 
int i, j, k, v;

126 
double partial_product;

127 
double tmp[64]; 
128  
129 
for (i=0; i<8; i++) 
130 
for (j=0; j<8; j++) 
131 
{ 
132 
partial_product = 0.0; 
133  
134 
for (k=0; k<8; k++) 
135 
partial_product+= c[k][j]*block[8*i+k];

136  
137 
tmp[8*i+j] = partial_product;

138 
} 
139  
140 
/* Transpose operation is integrated into address mapping by switching

141 
loop order of i and j */

142  
143 
for (j=0; j<8; j++) 
144 
for (i=0; i<8; i++) 
145 
{ 
146 
partial_product = 0.0; 
147  
148 
for (k=0; k<8; k++) 
149 
partial_product+= c[k][i]*tmp[8*k+j];

150  
151 
v = (int) floor(partial_product+0.5); 
152 
block[8*i+j] = v;

153 
} 
154 
} 