git.sesse.net Git - vlc/blob - modules/codec/wmafixed/mdct.c

   1 /*
   2  * WMA compatible decoder
   3  * Copyright (c) 2002 The FFmpeg Project.
   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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  18  */
  19
  20 #include <string.h>
  21 #include "wmafixed.h"
  22 #include "mdct.h"
  23
  24 /*these are the sin and cos rotations used by the MDCT*/
  25
  26 /*accessed too infrequently to give much speedup in IRAM*/
  27
  28 int32_t *tcosarray[5], *tsinarray[5];
  29 int32_t tcos0[1024], tcos1[512], tcos2[256], tcos3[128], tcos4[64];
  30 int32_t tsin0[1024], tsin1[512], tsin2[256], tsin3[128], tsin4[64];
  31
  32 uint16_t revtab0[1024];
  33
  34 /**
  35  * init MDCT or IMDCT computation.
  36  */
  37 int ff_mdct_init(MDCTContext *s, int nbits, int inverse)
  38 {
  39     int n, n4, i;
  40
  41     memset(s, 0, sizeof(*s));
  42     n = 1 << nbits;            /* nbits ranges from 12 to 8 inclusive */
  43     s->nbits = nbits;
  44     s->n = n;
  45     n4 = n >> 2;
  46     s->tcos = tcosarray[12-nbits];
  47     s->tsin = tsinarray[12-nbits];
  48     for(i=0;i<n4;i++)
  49     {
  50         int32_t ip = itofix32(i) + 0x2000;
  51         ip = ip >> nbits;
  52
  53         /*I can't remember why this works, but it seems
  54           to agree for ~24 bits, maybe more!*/
  55         s->tsin[i] = - fsincos(ip<<16, &(s->tcos[i]));
  56         s->tcos[i] *=-1;
  57     }
  58
  59     (&s->fft)->nbits = nbits-2;
  60     (&s->fft)->inverse = inverse;
  61
  62     return 0;
  63
  64 }
  65
  66 /**
  67  * Compute inverse MDCT of size N = 2^nbits
  68  * @param output N samples
  69  * @param input N/2 samples
  70  * @param tmp N/2 samples
  71  */
  72 void ff_imdct_calc(MDCTContext *s,
  73                    int32_t *output,
  74                    int32_t *input)
  75 {
  76     int k, n8, n4, n2, n, j,scale;
  77     const int32_t *tcos = s->tcos;
  78     const int32_t *tsin = s->tsin;
  79     const int32_t *in1, *in2;
  80     FFTComplex *z1 = (FFTComplex *)output;
  81     FFTComplex *z2 = (FFTComplex *)input;
  82     int revtabshift = 12 - s->nbits;
  83
  84     n = 1 << s->nbits;
  85
  86     n2 = n >> 1;
  87     n4 = n >> 2;
  88     n8 = n >> 3;
  89
  90     /* pre rotation */
  91     in1 = input;
  92     in2 = input + n2 - 1;
  93
  94     for(k = 0; k < n4; k++)
  95     {
  96         j=revtab0[k<<revtabshift];
  97         CMUL(&z1[j].re, &z1[j].im, *in2, *in1, tcos[k], tsin[k]);
  98         in1 += 2;
  99         in2 -= 2;
 100     }
 101
 102     scale = fft_calc_unscaled(&s->fft, z1);
 103
 104     /* post rotation + reordering */
 105     for(k = 0; k < n4; k++)
 106     {
 107         CMUL(&z2[k].re, &z2[k].im, (z1[k].re), (z1[k].im), tcos[k], tsin[k]);
 108     }
 109
 110     for(k = 0; k < n8; k++)
 111     {
 112         int32_t r1,r2,r3,r4,r1n,r2n,r3n;
 113
 114         r1 = z2[n8 + k].im;
 115         r1n = r1 * -1;
 116         r2 = z2[n8-1-k].re;
 117         r2n = r2 * -1;
 118         r3 = z2[k+n8].re;
 119         r3n = r3 * -1;
 120         r4 = z2[n8-k-1].im;
 121
 122         output[2*k] = r1n;
 123         output[n2-1-2*k] = r1;
 124
 125         output[2*k+1] = r2;
 126         output[n2-1-2*k-1] = r2n;
 127
 128         output[n2 + 2*k]= r3n;
 129         output[n-1- 2*k]= r3n;
 130
 131         output[n2 + 2*k+1]= r4;
 132         output[n-2 - 2 * k] = r4;
 133     }
 134 }
 135
 136 /* init MDCT */
 137
 138 int mdct_init_global(void)
 139 {
 140     int i,j,m;
 141
 142     /* although seemingly degenerate, these cannot actually be merged together without
 143        a substantial increase in error which is unjustified by the tiny memory savings*/
 144
 145     tcosarray[0] = tcos0; tcosarray[1] = tcos1; tcosarray[2] = tcos2; tcosarray[3] = tcos3;tcosarray[4] = tcos4;
 146     tsinarray[0] = tsin0; tsinarray[1] = tsin1; tsinarray[2] = tsin2; tsinarray[3] = tsin3;tsinarray[4] = tsin4;
 147
 148     /* init the MDCT bit reverse table here rather then in fft_init */
 149
 150     for(i=0;i<1024;i++)           /*hard coded to a 2048 bit rotation*/
 151     {                             /*smaller sizes can reuse the largest*/
 152         m=0;
 153         for(j=0;j<10;j++)
 154         {
 155             m |= ((i >> j) & 1) << (10-j-1);
 156         }
 157
 158        revtab0[i]=m;
 159     }
 160
 161     fft_init_global();
 162
 163     return 0;
 164 }