www.pudn.com > lame3(mp3编码程序和资料).zip > fft.c
/*
** FFT and FHT routines
** Copyright 1988, 1993; Ron Mayer
**
** 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
*/
#include
#include "util.h"
#include "psymodel.h"
#include "lame.h"
#define TRI_SIZE (5-1) /* 1024 = 4**5 */
static FLOAT costab[TRI_SIZE*2];
static FLOAT window[BLKSIZE / 2], window_s[BLKSIZE_s / 2];
static INLINE void fht(FLOAT *fz, int n)
{
short k4;
FLOAT *fi, *fn, *gi;
FLOAT *tri;
fn = fz + n;
tri = &costab[0];
k4 = 4;
do {
FLOAT s1, c1;
short 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= 0);
} else if (chn == 2) {
do {
FLOAT f0,f1,f2,f3, w;
i = rv_tbl[j << 2];
f0 = ms00(ch2); w = ms10(ch2); f1 = f0 - w; f0 = f0 + w;
f2 = ms20(ch2); w = ms30(ch2); 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(ch2); w = ms11(ch2); f1 = f0 - w; f0 = f0 + w;
f2 = ms21(ch2); w = ms31(ch2); 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);
} else {
do {
FLOAT f0,f1,f2,f3, w;
i = rv_tbl[j << 2];
f0 = ms00(ch3); w = ms10(ch3); f1 = f0 - w; f0 = f0 + w;
f2 = ms20(ch3); w = ms30(ch3); 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(ch3); w = ms11(ch3); f1 = f0 - w; f0 = f0 + w;
f2 = ms21(ch3); w = ms31(ch3); 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);
}
fht(x, BLKSIZE_s);
}
}
void fft_long(
FLOAT x[BLKSIZE], int chn, short *buffer[2])
{
short i,jj = BLKSIZE / 8 - 1;
x += BLKSIZE / 2;
if (chn < 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);
} else if (chn == 2) {
do {
FLOAT f0,f1,f2,f3, w;
i = rv_tbl[jj];
f0 = ml00(ch2); w = ml10(ch2); f1 = f0 - w; f0 = f0 + w;
f2 = ml20(ch2); w = ml30(ch2); 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(ch2); w = ml11(ch2); f1 = f0 - w; f0 = f0 + w;
f2 = ml21(ch2); w = ml31(ch2); 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);
} else {
do {
FLOAT f0,f1,f2,f3, w;
i = rv_tbl[jj];
f0 = ml00(ch3); w = ml10(ch3); f1 = f0 - w; f0 = f0 + w;
f2 = ml20(ch3); w = ml30(ch3); 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(ch3); w = ml11(ch3); f1 = f0 - w; f0 = f0 + w;
f2 = ml21(ch3); w = ml31(ch3); 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);
}
fht(x, BLKSIZE);
}
void init_fft(void)
{
int i;
FLOAT r = PI*0.125;
for (i = 0; i < TRI_SIZE; i++) {
costab[i*2 ] = cos(r);
costab[i*2+1] = sin(r);
r *= 0.25;
}
/*
* calculate HANN window coefficients
*/
for (i = 0; i < BLKSIZE / 2; i++)
window[i] = 0.5 * (1.0 - cos(2.0 * PI * (i + 0.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));
}