2 * MDCT/IMDCT transforms
3 * Copyright (c) 2002 Fabrice Bellard
5 * This file is part of Libav.
7 * Libav 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 * Libav 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 Libav; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 #include "libavutil/common.h"
25 #include "libavutil/mathematics.h"
27 #include "fft-internal.h"
31 * MDCT/IMDCT transforms.
35 # define RSCALE(x) (x)
37 # define RSCALE(x) ((x) >> 1)
41 * init MDCT or IMDCT computation.
43 av_cold
int ff_mdct_init(FFTContext
*s
, int nbits
, int inverse
, double scale
)
49 memset(s
, 0, sizeof(*s
));
54 s
->mdct_permutation
= FF_MDCT_PERM_NONE
;
56 if (ff_fft_init(s
, s
->mdct_bits
- 2, inverse
) < 0)
59 s
->tcos
= av_malloc(n
/2 * sizeof(FFTSample
));
63 switch (s
->mdct_permutation
) {
64 case FF_MDCT_PERM_NONE
:
65 s
->tsin
= s
->tcos
+ n4
;
68 case FF_MDCT_PERM_INTERLEAVE
:
69 s
->tsin
= s
->tcos
+ 1;
76 theta
= 1.0 / 8.0 + (scale
< 0 ? n4
: 0);
77 scale
= sqrt(fabs(scale
));
79 alpha
= 2 * M_PI
* (i
+ theta
) / n
;
80 s
->tcos
[i
*tstep
] = FIX15(-cos(alpha
) * scale
);
81 s
->tsin
[i
*tstep
] = FIX15(-sin(alpha
) * scale
);
90 * Compute the middle half of the inverse MDCT of size N = 2^nbits,
91 * thus excluding the parts that can be derived by symmetry
92 * @param output N/2 samples
93 * @param input N/2 samples
95 void ff_imdct_half_c(FFTContext
*s
, FFTSample
*output
, const FFTSample
*input
)
97 int k
, n8
, n4
, n2
, n
, j
;
98 const uint16_t *revtab
= s
->revtab
;
99 const FFTSample
*tcos
= s
->tcos
;
100 const FFTSample
*tsin
= s
->tsin
;
101 const FFTSample
*in1
, *in2
;
102 FFTComplex
*z
= (FFTComplex
*)output
;
104 n
= 1 << s
->mdct_bits
;
111 in2
= input
+ n2
- 1;
112 for(k
= 0; k
< n4
; k
++) {
114 CMUL(z
[j
].re
, z
[j
].im
, *in2
, *in1
, tcos
[k
], tsin
[k
]);
120 /* post rotation + reordering */
121 for(k
= 0; k
< n8
; k
++) {
122 FFTSample r0
, i0
, r1
, i1
;
123 CMUL(r0
, i1
, z
[n8
-k
-1].im
, z
[n8
-k
-1].re
, tsin
[n8
-k
-1], tcos
[n8
-k
-1]);
124 CMUL(r1
, i0
, z
[n8
+k
].im
, z
[n8
+k
].re
, tsin
[n8
+k
], tcos
[n8
+k
]);
133 * Compute inverse MDCT of size N = 2^nbits
134 * @param output N samples
135 * @param input N/2 samples
137 void ff_imdct_calc_c(FFTContext
*s
, FFTSample
*output
, const FFTSample
*input
)
140 int n
= 1 << s
->mdct_bits
;
144 ff_imdct_half_c(s
, output
+n4
, input
);
146 for(k
= 0; k
< n4
; k
++) {
147 output
[k
] = -output
[n2
-k
-1];
148 output
[n
-k
-1] = output
[n2
+k
];
153 * Compute MDCT of size N = 2^nbits
154 * @param input N samples
155 * @param out N/2 samples
157 void ff_mdct_calc_c(FFTContext
*s
, FFTSample
*out
, const FFTSample
*input
)
159 int i
, j
, n
, n8
, n4
, n2
, n3
;
161 const uint16_t *revtab
= s
->revtab
;
162 const FFTSample
*tcos
= s
->tcos
;
163 const FFTSample
*tsin
= s
->tsin
;
164 FFTComplex
*x
= (FFTComplex
*)out
;
166 n
= 1 << s
->mdct_bits
;
174 re
= RSCALE(-input
[2*i
+n3
] - input
[n3
-1-2*i
]);
175 im
= RSCALE(-input
[n4
+2*i
] + input
[n4
-1-2*i
]);
177 CMUL(x
[j
].re
, x
[j
].im
, re
, im
, -tcos
[i
], tsin
[i
]);
179 re
= RSCALE( input
[2*i
] - input
[n2
-1-2*i
]);
180 im
= RSCALE(-input
[n2
+2*i
] - input
[ n
-1-2*i
]);
182 CMUL(x
[j
].re
, x
[j
].im
, re
, im
, -tcos
[n8
+ i
], tsin
[n8
+ i
]);
189 FFTSample r0
, i0
, r1
, i1
;
190 CMUL(i1
, r0
, x
[n8
-i
-1].re
, x
[n8
-i
-1].im
, -tsin
[n8
-i
-1], -tcos
[n8
-i
-1]);
191 CMUL(i0
, r1
, x
[n8
+i
].re
, x
[n8
+i
].im
, -tsin
[n8
+i
], -tcos
[n8
+i
]);
199 av_cold
void ff_mdct_end(FFTContext
*s
)