2 * MDCT/IMDCT transforms
3 * Copyright (c) 2002 Fabrice Bellard
5 * This file is part of FFmpeg.
7 * FFmpeg is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
12 * FFmpeg is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with FFmpeg; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 * @file libavcodec/mdct.c
25 * MDCT/IMDCT transforms.
28 // Generate a Kaiser-Bessel Derived Window.
29 #define BESSEL_I0_ITER 50 // default: 50 iterations of Bessel I0 approximation
30 av_cold
void ff_kbd_window_init(float *window
, float alpha
, int n
)
33 double sum
= 0.0, bessel
, tmp
;
34 double local_window
[n
];
35 double alpha2
= (alpha
* M_PI
/ n
) * (alpha
* M_PI
/ n
);
37 for (i
= 0; i
< n
; i
++) {
38 tmp
= i
* (n
- i
) * alpha2
;
40 for (j
= BESSEL_I0_ITER
; j
> 0; j
--)
41 bessel
= bessel
* tmp
/ (j
* j
) + 1;
43 local_window
[i
] = sum
;
47 for (i
= 0; i
< n
; i
++)
48 window
[i
] = sqrt(local_window
[i
] / sum
);
51 DECLARE_ALIGNED(16, float, ff_sine_32
[ 32]);
52 DECLARE_ALIGNED(16, float, ff_sine_64
[ 64]);
53 DECLARE_ALIGNED(16, float, ff_sine_128
[ 128]);
54 DECLARE_ALIGNED(16, float, ff_sine_256
[ 256]);
55 DECLARE_ALIGNED(16, float, ff_sine_512
[ 512]);
56 DECLARE_ALIGNED(16, float, ff_sine_1024
[1024]);
57 DECLARE_ALIGNED(16, float, ff_sine_2048
[2048]);
58 DECLARE_ALIGNED(16, float, ff_sine_4096
[4096]);
59 float * const ff_sine_windows
[] = {
60 NULL
, NULL
, NULL
, NULL
, NULL
, // unused
61 ff_sine_32
, ff_sine_64
,
62 ff_sine_128
, ff_sine_256
, ff_sine_512
, ff_sine_1024
, ff_sine_2048
, ff_sine_4096
65 // Generate a sine window.
66 av_cold
void ff_sine_window_init(float *window
, int n
) {
68 for(i
= 0; i
< n
; i
++)
69 window
[i
] = sinf((i
+ 0.5) * (M_PI
/ (2.0 * n
)));
73 * init MDCT or IMDCT computation.
75 av_cold
int ff_mdct_init(FFTContext
*s
, int nbits
, int inverse
, double scale
)
81 memset(s
, 0, sizeof(*s
));
86 s
->permutation
= FF_MDCT_PERM_NONE
;
88 if (ff_fft_init(s
, s
->mdct_bits
- 2, inverse
) < 0)
91 s
->tcos
= av_malloc(n
/2 * sizeof(FFTSample
));
95 switch (s
->permutation
) {
96 case FF_MDCT_PERM_NONE
:
97 s
->tsin
= s
->tcos
+ n4
;
100 case FF_MDCT_PERM_INTERLEAVE
:
101 s
->tsin
= s
->tcos
+ 1;
108 theta
= 1.0 / 8.0 + (scale
< 0 ? n4
: 0);
109 scale
= sqrt(fabs(scale
));
111 alpha
= 2 * M_PI
* (i
+ theta
) / n
;
112 s
->tcos
[i
*tstep
] = -cos(alpha
) * scale
;
113 s
->tsin
[i
*tstep
] = -sin(alpha
) * scale
;
121 /* complex multiplication: p = a * b */
122 #define CMUL(pre, pim, are, aim, bre, bim) \
124 FFTSample _are = (are);\
125 FFTSample _aim = (aim);\
126 FFTSample _bre = (bre);\
127 FFTSample _bim = (bim);\
128 (pre) = _are * _bre - _aim * _bim;\
129 (pim) = _are * _bim + _aim * _bre;\
133 * Compute the middle half of the inverse MDCT of size N = 2^nbits,
134 * thus excluding the parts that can be derived by symmetry
135 * @param output N/2 samples
136 * @param input N/2 samples
138 void ff_imdct_half_c(FFTContext
*s
, FFTSample
*output
, const FFTSample
*input
)
140 int k
, n8
, n4
, n2
, n
, j
;
141 const uint16_t *revtab
= s
->revtab
;
142 const FFTSample
*tcos
= s
->tcos
;
143 const FFTSample
*tsin
= s
->tsin
;
144 const FFTSample
*in1
, *in2
;
145 FFTComplex
*z
= (FFTComplex
*)output
;
147 n
= 1 << s
->mdct_bits
;
154 in2
= input
+ n2
- 1;
155 for(k
= 0; k
< n4
; k
++) {
157 CMUL(z
[j
].re
, z
[j
].im
, *in2
, *in1
, tcos
[k
], tsin
[k
]);
163 /* post rotation + reordering */
164 for(k
= 0; k
< n8
; k
++) {
165 FFTSample r0
, i0
, r1
, i1
;
166 CMUL(r0
, i1
, z
[n8
-k
-1].im
, z
[n8
-k
-1].re
, tsin
[n8
-k
-1], tcos
[n8
-k
-1]);
167 CMUL(r1
, i0
, z
[n8
+k
].im
, z
[n8
+k
].re
, tsin
[n8
+k
], tcos
[n8
+k
]);
176 * Compute inverse MDCT of size N = 2^nbits
177 * @param output N samples
178 * @param input N/2 samples
180 void ff_imdct_calc_c(FFTContext
*s
, FFTSample
*output
, const FFTSample
*input
)
183 int n
= 1 << s
->mdct_bits
;
187 ff_imdct_half_c(s
, output
+n4
, input
);
189 for(k
= 0; k
< n4
; k
++) {
190 output
[k
] = -output
[n2
-k
-1];
191 output
[n
-k
-1] = output
[n2
+k
];
196 * Compute MDCT of size N = 2^nbits
197 * @param input N samples
198 * @param out N/2 samples
200 void ff_mdct_calc_c(FFTContext
*s
, FFTSample
*out
, const FFTSample
*input
)
202 int i
, j
, n
, n8
, n4
, n2
, n3
;
204 const uint16_t *revtab
= s
->revtab
;
205 const FFTSample
*tcos
= s
->tcos
;
206 const FFTSample
*tsin
= s
->tsin
;
207 FFTComplex
*x
= (FFTComplex
*)out
;
209 n
= 1 << s
->mdct_bits
;
217 re
= -input
[2*i
+3*n4
] - input
[n3
-1-2*i
];
218 im
= -input
[n4
+2*i
] + input
[n4
-1-2*i
];
220 CMUL(x
[j
].re
, x
[j
].im
, re
, im
, -tcos
[i
], tsin
[i
]);
222 re
= input
[2*i
] - input
[n2
-1-2*i
];
223 im
= -(input
[n2
+2*i
] + input
[n
-1-2*i
]);
225 CMUL(x
[j
].re
, x
[j
].im
, re
, im
, -tcos
[n8
+ i
], tsin
[n8
+ i
]);
232 FFTSample r0
, i0
, r1
, i1
;
233 CMUL(i1
, r0
, x
[n8
-i
-1].re
, x
[n8
-i
-1].im
, -tsin
[n8
-i
-1], -tcos
[n8
-i
-1]);
234 CMUL(i0
, r1
, x
[n8
+i
].re
, x
[n8
+i
].im
, -tsin
[n8
+i
], -tcos
[n8
+i
]);
242 av_cold
void ff_mdct_end(FFTContext
*s
)