ref: 0876f43acf4ac42f7bd3730a7cbed996cdc4eee6
parent: d79e7721a56ebc0237d8bf8aff7f2ae90d63fb58
author: oxygene2000 <oxygene2000>
date: Wed Feb 9 10:44:41 EST 2000
Added additional implementation of rdft using fftw (from rank 1 to rank >20 in my profiler = a lot faster)
--- a/encoder.c
+++ b/encoder.c
@@ -13,7 +13,7 @@
#include "rateconv.h"
#include "enc.h" /* encoder cores */
-
+#include "rdft.h"
#define BITHEADERBUFSIZE (65536)
/* ---------- functions ---------- */
@@ -586,6 +586,7 @@
printf("Cut-off frequency: %dHz.\n", cut_off); printf("\n");+ fftw_init();
for (i = 0; i < FileCount; i++) {char aac_fn[255];
char *fnp;
@@ -717,7 +718,7 @@
nSecs = nTotSecs - (60*nMins);
printf("Encoding %s took:\t%d:%.2d\t\n", FileNames[i], nMins, nSecs);}
-
+ fftw_destroy();
return FNO_ERROR;
}
--- /dev/null
+++ b/rdft.h
@@ -1,0 +1,6 @@
+#include <rfftw.h>
+rfftw_plan rdft_plan11;
+rfftw_plan rdft_plan8;
+
+void fftw_init();
+void fftw_destroy();
--- a/rdft_spectrum.c
+++ b/rdft_spectrum.c
@@ -33,6 +33,9 @@
#define CSQR4(a, b, u, v) { u=a*a-b*b; v=a*b; v+=v; } #define CMULT6(c, s, c1, s1, u, v) { u=c1*c-s1*s; v=c1*s+s1*c; }+void fftw_init(){}+void fftw_destroy(){}+
void rdft( double *fr, unsigned lg2n )
/* real discrete Fourier transform; not recursive */
{@@ -153,6 +156,7 @@
fi = fr + i;
gi = fr + k1 - i;
+//most timeconsuming part of rdft (tuning?)
do
{double a, b, g0, f0, f1, g1, f2, g2, f3, g3;
--- /dev/null
+++ b/rdft_spectrum2.c
@@ -1,0 +1,112 @@
+
+/* rdft_spectrum.c by James Salsman 25 April 1999 for public domain */
+
+/* spectrum( double *f, unsigned lg2n ):
+ input f[0...2^lg2n-1] amplitudes, modified in-place;
+ output f[0...2^(lg2n-1)-1] magnitudes (scaling depends on lg2n;
+ f[0] ill-defined);
+ complspectrum( double *f, unsigned lg2n ): same as above but
+ f[2^(lg2n-1)+1...2^lg2n-1] are also phases in radians
+ (f[2^(lg2n-1)] meaningless.) */
+
+/* computes a power spectrum using a real discrete Fourier transform,
+ faster than a full FFT but not exactly a Hartley transform. */
+
+/* adapted from Joerg Arndt's code found around
+ http://www.spektracom.de/~arndt/joerg.html
+ in turn from Ron Mayer's 1988 Stanford project, and Ron credits
+ Bracewell, Hartley, Gauss, and Euler. Don't forget Pappus,
+ Pythagoras, and Euclid! */
+
+#include <math.h> /* we need sqrt(), sin(), cos(), and atan2(),
+ and we want the sincos(x,*cos) FP primitive, but if we can't
+ have that, we should use floor() and fmod(,) instead of cos() */
+
+#include "rdft.h"
+#include <rfftw.h>
+
+
+void fftw_init(){+rdft_plan11=rfftw_create_plan(11,FFTW_REAL_TO_COMPLEX,FFTW_MEASURE);
+rdft_plan8=rfftw_create_plan(8,FFTW_REAL_TO_COMPLEX,FFTW_MEASURE);
+}
+
+void fftw_destroy(){+rfftw_destroy_plan(rdft_plan11);
+rfftw_destroy_plan(rdft_plan8);
+}
+
+void rdft( double *fr, unsigned lg2n )
+/* real discrete Fourier transform; not recursive */
+{+rfftw_plan rdft_plan;
+double fo[lg2n];
+switch(lg2n) {+ case 11: rfftw_one(rdft_plan11,fr,fo);
+ break;
+ case 8: rfftw_one(rdft_plan8,fr,fo);
+ break;
+ default: rdft_plan=rfftw_create_plan(lg2n,FFTW_REAL_TO_COMPLEX,
+ FFTW_ESTIMATE);
+ rfftw_one(rdft_plan,fr,fo);
+ rfftw_destroy_plan(rdft_plan);
+ printf("ERROR: rdft with size %i",lg2n);+}
+ memcpy(fr,fo,sizeof(fr));
+}
+
+void spectrum( double *f, unsigned lg2n )
+/* magnitude power spectrum from rdft */
+{+ const int n=(1<<lg2n), k=(n>>1); /* k=n/2 */
+ int i, j;
+
+ rdft(f, lg2n);
+
+ for ( i=1,j=n-1; i<k; i++,j-- )
+ {+ double a, b;
+
+ a = f[i] + f[j]; /* real part */
+ b = f[i] - f[j]; /* imag part */
+
+ f[i] = sqrt(a*a + b*b); /* magnitude -- please note that the
+ scaling depends on size n */
+ }
+
+ f[0]=sqrt(f[0]*f[0]+f[n/2]*f[n/2]); /* ill-defined bin */
+ /* Reciprocated division by zero precludes any meaning
+ of "0 Hz" so please avoid using the f[0] output! */
+}
+
+void complspectrum( double *f, unsigned lg2n )
+/* polar complex power spectrum from rdft */
+{+ const int n=(1<<lg2n), k=(n>>1); /* k=n/2 */
+ int i, j;
+
+ rdft(f, lg2n);
+
+ for ( i=1,j=n-1; i<k; i++,j-- )
+ {+ double a, b;
+
+ a = f[i] + f[j]; /* real part */
+ b = f[i] - f[j]; /* imag part */
+
+ f[i] = sqrt(a*a + b*b); /* magnitude -- please note that the
+ scaling depends on size n */
+ /* complex part: phase */
+
+ if (a != 0.0) f[j] = atan2(b, a); /* magnitude f[i] has phase f[n-i] */
+ else f[j] = 0.0;
+ }
+
+ f[0]=sqrt(f[0]*f[0]+f[n/2]*f[n/2]); /* ill-defined bin */
+
+ /* The corresponding phase, f[k], should really be set to
+ some kind of not-a-number value because it is completely
+ meaningless, instead of just tainted like f[0]. */
+}
+
+/* end of rdft_spectrum.c :jps 26 April 1999 */#include <rfftw.h>