## ffmpeg / libavcodec / fdctref.c @ bb7d4939

History | View | Annotate | Download (3.87 KB)

1 | de6d9b64 | Fabrice Bellard | ```
/* 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 copyright-holder 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 third-party
``` |
||

21 | ```
* patents.
``` |
||

22 | ```
*
``` |
||

23 | ```
* Commercial implementations of MPEG-1 and MPEG-2 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 | dc541ee7 | Fabrice Bellard | |

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 | } |