libavcodec/fft.c

   1 /*
   2  * FFT/IFFT transforms
   3  * Copyright (c) 2002 Fabrice Bellard.
   4  *
   5  * This file is part of FFmpeg.
   6  *
   7  * FFmpeg is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU Lesser General Public
   9  * License as published by the Free Software Foundation; either
  10  * version 2.1 of the License, or (at your option) any later version.
  11  *
  12  * FFmpeg is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15  * Lesser General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU Lesser General Public
  18  * License along with FFmpeg; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  20  */
  21
  22 /**
  23  * @file fft.c
  24  * FFT/IFFT transforms.
  25  */
  26
  27 #include "dsputil.h"
  28
  29 /**
  30  * The size of the FFT is 2^nbits. If inverse is TRUE, inverse FFT is
  31  * done
  32  */
  33 int ff_fft_init(FFTContext *s, int nbits, int inverse)
  34 {
  35     int i, j, m, n;
  36     float alpha, c1, s1, s2;
  37
  38     s->nbits = nbits;
  39     n = 1 << nbits;
  40
  41     s->exptab = av_malloc((n / 2) * sizeof(FFTComplex));
  42     if (!s->exptab)
  43         goto fail;
  44     s->revtab = av_malloc(n * sizeof(uint16_t));
  45     if (!s->revtab)
  46         goto fail;
  47     s->inverse = inverse;
  48
  49     s2 = inverse ? 1.0 : -1.0;
  50
  51     for(i=0;i<(n/2);i++) {
  52         alpha = 2 * M_PI * (float)i / (float)n;
  53         c1 = cos(alpha);
  54         s1 = sin(alpha) * s2;
  55         s->exptab[i].re = c1;
  56         s->exptab[i].im = s1;
  57     }
  58     s->fft_calc = ff_fft_calc_c;
  59     s->imdct_calc = ff_imdct_calc;
  60     s->exptab1 = NULL;
  61
  62     /* compute constant table for HAVE_SSE version */
  63 #if defined(HAVE_MMX) \
  64     || (defined(HAVE_ALTIVEC) && !defined(ALTIVEC_USE_REFERENCE_C_CODE))
  65     {
  66         int has_vectors = mm_support();
  67
  68         if (has_vectors) {
  69 #if defined(HAVE_MMX)
  70             if (has_vectors & MM_3DNOWEXT) {
  71                 /* 3DNowEx for K7/K8 */
  72                 s->imdct_calc = ff_imdct_calc_3dn2;
  73                 s->fft_calc = ff_fft_calc_3dn2;
  74             } else if (has_vectors & MM_3DNOW) {
  75                 /* 3DNow! for K6-2/3 */
  76                 s->fft_calc = ff_fft_calc_3dn;
  77             } else if (has_vectors & MM_SSE) {
  78                 /* SSE for P3/P4 */
  79                 s->imdct_calc = ff_imdct_calc_sse;
  80                 s->fft_calc = ff_fft_calc_sse;
  81             }
  82 #else /* HAVE_MMX */
  83             if (has_vectors & MM_ALTIVEC)
  84                 s->fft_calc = ff_fft_calc_altivec;
  85 #endif
  86         }
  87         if (s->fft_calc != ff_fft_calc_c) {
  88             int np, nblocks, np2, l;
  89             FFTComplex *q;
  90
  91             np = 1 << nbits;
  92             nblocks = np >> 3;
  93             np2 = np >> 1;
  94             s->exptab1 = av_malloc(np * 2 * sizeof(FFTComplex));
  95             if (!s->exptab1)
  96                 goto fail;
  97             q = s->exptab1;
  98             do {
  99                 for(l = 0; l < np2; l += 2 * nblocks) {
 100                     *q++ = s->exptab[l];
 101                     *q++ = s->exptab[l + nblocks];
 102
 103                     q->re = -s->exptab[l].im;
 104                     q->im = s->exptab[l].re;
 105                     q++;
 106                     q->re = -s->exptab[l + nblocks].im;
 107                     q->im = s->exptab[l + nblocks].re;
 108                     q++;
 109                 }
 110                 nblocks = nblocks >> 1;
 111             } while (nblocks != 0);
 112             av_freep(&s->exptab);
 113         }
 114     }
 115 #endif
 116
 117     /* compute bit reverse table */
 118
 119     for(i=0;i<n;i++) {
 120         m=0;
 121         for(j=0;j<nbits;j++) {
 122             m |= ((i >> j) & 1) << (nbits-j-1);
 123         }
 124         s->revtab[i]=m;
 125     }
 126     return 0;
 127  fail:
 128     av_freep(&s->revtab);
 129     av_freep(&s->exptab);
 130     av_freep(&s->exptab1);
 131     return -1;
 132 }
 133
 134 /* butter fly op */
 135 #define BF(pre, pim, qre, qim, pre1, pim1, qre1, qim1) \
 136 {\
 137   FFTSample ax, ay, bx, by;\
 138   bx=pre1;\
 139   by=pim1;\
 140   ax=qre1;\
 141   ay=qim1;\
 142   pre = (bx + ax);\
 143   pim = (by + ay);\
 144   qre = (bx - ax);\
 145   qim = (by - ay);\
 146 }
 147
 148 #define MUL16(a,b) ((a) * (b))
 149
 150 #define CMUL(pre, pim, are, aim, bre, bim) \
 151 {\
 152    pre = (MUL16(are, bre) - MUL16(aim, bim));\
 153    pim = (MUL16(are, bim) + MUL16(bre, aim));\
 154 }
 155
 156 /**
 157  * Do a complex FFT with the parameters defined in ff_fft_init(). The
 158  * input data must be permuted before with s->revtab table. No
 159  * 1.0/sqrt(n) normalization is done.
 160  */
 161 void ff_fft_calc_c(FFTContext *s, FFTComplex *z)
 162 {
 163     int ln = s->nbits;
 164     int j, np, np2;
 165     int nblocks, nloops;
 166     register FFTComplex *p, *q;
 167     FFTComplex *exptab = s->exptab;
 168     int l;
 169     FFTSample tmp_re, tmp_im;
 170
 171     np = 1 << ln;
 172
 173     /* pass 0 */
 174
 175     p=&z[0];
 176     j=(np >> 1);
 177     do {
 178         BF(p[0].re, p[0].im, p[1].re, p[1].im,
 179            p[0].re, p[0].im, p[1].re, p[1].im);
 180         p+=2;
 181     } while (--j != 0);
 182
 183     /* pass 1 */
 184
 185
 186     p=&z[0];
 187     j=np >> 2;
 188     if (s->inverse) {
 189         do {
 190             BF(p[0].re, p[0].im, p[2].re, p[2].im,
 191                p[0].re, p[0].im, p[2].re, p[2].im);
 192             BF(p[1].re, p[1].im, p[3].re, p[3].im,
 193                p[1].re, p[1].im, -p[3].im, p[3].re);
 194             p+=4;
 195         } while (--j != 0);
 196     } else {
 197         do {
 198             BF(p[0].re, p[0].im, p[2].re, p[2].im,
 199                p[0].re, p[0].im, p[2].re, p[2].im);
 200             BF(p[1].re, p[1].im, p[3].re, p[3].im,
 201                p[1].re, p[1].im, p[3].im, -p[3].re);
 202             p+=4;
 203         } while (--j != 0);
 204     }
 205     /* pass 2 .. ln-1 */
 206
 207     nblocks = np >> 3;
 208     nloops = 1 << 2;
 209     np2 = np >> 1;
 210     do {
 211         p = z;
 212         q = z + nloops;
 213         for (j = 0; j < nblocks; ++j) {
 214             BF(p->re, p->im, q->re, q->im,
 215                p->re, p->im, q->re, q->im);
 216
 217             p++;
 218             q++;
 219             for(l = nblocks; l < np2; l += nblocks) {
 220                 CMUL(tmp_re, tmp_im, exptab[l].re, exptab[l].im, q->re, q->im);
 221                 BF(p->re, p->im, q->re, q->im,
 222                    p->re, p->im, tmp_re, tmp_im);
 223                 p++;
 224                 q++;
 225             }
 226
 227             p += nloops;
 228             q += nloops;
 229         }
 230         nblocks = nblocks >> 1;
 231         nloops = nloops << 1;
 232     } while (nblocks != 0);
 233 }
 234
 235 /**
 236  * Do the permutation needed BEFORE calling ff_fft_calc()
 237  */
 238 void ff_fft_permute(FFTContext *s, FFTComplex *z)
 239 {
 240     int j, k, np;
 241     FFTComplex tmp;
 242     const uint16_t *revtab = s->revtab;
 243
 244     /* reverse */
 245     np = 1 << s->nbits;
 246     for(j=0;j<np;j++) {
 247         k = revtab[j];
 248         if (k < j) {
 249             tmp = z[k];
 250             z[k] = z[j];
 251             z[j] = tmp;
 252         }
 253     }
 254 }
 255
 256 void ff_fft_end(FFTContext *s)
 257 {
 258     av_freep(&s->revtab);
 259     av_freep(&s->exptab);
 260     av_freep(&s->exptab1);
 261 }
 262