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
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 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 // Generate a sine window.
52 void ff_sine_window_init(float *window
, int n
) {
54 for(i
= 0; i
< n
; i
++)
55 window
[i
] = sin((i
+ 0.5) / (2 * n
) * M_PI
);
59 * init MDCT or IMDCT computation.
61 int ff_mdct_init(MDCTContext
*s
, int nbits
, int inverse
)
66 memset(s
, 0, sizeof(*s
));
71 s
->tcos
= av_malloc(n4
* sizeof(FFTSample
));
74 s
->tsin
= av_malloc(n4
* sizeof(FFTSample
));
79 alpha
= 2 * M_PI
* (i
+ 1.0 / 8.0) / n
;
80 s
->tcos
[i
] = -cos(alpha
);
81 s
->tsin
[i
] = -sin(alpha
);
83 if (ff_fft_init(&s
->fft
, s
->nbits
- 2, inverse
) < 0)
92 /* complex multiplication: p = a * b */
93 #define CMUL(pre, pim, are, aim, bre, bim) \
99 (pre) = _are * _bre - _aim * _bim;\
100 (pim) = _are * _bim + _aim * _bre;\
103 static void imdct_c(MDCTContext
*s
, const FFTSample
*input
, FFTSample
*tmp
)
106 const uint16_t *revtab
= s
->fft
.revtab
;
107 const FFTSample
*tcos
= s
->tcos
;
108 const FFTSample
*tsin
= s
->tsin
;
109 const FFTSample
*in1
, *in2
;
110 FFTComplex
*z
= (FFTComplex
*)tmp
;
118 in2
= input
+ n2
- 1;
119 for(k
= 0; k
< n4
; k
++) {
121 CMUL(z
[j
].re
, z
[j
].im
, *in2
, *in1
, tcos
[k
], tsin
[k
]);
125 ff_fft_calc(&s
->fft
, z
);
127 /* post rotation + reordering */
129 for(k
= 0; k
< n4
; k
++) {
130 CMUL(z
[k
].re
, z
[k
].im
, z
[k
].re
, z
[k
].im
, tcos
[k
], tsin
[k
]);
135 * Compute inverse MDCT of size N = 2^nbits
136 * @param output N samples
137 * @param input N/2 samples
138 * @param tmp N/2 samples
140 void ff_imdct_calc(MDCTContext
*s
, FFTSample
*output
,
141 const FFTSample
*input
, FFTSample
*tmp
)
144 FFTComplex
*z
= (FFTComplex
*)tmp
;
149 imdct_c(s
, input
, tmp
);
151 for(k
= 0; k
< n8
; k
++) {
152 output
[2*k
] = -z
[n8
+ k
].im
;
153 output
[n2
-1-2*k
] = z
[n8
+ k
].im
;
155 output
[2*k
+1] = z
[n8
-1-k
].re
;
156 output
[n2
-1-2*k
-1] = -z
[n8
-1-k
].re
;
158 output
[n2
+ 2*k
]=-z
[k
+n8
].re
;
159 output
[n
-1- 2*k
]=-z
[k
+n8
].re
;
161 output
[n2
+ 2*k
+1]=z
[n8
-k
-1].im
;
162 output
[n
-2 - 2 * k
] = z
[n8
-k
-1].im
;
167 * Compute the middle half of the inverse MDCT of size N = 2^nbits,
168 * thus excluding the parts that can be derived by symmetry
169 * @param output N/2 samples
170 * @param input N/2 samples
171 * @param tmp N/2 samples
173 void ff_imdct_half(MDCTContext
*s
, FFTSample
*output
,
174 const FFTSample
*input
, FFTSample
*tmp
)
177 FFTComplex
*z
= (FFTComplex
*)tmp
;
182 imdct_c(s
, input
, tmp
);
184 for(k
= 0; k
< n8
; k
++) {
185 output
[n4
-1-2*k
] = z
[n8
+k
].im
;
186 output
[n4
-1-2*k
-1] = -z
[n8
-k
-1].re
;
187 output
[n4
+ 2*k
] = -z
[n8
+k
].re
;
188 output
[n4
+ 2*k
+1] = z
[n8
-k
-1].im
;
193 * Compute MDCT of size N = 2^nbits
194 * @param input N samples
195 * @param out N/2 samples
196 * @param tmp temporary storage of N/2 samples
198 void ff_mdct_calc(MDCTContext
*s
, FFTSample
*out
,
199 const FFTSample
*input
, FFTSample
*tmp
)
201 int i
, j
, n
, n8
, n4
, n2
, n3
;
202 FFTSample re
, im
, re1
, im1
;
203 const uint16_t *revtab
= s
->fft
.revtab
;
204 const FFTSample
*tcos
= s
->tcos
;
205 const FFTSample
*tsin
= s
->tsin
;
206 FFTComplex
*x
= (FFTComplex
*)tmp
;
216 re
= -input
[2*i
+3*n4
] - input
[n3
-1-2*i
];
217 im
= -input
[n4
+2*i
] + input
[n4
-1-2*i
];
219 CMUL(x
[j
].re
, x
[j
].im
, re
, im
, -tcos
[i
], tsin
[i
]);
221 re
= input
[2*i
] - input
[n2
-1-2*i
];
222 im
= -(input
[n2
+2*i
] + input
[n
-1-2*i
]);
224 CMUL(x
[j
].re
, x
[j
].im
, re
, im
, -tcos
[n8
+ i
], tsin
[n8
+ i
]);
227 ff_fft_calc(&s
->fft
, x
);
233 CMUL(re1
, im1
, re
, im
, -tsin
[i
], -tcos
[i
]);
239 void ff_mdct_end(MDCTContext
*s
)