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