diff options
Diffstat (limited to 'libmp3lame/fft.c')
| -rw-r--r-- | libmp3lame/fft.c | 322 |
1 files changed, 322 insertions, 0 deletions
diff --git a/libmp3lame/fft.c b/libmp3lame/fft.c new file mode 100644 index 0000000..ec036bf --- /dev/null +++ b/libmp3lame/fft.c @@ -0,0 +1,322 @@ +/* +** FFT and FHT routines +** Copyright 1988, 1993; Ron Mayer +** Copyright (c) 1999-2000 Takehiro Tominaga +** +** fht(fz,n); +** Does a hartley transform of "n" points in the array "fz". +** +** NOTE: This routine uses at least 2 patented algorithms, and may be +** under the restrictions of a bunch of different organizations. +** Although I wrote it completely myself; it is kind of a derivative +** of a routine I once authored and released under the GPL, so it +** may fall under the free software foundation's restrictions; +** it was worked on as a Stanford Univ project, so they claim +** some rights to it; it was further optimized at work here, so +** I think this company claims parts of it. The patents are +** held by R. Bracewell (the FHT algorithm) and O. Buneman (the +** trig generator), both at Stanford Univ. +** If it were up to me, I'd say go do whatever you want with it; +** but it would be polite to give credit to the following people +** if you use this anywhere: +** Euler - probable inventor of the fourier transform. +** Gauss - probable inventor of the FFT. +** Hartley - probable inventor of the hartley transform. +** Buneman - for a really cool trig generator +** Mayer(me) - for authoring this particular version and +** including all the optimizations in one package. +** Thanks, +** Ron Mayer; mayer@acuson.com +** and added some optimization by +** Mather - idea of using lookup table +** Takehiro - some dirty hack for speed up +*/ + +/* $Id: fft.c,v 1.33 2008/04/12 18:18:06 robert Exp $ */ + +#ifdef HAVE_CONFIG_H +# include <config.h> +#endif + +#include "lame.h" +#include "machine.h" +#include "encoder.h" +#include "util.h" +#include "fft.h" + + + + + +#define TRI_SIZE (5-1) /* 1024 = 4**5 */ + +/* fft.c */ +static FLOAT window[BLKSIZE], window_s[BLKSIZE_s / 2]; + +static const FLOAT costab[TRI_SIZE * 2] = { + 9.238795325112867e-01, 3.826834323650898e-01, + 9.951847266721969e-01, 9.801714032956060e-02, + 9.996988186962042e-01, 2.454122852291229e-02, + 9.999811752826011e-01, 6.135884649154475e-03 +}; + +static void +fht(FLOAT * fz, int n) +{ + const FLOAT *tri = costab; + int k4; + FLOAT *fi, *gi; + FLOAT const *fn; + + n <<= 1; /* to get BLKSIZE, because of 3DNow! ASM routine */ + fn = fz + n; + k4 = 4; + do { + FLOAT s1, c1; + int i, k1, k2, k3, kx; + kx = k4 >> 1; + k1 = k4; + k2 = k4 << 1; + k3 = k2 + k1; + k4 = k2 << 1; + fi = fz; + gi = fi + kx; + do { + FLOAT f0, f1, f2, f3; + f1 = fi[0] - fi[k1]; + f0 = fi[0] + fi[k1]; + f3 = fi[k2] - fi[k3]; + f2 = fi[k2] + fi[k3]; + fi[k2] = f0 - f2; + fi[0] = f0 + f2; + fi[k3] = f1 - f3; + fi[k1] = f1 + f3; + f1 = gi[0] - gi[k1]; + f0 = gi[0] + gi[k1]; + f3 = SQRT2 * gi[k3]; + f2 = SQRT2 * gi[k2]; + gi[k2] = f0 - f2; + gi[0] = f0 + f2; + gi[k3] = f1 - f3; + gi[k1] = f1 + f3; + gi += k4; + fi += k4; + } while (fi < fn); + c1 = tri[0]; + s1 = tri[1]; + for (i = 1; i < kx; i++) { + FLOAT c2, s2; + c2 = 1 - (2 * s1) * s1; + s2 = (2 * s1) * c1; + fi = fz + i; + gi = fz + k1 - i; + do { + FLOAT a, b, g0, f0, f1, g1, f2, g2, f3, g3; + b = s2 * fi[k1] - c2 * gi[k1]; + a = c2 * fi[k1] + s2 * gi[k1]; + f1 = fi[0] - a; + f0 = fi[0] + a; + g1 = gi[0] - b; + g0 = gi[0] + b; + b = s2 * fi[k3] - c2 * gi[k3]; + a = c2 * fi[k3] + s2 * gi[k3]; + f3 = fi[k2] - a; + f2 = fi[k2] + a; + g3 = gi[k2] - b; + g2 = gi[k2] + b; + b = s1 * f2 - c1 * g3; + a = c1 * f2 + s1 * g3; + fi[k2] = f0 - a; + fi[0] = f0 + a; + gi[k3] = g1 - b; + gi[k1] = g1 + b; + b = c1 * g2 - s1 * f3; + a = s1 * g2 + c1 * f3; + gi[k2] = g0 - a; + gi[0] = g0 + a; + fi[k3] = f1 - b; + fi[k1] = f1 + b; + gi += k4; + fi += k4; + } while (fi < fn); + c2 = c1; + c1 = c2 * tri[0] - s1 * tri[1]; + s1 = c2 * tri[1] + s1 * tri[0]; + } + tri += 2; + } while (k4 < n); +} + +static const unsigned char rv_tbl[] = { + 0x00, 0x80, 0x40, 0xc0, 0x20, 0xa0, 0x60, 0xe0, + 0x10, 0x90, 0x50, 0xd0, 0x30, 0xb0, 0x70, 0xf0, + 0x08, 0x88, 0x48, 0xc8, 0x28, 0xa8, 0x68, 0xe8, + 0x18, 0x98, 0x58, 0xd8, 0x38, 0xb8, 0x78, 0xf8, + 0x04, 0x84, 0x44, 0xc4, 0x24, 0xa4, 0x64, 0xe4, + 0x14, 0x94, 0x54, 0xd4, 0x34, 0xb4, 0x74, 0xf4, + 0x0c, 0x8c, 0x4c, 0xcc, 0x2c, 0xac, 0x6c, 0xec, + 0x1c, 0x9c, 0x5c, 0xdc, 0x3c, 0xbc, 0x7c, 0xfc, + 0x02, 0x82, 0x42, 0xc2, 0x22, 0xa2, 0x62, 0xe2, + 0x12, 0x92, 0x52, 0xd2, 0x32, 0xb2, 0x72, 0xf2, + 0x0a, 0x8a, 0x4a, 0xca, 0x2a, 0xaa, 0x6a, 0xea, + 0x1a, 0x9a, 0x5a, 0xda, 0x3a, 0xba, 0x7a, 0xfa, + 0x06, 0x86, 0x46, 0xc6, 0x26, 0xa6, 0x66, 0xe6, + 0x16, 0x96, 0x56, 0xd6, 0x36, 0xb6, 0x76, 0xf6, + 0x0e, 0x8e, 0x4e, 0xce, 0x2e, 0xae, 0x6e, 0xee, + 0x1e, 0x9e, 0x5e, 0xde, 0x3e, 0xbe, 0x7e, 0xfe +}; + +#define ch01(index) (buffer[chn][index]) + +#define ml00(f) (window[i ] * f(i)) +#define ml10(f) (window[i + 0x200] * f(i + 0x200)) +#define ml20(f) (window[i + 0x100] * f(i + 0x100)) +#define ml30(f) (window[i + 0x300] * f(i + 0x300)) + +#define ml01(f) (window[i + 0x001] * f(i + 0x001)) +#define ml11(f) (window[i + 0x201] * f(i + 0x201)) +#define ml21(f) (window[i + 0x101] * f(i + 0x101)) +#define ml31(f) (window[i + 0x301] * f(i + 0x301)) + +#define ms00(f) (window_s[i ] * f(i + k)) +#define ms10(f) (window_s[0x7f - i] * f(i + k + 0x80)) +#define ms20(f) (window_s[i + 0x40] * f(i + k + 0x40)) +#define ms30(f) (window_s[0x3f - i] * f(i + k + 0xc0)) + +#define ms01(f) (window_s[i + 0x01] * f(i + k + 0x01)) +#define ms11(f) (window_s[0x7e - i] * f(i + k + 0x81)) +#define ms21(f) (window_s[i + 0x41] * f(i + k + 0x41)) +#define ms31(f) (window_s[0x3e - i] * f(i + k + 0xc1)) + + +void +fft_short(lame_internal_flags const *const gfc, + FLOAT x_real[3][BLKSIZE_s], int chn, const sample_t * buffer[2]) +{ + int i; + int j; + int b; + + for (b = 0; b < 3; b++) { + FLOAT *x = &x_real[b][BLKSIZE_s / 2]; + short const k = (576 / 3) * (b + 1); + j = BLKSIZE_s / 8 - 1; + do { + FLOAT f0, f1, f2, f3, w; + + i = rv_tbl[j << 2]; + + f0 = ms00(ch01); + w = ms10(ch01); + f1 = f0 - w; + f0 = f0 + w; + f2 = ms20(ch01); + w = ms30(ch01); + f3 = f2 - w; + f2 = f2 + w; + + x -= 4; + x[0] = f0 + f2; + x[2] = f0 - f2; + x[1] = f1 + f3; + x[3] = f1 - f3; + + f0 = ms01(ch01); + w = ms11(ch01); + f1 = f0 - w; + f0 = f0 + w; + f2 = ms21(ch01); + w = ms31(ch01); + f3 = f2 - w; + f2 = f2 + w; + + x[BLKSIZE_s / 2 + 0] = f0 + f2; + x[BLKSIZE_s / 2 + 2] = f0 - f2; + x[BLKSIZE_s / 2 + 1] = f1 + f3; + x[BLKSIZE_s / 2 + 3] = f1 - f3; + } while (--j >= 0); + + gfc->fft_fht(x, BLKSIZE_s / 2); + /* BLKSIZE_s/2 because of 3DNow! ASM routine */ + } +} + +void +fft_long(lame_internal_flags const *const gfc, + FLOAT x[BLKSIZE], int chn, const sample_t * buffer[2]) +{ + int i; + int jj = BLKSIZE / 8 - 1; + x += BLKSIZE / 2; + + do { + FLOAT f0, f1, f2, f3, w; + + i = rv_tbl[jj]; + f0 = ml00(ch01); + w = ml10(ch01); + f1 = f0 - w; + f0 = f0 + w; + f2 = ml20(ch01); + w = ml30(ch01); + f3 = f2 - w; + f2 = f2 + w; + + x -= 4; + x[0] = f0 + f2; + x[2] = f0 - f2; + x[1] = f1 + f3; + x[3] = f1 - f3; + + f0 = ml01(ch01); + w = ml11(ch01); + f1 = f0 - w; + f0 = f0 + w; + f2 = ml21(ch01); + w = ml31(ch01); + f3 = f2 - w; + f2 = f2 + w; + + x[BLKSIZE / 2 + 0] = f0 + f2; + x[BLKSIZE / 2 + 2] = f0 - f2; + x[BLKSIZE / 2 + 1] = f1 + f3; + x[BLKSIZE / 2 + 3] = f1 - f3; + } while (--jj >= 0); + + gfc->fft_fht(x, BLKSIZE / 2); + /* BLKSIZE/2 because of 3DNow! ASM routine */ +} + +#ifdef HAVE_NASM +extern void fht_3DN(FLOAT * fz, int n); +extern void fht_SSE(FLOAT * fz, int n); +#endif + +void +init_fft(lame_internal_flags * const gfc) +{ + int i; + + /* The type of window used here will make no real difference, but */ + /* in the interest of merging nspsytune stuff - switch to blackman window */ + for (i = 0; i < BLKSIZE; i++) + /* blackman window */ + window[i] = 0.42 - 0.5 * cos(2 * PI * (i + .5) / BLKSIZE) + + 0.08 * cos(4 * PI * (i + .5) / BLKSIZE); + + for (i = 0; i < BLKSIZE_s / 2; i++) + window_s[i] = 0.5 * (1.0 - cos(2.0 * PI * (i + 0.5) / BLKSIZE_s)); + +#ifdef HAVE_NASM + if (gfc->CPU_features.AMD_3DNow) { + gfc->fft_fht = fht_3DN; + } + else if (gfc->CPU_features.SSE) { + gfc->fft_fht = fht_SSE; + } + else +#endif + { + gfc->fft_fht = fht; + } +} |