ref: 3296912cd15b69be806275733ef75ab74df5194c
parent: 4367273bd26e526e94932b963702d929f6780537
author: lenox <lenox>
date: Mon Dec 20 05:24:28 EST 1999
some cleanup
--- a/fastfft.c
+++ b/fastfft.c
@@ -1,6 +1,4 @@
#include "transfo.h"
-#ifndef USE_ORIG_FFT
-
#include "fastfft.h"
int unscambled64[64];
@@ -1644,5 +1642,3 @@
unscambled512[i]=PFFTW(permutation_512)(i);
}
}
-
-#endif
\ No newline at end of file
--- a/fastfft.h
+++ b/fastfft.h
@@ -5,11 +5,7 @@
#define CONCAT(a, b) CONCAT_AUX(a,b)
#define PFFTW(name) CONCAT(pfftw_d_, name)
-#define PFFTWI(name) CONCAT(pfftwi_d_, name)
#define FFTW_KONST(x) ((fftw_real) x)
-void PFFTW(twiddle_2)(fftw_complex *A, const fftw_complex *W, int iostride);
void PFFTW(twiddle_4)(fftw_complex *A, const fftw_complex *W, int iostride);
-void PFFTWI(twiddle_2)(fftw_complex *A, const fftw_complex *W, int iostride);
-void PFFTWI(twiddle_4)(fftw_complex *A, const fftw_complex *W, int iostride);
--- a/transfo.c
+++ b/transfo.c
@@ -1,14 +1,9 @@
#include <math.h>
#include "transfo.h"
-#ifdef USE_ORIG_FFT
- double FFTarray [1024]; /* the array for in-place FFT */
-#else
- fftw_complex_d FFTarray[512]; /* the array for in-place FFT */
- extern int unscambled64[64]; /* the permutation array for FFT64*/
- extern int unscambled512[512]; /* the permutation array for FFT512*/
- int unscambled;
-#endif
+fftw_complex_d FFTarray[512]; /* the array for in-place FFT */
+extern int unscambled64[64]; /* the permutation array for FFT64*/
+extern int unscambled512[512]; /* the permutation array for FFT512*/
#ifndef M_PI
#define M_PI 3.14159265358979323846
@@ -24,9 +19,7 @@
void MDCT (double *data, int N) {-#ifndef USE_ORIG_FFT
static int init = 1;
-#endif
double tempr, tempi, c, s, cold, cfreq, sfreq; /* temps for pre and post twiddle */
double freq = 2.0 * M_PI / N;
double fac,cosfreq8,sinfreq8;
@@ -39,13 +32,12 @@
int a2 = a >> 1;
int a4 = a >> 2;
int b4 = b >> 2;
+ int unscambled;
-#ifndef USE_ORIG_FFT
if (init) {init = 0;
MakeFFTOrder();
}
-#endif
/* Choosing to allocate 2/N factor to Inverse Xform! */
fac = 2.; /* 2 from MDCT inverse to forward */
@@ -75,13 +67,8 @@
}
/* calculate pre-twiddled FFT input */
-#ifdef USE_ORIG_FFT
- FFTarray [2 * i] = tempr * c + tempi * s;
- FFTarray [2 * i + 1] = tempi * c - tempr * s;
-#else
FFTarray [i].re = tempr * c + tempi * s;
FFTarray [i].im = tempi * c - tempr * s;
-#endif
/* use recurrence to prepare cosine and sine for next value of i */
cold = c;
@@ -90,15 +77,11 @@
}
/* Perform in-place complex FFT of length N/4 */
-#ifdef USE_ORIG_FFT
- FFT (FFTarray, N4, -1);
-#else
switch (N) {case 256: pfftw_d_64(FFTarray);
break;
case 2048:pfftw_d_512(FFTarray);
}
-#endif
/* prepare for recurrence relations in post-twiddle */
c = cosfreq8;
@@ -109,10 +92,6 @@
/* get post-twiddled FFT output */
/* Note: fac allocates 4/N factor from IFFT to forward and inverse */
-#ifdef USE_ORIG_FFT
- tempr = fac * (FFTarray [2 * i] * c + FFTarray [2 * i + 1] * s);
- tempi = fac * (FFTarray [2 * i + 1] * c - FFTarray [2 * i] * s);
-#else
switch (N) {case 256:
unscambled = unscambled64[i];
@@ -122,7 +101,6 @@
}
tempr = fac * (FFTarray [unscambled].re * c + FFTarray [unscambled].im * s);
tempi = fac * (FFTarray [unscambled].im * c - FFTarray [unscambled].re * s);
-#endif
/* fill in output values */
data [2 * i] = -tempr; /* first half even */
@@ -136,91 +114,3 @@
s = s * cfreq + cold * sfreq;
}
}
-
-#ifdef USE_ORIG_FFT /* OLD FFT */
-
-#define SWAP(a,b) tempr=a;a=b;b=tempr
-
-void FFT (double *data, int nn, int isign) {-/* Varient of Numerical Recipes code from off the internet. It takes nn
-interleaved complex input data samples in the array data and returns nn interleaved
-complex data samples in place where the output is the FFT of input if isign==1 and it
-is nn times the IFFT of the input if isign==-1.
-(Note: it doesn't renormalize by 1/N when doing the inverse transform!!!)
-(Note: this follows physicists convention of +i, not EE of -j in forward
-transform!!!!)
- */
-
-/* Press, Flannery, Teukolsky, Vettering "Numerical
- * Recipes in C" tuned up ; Code works only when nn is
- * a power of 2 */
-
- int n, mmax, m, j, i;
- double wtemp, wr, wpr, wpi, wi, theta, wpin;
- double tempr, tempi, datar, datai;
- double data1r, data1i;
-
- n = nn << 1;
-
-/* bit reversal */
-
- j = 0;
- for (i = 0; i < n; i += 2) {- if (j > i) { /* could use j>i+1 to help compiler analysis */- SWAP (data [j], data [i]);
- SWAP (data [j + 1], data [i + 1]);
- }
- m = nn;
- while (m >= 2 && j >= m) {- j -= m;
- m >>= 1;
- }
- j += m;
- }
-
-// theta = 3.141592653589795 * 0.5;
- theta = M_PI_2;
- if (isign < 0)
- theta = -theta;
- wpin = 0; /* sin(+-PI) */
- for (mmax = 2; n > mmax; mmax <<= 1) {- wpi = wpin;
- wpin = sin (theta);
- wpr = 1 - wpin * wpin - wpin * wpin; /* cos(theta*2) */
- theta *= .5;
- wr = 1;
- wi = 0;
- for (m = 0; m < mmax; m += 2) {- j = m + mmax;
- tempr = wr * (data1r = data [j]);
- tempi = wi * (data1i = data [j + 1]);
- for (i = m; i < n - mmax * 2; i += mmax * 2) {-/* mixed precision not significantly more
- * accurate here; if removing float casts,
- * tempr and tempi should be double */
- tempr -= tempi;
- tempi = wr *data1i + wi *data1r;
-/* don't expect compiler to analyze j > i + 1 */
- data1r = data [j + mmax * 2];
- data1i = data [j + mmax * 2 + 1];
- data [i] = (datar = data [i]) + tempr;
- data [i + 1] = (datai = data [i + 1]) + tempi;
- data [j] = datar - tempr;
- data [j + 1] = datai - tempi;
- tempr = wr *data1r;
- tempi = wi *data1i;
- j += mmax * 2;
- }
- tempr -= tempi;
- tempi = wr * data1i + wi *data1r;
- data [i] = (datar = data [i]) + tempr;
- data [i + 1] = (datai = data [i + 1]) + tempi;
- data [j] = datar - tempr;
- data [j + 1] = datai - tempi;
- wr = (wtemp = wr) * wpr - wi * wpi;
- wi = wtemp * wpi + wi * wpr;
- }
- }
-}
-
-#endif /* OLD FFT*/
\ No newline at end of file
--- a/transfo.h
+++ b/transfo.h
@@ -1,13 +1,9 @@
#ifndef TRANSFORM_H
#define TRANSFORM_H
-/* Define this to use original FFT*/
-//#define USE_ORIG_FFT
-
void MDCT(double* data, int N);
void FFT(double *data, int nn, int isign);
-#ifndef USE_ORIG_FFT
#define c_re(c) ((c).re)
#define c_im(c) ((c).im)
@@ -17,7 +13,6 @@
#define DEFINE_PFFTW(size) \
void pfftw_d_##size(fftw_complex_d *input); \
- void pfftwi_d_##size(fftw_complex_d *input); \
int pfftw_d_permutation_##size(int i);
DEFINE_PFFTW(16)
@@ -30,4 +25,3 @@
#endif
-#endif
\ No newline at end of file