git.sesse.net Git - vlc/blob - modules/codec/wmafixed/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 <inttypes.h>
  21 #include "fft.h"
  22 #include "wmafixed.h"
  23
  24 #define IBSS_ATTR
  25 #define ICONST_ATTR
  26 #define ICODE_ATTR
  27
  28 FFTComplex  exptab0[512] IBSS_ATTR;
  29
  30 /* butter fly op */
  31 #define BF(pre, pim, qre, qim, pre1, pim1, qre1, qim1) \
  32 { \
  33   int32_t ax, ay, bx, by; \
  34   bx=pre1; \
  35   by=pim1; \
  36   ax=qre1; \
  37   ay=qim1; \
  38   pre = (bx + ax); \
  39   pim = (by + ay); \
  40   qre = (bx - ax); \
  41   qim = (by - ay); \
  42 }
  43
  44
  45 int fft_calc_unscaled(FFTContext *s, FFTComplex *z)
  46 {
  47     int ln = s->nbits;
  48     int j, np, np2;
  49     int nblocks, nloops;
  50     register FFTComplex *p, *q;
  51     int l;
  52     int32_t tmp_re, tmp_im;
  53     int tabshift = 10-ln;
  54
  55     np = 1 << ln;
  56
  57     /* pass 0 */
  58     p=&z[0];
  59     j=(np >> 1);
  60     do
  61     {
  62         BF(p[0].re, p[0].im, p[1].re, p[1].im,
  63            p[0].re, p[0].im, p[1].re, p[1].im);
  64         p+=2;
  65     }
  66     while (--j != 0);
  67
  68     /* pass 1 */
  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
  96     /* pass 2 .. ln-1 */
  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
 129     return 0;
 130 }
 131
 132 int fft_init_global(void)
 133 {
 134     int i, n;
 135     int32_t c1, s1, s2;
 136
 137     n=1<<10;
 138     s2 = 1 ? 1 : -1;
 139
 140     for(i=0;i<(n/2);++i)
 141     {
 142         int32_t ifix = itofix32(i);
 143         int32_t nfix = itofix32(n);
 144         int32_t res = fixdiv32(ifix,nfix);
 145
 146         s1 = fsincos(res<<16, &c1);
 147
 148         exptab0[i].re = c1;
 149         exptab0[i].im = s1*s2;
 150     }
 151
 152     return 0;
 153 }