1/* 2 * MDCT/IMDCT transforms 3 * Copyright (c) 2002 Fabrice Bellard 4 * 5 * This file is part of FFmpeg. 6 * 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. 11 * 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. 16 * 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 20 */ 21 22#include <stdlib.h> 23#include <string.h> 24#include "libavutil/common.h" 25#include "libavutil/mathematics.h" 26#include "fft.h" 27#include "fft-internal.h" 28 29/** 30 * @file 31 * MDCT/IMDCT transforms. 32 */ 33 34#if FFT_FLOAT 35# define RSCALE(x) (x) 36#else 37#if FFT_FIXED_32 38# define RSCALE(x) (((x) + 32) >> 6) 39#else /* FFT_FIXED_32 */ 40# define RSCALE(x) ((x) >> 1) 41#endif /* FFT_FIXED_32 */ 42#endif 43 44/** 45 * init MDCT or IMDCT computation. 46 */ 47av_cold int ff_mdct_init(FFTContext *s, int nbits, int inverse, double scale) 48{ 49 int n, n4, i; 50 double alpha, theta; 51 int tstep; 52 53 memset(s, 0, sizeof(*s)); 54 n = 1 << nbits; 55 s->mdct_bits = nbits; 56 s->mdct_size = n; 57 n4 = n >> 2; 58 s->mdct_permutation = FF_MDCT_PERM_NONE; 59 60 if (ff_fft_init(s, s->mdct_bits - 2, inverse) < 0) 61 goto fail; 62 63 s->tcos = av_malloc_array(n/2, sizeof(FFTSample)); 64 if (!s->tcos) 65 goto fail; 66 67 switch (s->mdct_permutation) { 68 case FF_MDCT_PERM_NONE: 69 s->tsin = s->tcos + n4; 70 tstep = 1; 71 break; 72 case FF_MDCT_PERM_INTERLEAVE: 73 s->tsin = s->tcos + 1; 74 tstep = 2; 75 break; 76 default: 77 goto fail; 78 } 79 80 theta = 1.0 / 8.0 + (scale < 0 ? n4 : 0); 81 scale = sqrt(fabs(scale)); 82 for(i=0;i<n4;i++) { 83 alpha = 2 * M_PI * (i + theta) / n; 84 s->tcos[i*tstep] = FIX15(-cos(alpha) * scale); 85 s->tsin[i*tstep] = FIX15(-sin(alpha) * scale); 86 } 87 return 0; 88 fail: 89 ff_mdct_end(s); 90 return -1; 91} 92 93/** 94 * Compute the middle half of the inverse MDCT of size N = 2^nbits, 95 * thus excluding the parts that can be derived by symmetry 96 * @param output N/2 samples 97 * @param input N/2 samples 98 */ 99void ff_imdct_half_c(FFTContext *s, FFTSample *output, const FFTSample *input) 100{ 101 int k, n8, n4, n2, n, j; 102 const uint16_t *revtab = s->revtab; 103 const FFTSample *tcos = s->tcos; 104 const FFTSample *tsin = s->tsin; 105 const FFTSample *in1, *in2; 106 FFTComplex *z = (FFTComplex *)output; 107 108 n = 1 << s->mdct_bits; 109 n2 = n >> 1; 110 n4 = n >> 2; 111 n8 = n >> 3; 112 113 /* pre rotation */ 114 in1 = input; 115 in2 = input + n2 - 1; 116 for(k = 0; k < n4; k++) { 117 j=revtab[k]; 118 CMUL(z[j].re, z[j].im, *in2, *in1, tcos[k], tsin[k]); 119 in1 += 2; 120 in2 -= 2; 121 } 122 s->fft_calc(s, z); 123 124 /* post rotation + reordering */ 125 for(k = 0; k < n8; k++) { 126 FFTSample r0, i0, r1, i1; 127 CMUL(r0, i1, z[n8-k-1].im, z[n8-k-1].re, tsin[n8-k-1], tcos[n8-k-1]); 128 CMUL(r1, i0, z[n8+k ].im, z[n8+k ].re, tsin[n8+k ], tcos[n8+k ]); 129 z[n8-k-1].re = r0; 130 z[n8-k-1].im = i0; 131 z[n8+k ].re = r1; 132 z[n8+k ].im = i1; 133 } 134} 135 136/** 137 * Compute inverse MDCT of size N = 2^nbits 138 * @param output N samples 139 * @param input N/2 samples 140 */ 141void ff_imdct_calc_c(FFTContext *s, FFTSample *output, const FFTSample *input) 142{ 143 int k; 144 int n = 1 << s->mdct_bits; 145 int n2 = n >> 1; 146 int n4 = n >> 2; 147 148 ff_imdct_half_c(s, output+n4, input); 149 150 for(k = 0; k < n4; k++) { 151 output[k] = -output[n2-k-1]; 152 output[n-k-1] = output[n2+k]; 153 } 154} 155 156/** 157 * Compute MDCT of size N = 2^nbits 158 * @param input N samples 159 * @param out N/2 samples 160 */ 161void ff_mdct_calc_c(FFTContext *s, FFTSample *out, const FFTSample *input) 162{ 163 int i, j, n, n8, n4, n2, n3; 164 FFTDouble re, im; 165 const uint16_t *revtab = s->revtab; 166 const FFTSample *tcos = s->tcos; 167 const FFTSample *tsin = s->tsin; 168 FFTComplex *x = (FFTComplex *)out; 169 170 n = 1 << s->mdct_bits; 171 n2 = n >> 1; 172 n4 = n >> 2; 173 n8 = n >> 3; 174 n3 = 3 * n4; 175 176 /* pre rotation */ 177 for(i=0;i<n8;i++) { 178 re = RSCALE(-input[2*i+n3] - input[n3-1-2*i]); 179 im = RSCALE(-input[n4+2*i] + input[n4-1-2*i]); 180 j = revtab[i]; 181 CMUL(x[j].re, x[j].im, re, im, -tcos[i], tsin[i]); 182 183 re = RSCALE( input[2*i] - input[n2-1-2*i]); 184 im = RSCALE(-input[n2+2*i] - input[ n-1-2*i]); 185 j = revtab[n8 + i]; 186 CMUL(x[j].re, x[j].im, re, im, -tcos[n8 + i], tsin[n8 + i]); 187 } 188 189 s->fft_calc(s, x); 190 191 /* post rotation */ 192 for(i=0;i<n8;i++) { 193 FFTSample r0, i0, r1, i1; 194 CMUL(i1, r0, x[n8-i-1].re, x[n8-i-1].im, -tsin[n8-i-1], -tcos[n8-i-1]); 195 CMUL(i0, r1, x[n8+i ].re, x[n8+i ].im, -tsin[n8+i ], -tcos[n8+i ]); 196 x[n8-i-1].re = r0; 197 x[n8-i-1].im = i0; 198 x[n8+i ].re = r1; 199 x[n8+i ].im = i1; 200 } 201} 202 203av_cold void ff_mdct_end(FFTContext *s) 204{ 205 av_freep(&s->tcos); 206 ff_fft_end(s); 207} 208