lib/lame/fft.c

   1 /*
   2 ** FFT and FHT routines
   3 **  Copyright 1988, 1993; Ron Mayer
   4 **
   5 **  fht(fz,n);
   6 **      Does a hartley transform of "n" points in the array "fz".
   7 **
   8 ** NOTE: This routine uses at least 2 patented algorithms, and may be
   9 **       under the restrictions of a bunch of different organizations.
  10 **       Although I wrote it completely myself; it is kind of a derivative
  11 **       of a routine I once authored and released under the GPL, so it
  12 **       may fall under the free software foundation's restrictions;
  13 **       it was worked on as a Stanford Univ project, so they claim
  14 **       some rights to it; it was further optimized at work here, so
  15 **       I think this company claims parts of it.  The patents are
  16 **       held by R. Bracewell (the FHT algorithm) and O. Buneman (the
  17 **       trig generator), both at Stanford Univ.
  18 **       If it were up to me, I'd say go do whatever you want with it;
  19 **       but it would be polite to give credit to the following people
  20 **       if you use this anywhere:
  21 **           Euler     - probable inventor of the fourier transform.
  22 **           Gauss     - probable inventor of the FFT.
  23 **           Hartley   - probable inventor of the hartley transform.
  24 **           Buneman   - for a really cool trig generator
  25 **           Mayer(me) - for authoring this particular version and
  26 **                       including all the optimizations in one package.
  27 **       Thanks,
  28 **       Ron Mayer; mayer@acuson.com
  29 ** and added some optimization by
  30 **           Mather    - idea of using lookup table
  31 **           Takehiro  - some dirty hack for speed up
  32 */
  33
  34 /* $Id: fft.c,v 1.1 2002/04/28 17:30:18 kramm Exp $ */
  35
  36 #include "config_static.h"
  37
  38 #include <math.h>
  39 #include "util.h"
  40 #include "fft.h"
  41
  42
  43
  44
  45 #ifdef WITH_DMALLOC
  46 #include <dmalloc.h>
  47 #endif
  48
  49 #define TRI_SIZE (5-1) /* 1024 =  4**5 */
  50
  51 static const FLOAT costab[TRI_SIZE*2] = {
  52   9.238795325112867e-01, 3.826834323650898e-01,
  53   9.951847266721969e-01, 9.801714032956060e-02,
  54   9.996988186962042e-01, 2.454122852291229e-02,
  55   9.999811752826011e-01, 6.135884649154475e-03
  56 };
  57
  58 static void fht(FLOAT *fz, int n)
  59 {
  60     const FLOAT *tri = costab;
  61     int           k4;
  62     FLOAT *fi, *fn, *gi;
  63
  64     n <<= 1;        /* to get BLKSIZE, because of 3DNow! ASM routine */
  65     fn = fz + n;
  66     k4 = 4;
  67     do {
  68         FLOAT s1, c1;
  69         int   i, k1, k2, k3, kx;
  70         kx  = k4 >> 1;
  71         k1  = k4;
  72         k2  = k4 << 1;
  73         k3  = k2 + k1;
  74         k4  = k2 << 1;
  75         fi  = fz;
  76         gi  = fi + kx;
  77         do {
  78             FLOAT f0,f1,f2,f3;
  79             f1      = fi[0]  - fi[k1];
  80             f0      = fi[0]  + fi[k1];
  81             f3      = fi[k2] - fi[k3];
  82             f2      = fi[k2] + fi[k3];
  83             fi[k2]  = f0     - f2;
  84             fi[0 ]  = f0     + f2;
  85             fi[k3]  = f1     - f3;
  86             fi[k1]  = f1     + f3;
  87             f1      = gi[0]  - gi[k1];
  88             f0      = gi[0]  + gi[k1];
  89             f3      = SQRT2  * gi[k3];
  90             f2      = SQRT2  * gi[k2];
  91             gi[k2]  = f0     - f2;
  92             gi[0 ]  = f0     + f2;
  93             gi[k3]  = f1     - f3;
  94             gi[k1]  = f1     + f3;
  95             gi     += k4;
  96             fi     += k4;
  97         } while (fi<fn);
  98         c1 = tri[0];
  99         s1 = tri[1];
 100         for (i = 1; i < kx; i++) {
 101             FLOAT c2,s2;
 102             c2 = 1 - (2*s1)*s1;
 103             s2 = (2*s1)*c1;
 104             fi = fz + i;
 105             gi = fz + k1 - i;
 106             do {
 107                 FLOAT a,b,g0,f0,f1,g1,f2,g2,f3,g3;
 108                 b       = s2*fi[k1] - c2*gi[k1];
 109                 a       = c2*fi[k1] + s2*gi[k1];
 110                 f1      = fi[0 ]    - a;
 111                 f0      = fi[0 ]    + a;
 112                 g1      = gi[0 ]    - b;
 113                 g0      = gi[0 ]    + b;
 114                 b       = s2*fi[k3] - c2*gi[k3];
 115                 a       = c2*fi[k3] + s2*gi[k3];
 116                 f3      = fi[k2]    - a;
 117                 f2      = fi[k2]    + a;
 118                 g3      = gi[k2]    - b;
 119                 g2      = gi[k2]    + b;
 120                 b       = s1*f2     - c1*g3;
 121                 a       = c1*f2     + s1*g3;
 122                 fi[k2]  = f0        - a;
 123                 fi[0 ]  = f0        + a;
 124                 gi[k3]  = g1        - b;
 125                 gi[k1]  = g1        + b;
 126                 b       = c1*g2     - s1*f3;
 127                 a       = s1*g2     + c1*f3;
 128                 gi[k2]  = g0        - a;
 129                 gi[0 ]  = g0        + a;
 130                 fi[k3]  = f1        - b;
 131                 fi[k1]  = f1        + b;
 132                 gi     += k4;
 133                 fi     += k4;
 134             } while (fi<fn);
 135             c2 = c1;
 136             c1 = c2 * tri[0] - s1 * tri[1];
 137             s1 = c2 * tri[1] + s1 * tri[0];
 138         }
 139         tri += 2;
 140     } while (k4<n);
 141 }
 142
 143 static const unsigned char rv_tbl[] = {
 144     0x00,    0x80,    0x40,    0xc0,    0x20,    0xa0,    0x60,    0xe0,
 145     0x10,    0x90,    0x50,    0xd0,    0x30,    0xb0,    0x70,    0xf0,
 146     0x08,    0x88,    0x48,    0xc8,    0x28,    0xa8,    0x68,    0xe8,
 147     0x18,    0x98,    0x58,    0xd8,    0x38,    0xb8,    0x78,    0xf8,
 148     0x04,    0x84,    0x44,    0xc4,    0x24,    0xa4,    0x64,    0xe4,
 149     0x14,    0x94,    0x54,    0xd4,    0x34,    0xb4,    0x74,    0xf4,
 150     0x0c,    0x8c,    0x4c,    0xcc,    0x2c,    0xac,    0x6c,    0xec,
 151     0x1c,    0x9c,    0x5c,    0xdc,    0x3c,    0xbc,    0x7c,    0xfc,
 152     0x02,    0x82,    0x42,    0xc2,    0x22,    0xa2,    0x62,    0xe2,
 153     0x12,    0x92,    0x52,    0xd2,    0x32,    0xb2,    0x72,    0xf2,
 154     0x0a,    0x8a,    0x4a,    0xca,    0x2a,    0xaa,    0x6a,    0xea,
 155     0x1a,    0x9a,    0x5a,    0xda,    0x3a,    0xba,    0x7a,    0xfa,
 156     0x06,    0x86,    0x46,    0xc6,    0x26,    0xa6,    0x66,    0xe6,
 157     0x16,    0x96,    0x56,    0xd6,    0x36,    0xb6,    0x76,    0xf6,
 158     0x0e,    0x8e,    0x4e,    0xce,    0x2e,    0xae,    0x6e,    0xee,
 159     0x1e,    0x9e,    0x5e,    0xde,    0x3e,    0xbe,    0x7e,    0xfe
 160 };
 161
 162 #define ch01(index)  (buffer[chn][index])
 163
 164 #define ml00(f) (window[i        ] * f(i))
 165 #define ml10(f) (window[i + 0x200] * f(i + 0x200))
 166 #define ml20(f) (window[i + 0x100] * f(i + 0x100))
 167 #define ml30(f) (window[i + 0x300] * f(i + 0x300))
 168
 169 #define ml01(f) (window[i + 0x001] * f(i + 0x001))
 170 #define ml11(f) (window[i + 0x201] * f(i + 0x201))
 171 #define ml21(f) (window[i + 0x101] * f(i + 0x101))
 172 #define ml31(f) (window[i + 0x301] * f(i + 0x301))
 173
 174 #define ms00(f) (window_s[i       ] * f(i + k))
 175 #define ms10(f) (window_s[0x7f - i] * f(i + k + 0x80))
 176 #define ms20(f) (window_s[i + 0x40] * f(i + k + 0x40))
 177 #define ms30(f) (window_s[0x3f - i] * f(i + k + 0xc0))
 178
 179 #define ms01(f) (window_s[i + 0x01] * f(i + k + 0x01))
 180 #define ms11(f) (window_s[0x7e - i] * f(i + k + 0x81))
 181 #define ms21(f) (window_s[i + 0x41] * f(i + k + 0x41))
 182 #define ms31(f) (window_s[0x3e - i] * f(i + k + 0xc1))
 183
 184
 185 void fft_short(lame_internal_flags * const gfc,
 186                 FLOAT x_real[3][BLKSIZE_s], int chn, const sample_t *buffer[2])
 187 {
 188     const FLOAT*  window_s = (const FLOAT *)&gfc->window_s[0];
 189     int           i;
 190     int           j;
 191     int           b;
 192
 193     for (b = 0; b < 3; b++) {
 194         FLOAT *x = &x_real[b][BLKSIZE_s / 2];
 195         short k = (576 / 3) * (b + 1);
 196         j = BLKSIZE_s / 8 - 1;
 197         do {
 198           FLOAT f0,f1,f2,f3, w;
 199
 200           i = rv_tbl[j << 2];
 201
 202           f0 = ms00(ch01); w = ms10(ch01); f1 = f0 - w; f0 = f0 + w;
 203           f2 = ms20(ch01); w = ms30(ch01); f3 = f2 - w; f2 = f2 + w;
 204
 205           x -= 4;
 206           x[0] = f0 + f2;
 207           x[2] = f0 - f2;
 208           x[1] = f1 + f3;
 209           x[3] = f1 - f3;
 210
 211           f0 = ms01(ch01); w = ms11(ch01); f1 = f0 - w; f0 = f0 + w;
 212           f2 = ms21(ch01); w = ms31(ch01); f3 = f2 - w; f2 = f2 + w;
 213
 214           x[BLKSIZE_s / 2 + 0] = f0 + f2;
 215           x[BLKSIZE_s / 2 + 2] = f0 - f2;
 216           x[BLKSIZE_s / 2 + 1] = f1 + f3;
 217           x[BLKSIZE_s / 2 + 3] = f1 - f3;
 218         } while (--j >= 0);
 219
 220         gfc->fft_fht(x, BLKSIZE_s/2);
 221         /* BLKSIZE_s/2 because of 3DNow! ASM routine */
 222     }
 223 }
 224
 225 void fft_long(lame_internal_flags * const gfc,
 226                FLOAT x[BLKSIZE], int chn, const sample_t *buffer[2] )
 227 {
 228     const FLOAT*  window = (const FLOAT *)&gfc->window[0];
 229     int           i;
 230     int           jj = BLKSIZE / 8 - 1;
 231     x += BLKSIZE / 2;
 232
 233     do {
 234       FLOAT f0,f1,f2,f3, w;
 235
 236       i = rv_tbl[jj];
 237       f0 = ml00(ch01); w = ml10(ch01); f1 = f0 - w; f0 = f0 + w;
 238       f2 = ml20(ch01); w = ml30(ch01); f3 = f2 - w; f2 = f2 + w;
 239
 240       x -= 4;
 241       x[0] = f0 + f2;
 242       x[2] = f0 - f2;
 243       x[1] = f1 + f3;
 244       x[3] = f1 - f3;
 245
 246       f0 = ml01(ch01); w = ml11(ch01); f1 = f0 - w; f0 = f0 + w;
 247       f2 = ml21(ch01); w = ml31(ch01); f3 = f2 - w; f2 = f2 + w;
 248
 249       x[BLKSIZE / 2 + 0] = f0 + f2;
 250       x[BLKSIZE / 2 + 2] = f0 - f2;
 251       x[BLKSIZE / 2 + 1] = f1 + f3;
 252       x[BLKSIZE / 2 + 3] = f1 - f3;
 253     } while (--jj >= 0);
 254
 255     gfc->fft_fht(x, BLKSIZE/2);
 256     /* BLKSIZE/2 because of 3DNow! ASM routine */
 257 }
 258
 259
 260 void init_fft(lame_internal_flags * const gfc)
 261 {
 262     FLOAT *window   = &gfc->window[0];
 263     FLOAT *window_s = &gfc->window_s[0];
 264     int i;
 265
 266 #if 0
 267     if (gfc->nsPsy.use) {
 268       for (i = 0; i < BLKSIZE ; i++)
 269         /* blackman window */
 270         window[i] = 0.42-0.5*cos(2*PI*i/(BLKSIZE-1))+0.08*cos(4*PI*i/(BLKSIZE-1));
 271     } else {
 272       /*
 273        * calculate HANN window coefficients
 274        */
 275       for (i = 0; i < BLKSIZE ; i++)
 276         window[i] = 0.5 * (1.0 - cos(2.0 * PI * (i + 0.5) / BLKSIZE));
 277     }
 278 #endif
 279
 280     // The type of window used here will make no real difference, but
 281     // in the interest of merging nspsytune stuff - switch to blackman window
 282     for (i = 0; i < BLKSIZE ; i++)
 283       /* blackman window */
 284       window[i] = 0.42-0.5*cos(2*PI*(i+.5)/BLKSIZE)+
 285         0.08*cos(4*PI*(i+.5)/BLKSIZE);
 286
 287     for (i = 0; i < BLKSIZE_s/2 ; i++)
 288         window_s[i] = 0.5 * (1.0 - cos(2.0 * PI * (i + 0.5) / BLKSIZE_s));
 289
 290 #ifdef HAVE_NASM
 291     if (gfc->CPU_features.AMD_3DNow) {
 292         extern void fht_3DN(FLOAT *fz, int n);
 293         gfc->fft_fht = fht_3DN;
 294     } else
 295 #endif
 296 #ifdef USE_FFTSSE
 297     if (gfc->CPU_features.SIMD) {
 298         extern void fht_SSE(FLOAT *fz, int n);
 299         gfc->fft_fht = fht_SSE;
 300     } else
 301 #endif
 302 #ifdef USE_FFTFPU
 303     if (gfc->CPU_features.i387) {
 304         extern void fht_FPU(FLOAT *fz, int n);
 305         gfc->fft_fht = fht_FPU;
 306     } else
 307 #endif
 308         gfc->fft_fht = fht;
 309 }