apps/codecs/libwma/fft.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 "wmadec.h"
  21 #include "wmafixed.h"
  22
  23 FFTComplex  exptab0[512] IBSS_ATTR;
  24
  25 /* butter fly op */
  26 #define BF(pre, pim, qre, qim, pre1, pim1, qre1, qim1) \
  27 {\
  28   fixed32 ax, ay, bx, by;\
  29   bx=pre1;\
  30   by=pim1;\
  31   ax=qre1;\
  32   ay=qim1;\
  33   pre = (bx + ax);\
  34   pim = (by + ay);\
  35   qre = (bx - ax);\
  36   qim = (by - ay);\
  37 }
  38
  39
  40 int fft_calc_unscaled(FFTContext *s, FFTComplex *z)
  41 {
  42     int ln = s->nbits;
  43     int j, np, np2;
  44     int nblocks, nloops;
  45     register FFTComplex *p, *q;
  46    // FFTComplex *exptab = s->exptab;
  47     int l;
  48     fixed32 tmp_re, tmp_im;
  49     int tabshift = 10-ln;
  50
  51     np = 1 << ln;
  52
  53
  54     /* pass 0 */
  55
  56     p=&z[0];
  57     j=(np >> 1);
  58     do
  59     {
  60         BF(p[0].re, p[0].im, p[1].re, p[1].im,
  61            p[0].re, p[0].im, p[1].re, p[1].im);
  62         p+=2;
  63     }
  64     while (--j != 0);
  65
  66     /* pass 1 */
  67
  68
  69     p=&z[0];
  70     j=np >> 2;
  71     if (s->inverse)
  72     {
  73         do
  74         {
  75             BF(p[0].re, p[0].im, p[2].re, p[2].im,
  76                p[0].re, p[0].im, p[2].re, p[2].im);
  77             BF(p[1].re, p[1].im, p[3].re, p[3].im,
  78                p[1].re, p[1].im, -p[3].im, p[3].re);
  79             p+=4;
  80         }
  81         while (--j != 0);
  82     }
  83     else
  84     {
  85         do
  86         {
  87             BF(p[0].re, p[0].im, p[2].re, p[2].im,
  88                p[0].re, p[0].im, p[2].re, p[2].im);
  89             BF(p[1].re, p[1].im, p[3].re, p[3].im,
  90                p[1].re, p[1].im, p[3].im, -p[3].re);
  91             p+=4;
  92         }
  93         while (--j != 0);
  94     }
  95     /* pass 2 .. ln-1 */
  96
  97     nblocks = np >> 3;
  98     nloops = 1 << 2;
  99     np2 = np >> 1;
 100     do
 101     {
 102         p = z;
 103         q = z + nloops;
 104         for (j = 0; j < nblocks; ++j)
 105         {
 106             BF(p->re, p->im, q->re, q->im,
 107                p->re, p->im, q->re, q->im);
 108
 109             p++;
 110             q++;
 111             for(l = nblocks; l < np2; l += nblocks)
 112             {
 113                 CMUL(&tmp_re, &tmp_im, exptab0[(l<<tabshift)].re, exptab0[(l<<tabshift)].im, q->re, q->im);
 114                 //CMUL(&tmp_re, &tmp_im, exptab[l].re, exptab[l].im, q->re, q->im);
 115                 BF(p->re, p->im, q->re, q->im,
 116                    p->re, p->im, tmp_re, tmp_im);
 117                 p++;
 118                 q++;
 119             }
 120
 121             p += nloops;
 122             q += nloops;
 123         }
 124         nblocks = nblocks >> 1;
 125         nloops = nloops << 1;
 126     }
 127     while (nblocks != 0);
 128     return 0;
 129 }
 130
 131 int fft_init_global(void)
 132 {
 133     int i, n;
 134     fixed32 c1, s1, s2;
 135
 136     n=1<<10;
 137     s2 = 1 ? 1 : -1;
 138     for(i=0;i<(n/2);++i)
 139     {
 140         fixed32 ifix = itofix32(i);
 141         fixed32 nfix = itofix32(n);
 142         fixed32 res = fixdiv32(ifix,nfix);
 143
 144         s1 = fsincos(res<<16, &c1);
 145
 146         exptab0[i].re = c1;
 147         exptab0[i].im = s1*s2;
 148     }
 149
 150     return 0;
 151 }
 152