ffmpeg / libavcodec / mdct.c @ 94d85eaf
History  View  Annotate  Download (4.43 KB)
1 
/*


2 
* MDCT/IMDCT transforms

3 
* Copyright (c) 2002 Fabrice Bellard.

4 
*

5 
* This library is free software; you can redistribute it and/or

6 
* modify it under the terms of the GNU Lesser General Public

7 
* License as published by the Free Software Foundation; either

8 
* version 2 of the License, or (at your option) any later version.

9 
*

10 
* This library is distributed in the hope that it will be useful,

11 
* but WITHOUT ANY WARRANTY; without even the implied warranty of

12 
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU

13 
* Lesser General Public License for more details.

14 
*

15 
* You should have received a copy of the GNU Lesser General Public

16 
* License along with this library; if not, write to the Free Software

17 
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 021101301 USA

18 
*/

19 
#include "dsputil.h" 
20  
21 
/**

22 
* @file mdct.c

23 
* MDCT/IMDCT transforms.

24 
*/

25  
26 
/**

27 
* init MDCT or IMDCT computation.

28 
*/

29 
int ff_mdct_init(MDCTContext *s, int nbits, int inverse) 
30 
{ 
31 
int n, n4, i;

32 
float alpha;

33  
34 
memset(s, 0, sizeof(*s)); 
35 
n = 1 << nbits;

36 
s>nbits = nbits; 
37 
s>n = n; 
38 
n4 = n >> 2;

39 
s>tcos = av_malloc(n4 * sizeof(FFTSample));

40 
if (!s>tcos)

41 
goto fail;

42 
s>tsin = av_malloc(n4 * sizeof(FFTSample));

43 
if (!s>tsin)

44 
goto fail;

45  
46 
for(i=0;i<n4;i++) { 
47 
alpha = 2 * M_PI * (i + 1.0 / 8.0) / n; 
48 
s>tcos[i] = cos(alpha); 
49 
s>tsin[i] = sin(alpha); 
50 
} 
51 
if (ff_fft_init(&s>fft, s>nbits  2, inverse) < 0) 
52 
goto fail;

53 
return 0; 
54 
fail:

55 
av_freep(&s>tcos); 
56 
av_freep(&s>tsin); 
57 
return 1; 
58 
} 
59  
60 
/* complex multiplication: p = a * b */

61 
#define CMUL(pre, pim, are, aim, bre, bim) \

62 
{\ 
63 
float _are = (are);\

64 
float _aim = (aim);\

65 
float _bre = (bre);\

66 
float _bim = (bim);\

67 
(pre) = _are * _bre  _aim * _bim;\ 
68 
(pim) = _are * _bim + _aim * _bre;\ 
69 
} 
70  
71 
/**

72 
* Compute inverse MDCT of size N = 2^nbits

73 
* @param output N samples

74 
* @param input N/2 samples

75 
* @param tmp N/2 samples

76 
*/

77 
void ff_imdct_calc(MDCTContext *s, FFTSample *output,

78 
const FFTSample *input, FFTSample *tmp)

79 
{ 
80 
int k, n8, n4, n2, n, j;

81 
const uint16_t *revtab = s>fft.revtab;

82 
const FFTSample *tcos = s>tcos;

83 
const FFTSample *tsin = s>tsin;

84 
const FFTSample *in1, *in2;

85 
FFTComplex *z = (FFTComplex *)tmp; 
86  
87 
n = 1 << s>nbits;

88 
n2 = n >> 1;

89 
n4 = n >> 2;

90 
n8 = n >> 3;

91  
92 
/* pre rotation */

93 
in1 = input; 
94 
in2 = input + n2  1;

95 
for(k = 0; k < n4; k++) { 
96 
j=revtab[k]; 
97 
CMUL(z[j].re, z[j].im, *in2, *in1, tcos[k], tsin[k]); 
98 
in1 += 2;

99 
in2 = 2;

100 
} 
101 
ff_fft_calc(&s>fft, z); 
102  
103 
/* post rotation + reordering */

104 
/* XXX: optimize */

105 
for(k = 0; k < n4; k++) { 
106 
CMUL(z[k].re, z[k].im, z[k].re, z[k].im, tcos[k], tsin[k]); 
107 
} 
108 
for(k = 0; k < n8; k++) { 
109 
output[2*k] = z[n8 + k].im;

110 
output[n212*k] = z[n8 + k].im; 
111  
112 
output[2*k+1] = z[n81k].re; 
113 
output[n212*k1] = z[n81k].re; 
114  
115 
output[n2 + 2*k]=z[k+n8].re;

116 
output[n1 2*k]=z[k+n8].re; 
117  
118 
output[n2 + 2*k+1]=z[n8k1].im; 
119 
output[n2  2 * k] = z[n8k1].im; 
120 
} 
121 
} 
122  
123 
/**

124 
* Compute MDCT of size N = 2^nbits

125 
* @param input N samples

126 
* @param out N/2 samples

127 
* @param tmp temporary storage of N/2 samples

128 
*/

129 
void ff_mdct_calc(MDCTContext *s, FFTSample *out,

130 
const FFTSample *input, FFTSample *tmp)

131 
{ 
132 
int i, j, n, n8, n4, n2, n3;

133 
FFTSample re, im, re1, im1; 
134 
const uint16_t *revtab = s>fft.revtab;

135 
const FFTSample *tcos = s>tcos;

136 
const FFTSample *tsin = s>tsin;

137 
FFTComplex *x = (FFTComplex *)tmp; 
138  
139 
n = 1 << s>nbits;

140 
n2 = n >> 1;

141 
n4 = n >> 2;

142 
n8 = n >> 3;

143 
n3 = 3 * n4;

144  
145 
/* pre rotation */

146 
for(i=0;i<n8;i++) { 
147 
re = input[2*i+3*n4]  input[n312*i]; 
148 
im = input[n4+2*i] + input[n412*i]; 
149 
j = revtab[i]; 
150 
CMUL(x[j].re, x[j].im, re, im, tcos[i], tsin[i]); 
151  
152 
re = input[2*i]  input[n212*i]; 
153 
im = (input[n2+2*i] + input[n12*i]); 
154 
j = revtab[n8 + i]; 
155 
CMUL(x[j].re, x[j].im, re, im, tcos[n8 + i], tsin[n8 + i]); 
156 
} 
157  
158 
ff_fft_calc(&s>fft, x); 
159  
160 
/* post rotation */

161 
for(i=0;i<n4;i++) { 
162 
re = x[i].re; 
163 
im = x[i].im; 
164 
CMUL(re1, im1, re, im, tsin[i], tcos[i]); 
165 
out[2*i] = im1;

166 
out[n212*i] = re1; 
167 
} 
168 
} 
169  
170 
void ff_mdct_end(MDCTContext *s)

171 
{ 
172 
av_freep(&s>tcos); 
173 
av_freep(&s>tsin); 
174 
ff_fft_end(&s>fft); 
175 
} 