fft.c
上传用户:kjfoods
上传日期:2020-07-06
资源大小:29949k
文件大小:3k
源码类别:

midi

开发平台:

Unix_Linux

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