Commit 1eb5908b by kollo

Make use of fftw library if present

The internal fft algorithm has some strange flaws, but I keept it
for platforms without libfftw
Interface stays the same
parent d4ba5671
......@@ -8,3 +8,4 @@ Release notes for Version 1.25 (Sept 2016 -- 2017)
- fixed bug in EXIT (bytecode)
- added bluetooth support (linux)
- improved USB functions (linux)
- make use of fftw library when available
......@@ -139,6 +139,12 @@
/* mpz_primorial_ui in GMP */
#undef HAVE_PRIMORIAL_UI
/* fftw_make_plan in FFTW */
#undef HAVE_FFTW
/* rfftw_create_plan in RFFTW */
#undef HAVE_RFFTW
/* _SystemCycle in dynamical TINE lib */
#undef HAVE_TINE
......
......@@ -5668,6 +5668,106 @@ $as_echo "#define HAVE_PRIMORIAL_UI 1" >>confdefs.h
fi
# `-lfftw':
{ $as_echo "$as_me:${as_lineno-$LINENO}: checking for fftw_make_plan in -lfftw" >&5
$as_echo_n "checking for fftw_make_plan in -lfftw... " >&6; }
if ${ac_cv_lib_fftw_fftw_make_plan+:} false; then :
$as_echo_n "(cached) " >&6
else
ac_check_lib_save_LIBS=$LIBS
LIBS="-lfftw $libs $LIBS"
cat confdefs.h - <<_ACEOF >conftest.$ac_ext
/* end confdefs.h. */
/* Override any GCC internal prototype to avoid an error.
Use char because int might match the return type of a GCC
builtin and then its argument prototype would still apply. */
#ifdef __cplusplus
extern "C"
#endif
char fftw_make_plan ();
int
main ()
{
return fftw_make_plan ();
;
return 0;
}
_ACEOF
if ac_fn_c_try_link "$LINENO"; then :
ac_cv_lib_fftw_fftw_make_plan=yes
else
ac_cv_lib_fftw_fftw_make_plan=no
fi
rm -f core conftest.err conftest.$ac_objext \
conftest$ac_exeext conftest.$ac_ext
LIBS=$ac_check_lib_save_LIBS
fi
{ $as_echo "$as_me:${as_lineno-$LINENO}: result: $ac_cv_lib_fftw_fftw_make_plan" >&5
$as_echo "$ac_cv_lib_fftw_fftw_make_plan" >&6; }
if test "x$ac_cv_lib_fftw_fftw_make_plan" = xyes; then :
libs="$libs -lfftw"
$as_echo "#define HAVE_FFTW 1" >>confdefs.h
fi
# `-lrfftw':
{ $as_echo "$as_me:${as_lineno-$LINENO}: checking for rfftw_create_plan in -lrfftw" >&5
$as_echo_n "checking for rfftw_create_plan in -lrfftw... " >&6; }
if ${ac_cv_lib_rfftw_rfftw_create_plan+:} false; then :
$as_echo_n "(cached) " >&6
else
ac_check_lib_save_LIBS=$LIBS
LIBS="-lrfftw $libs $LIBS"
cat confdefs.h - <<_ACEOF >conftest.$ac_ext
/* end confdefs.h. */
/* Override any GCC internal prototype to avoid an error.
Use char because int might match the return type of a GCC
builtin and then its argument prototype would still apply. */
#ifdef __cplusplus
extern "C"
#endif
char rfftw_create_plan ();
int
main ()
{
return rfftw_create_plan ();
;
return 0;
}
_ACEOF
if ac_fn_c_try_link "$LINENO"; then :
ac_cv_lib_rfftw_rfftw_create_plan=yes
else
ac_cv_lib_rfftw_rfftw_create_plan=no
fi
rm -f core conftest.err conftest.$ac_objext \
conftest$ac_exeext conftest.$ac_ext
LIBS=$ac_check_lib_save_LIBS
fi
{ $as_echo "$as_me:${as_lineno-$LINENO}: result: $ac_cv_lib_rfftw_rfftw_create_plan" >&5
$as_echo "$ac_cv_lib_rfftw_rfftw_create_plan" >&6; }
if test "x$ac_cv_lib_rfftw_rfftw_create_plan" = xyes; then :
libs="-lrfftw $libs"
$as_echo "#define HAVE_RFFTW 1" >>confdefs.h
fi
if test "n$with_usb" = nyes; then
# `-lusb':
{ $as_echo "$as_me:${as_lineno-$LINENO}: checking for usb_open in -lusb" >&5
......
......@@ -274,6 +274,24 @@ libs="$libs -lgmp"
AC_DEFINE(HAVE_PRIMORIAL_UI,1,[mpz_primorial_ui in gmp])
],,[$libs])
# `-lfftw':
AC_CHECK_LIB([fftw], [fftw_make_plan], [
libs="$libs -lfftw"
AC_DEFINE(HAVE_FFTW,1,[fftw_make_plan in fftw])
],,[$libs])
# `-lrfftw':
AC_CHECK_LIB([rfftw], [rfftw_create_plan], [
libs="-lrfftw $libs"
AC_DEFINE(HAVE_RFFTW,1,[rfftw_create_plan in rfftw])
],,[$libs])
if test "n$with_usb" = nyes; then
# `-lusb':
AC_CHECK_LIB([usb], [usb_open], [
......
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "config.h"
#ifdef HAVE_RFFTW
#include <rfftw.h>
#include "defs.h"
#include "x11basic.h"
#include "mathematics.h"
#include "fft.h"
void fftw_realft(int N, fftw_real *in, fftw_real *out, int isign) {
rfftw_plan p;
if(isign) p=rfftw_create_plan(N, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
else p=rfftw_create_plan(N, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
rfftw_one(p, in, out);
rfftw_destroy_plan(p);
}
#else
/********************************************************************
* *
* THIS FILE IS PART OF THE OggSQUISH SOFTWARE CODEC SOURCE CODE. *
......@@ -2237,3 +2260,4 @@ wsave contains initialization calculations which must not be
["send index for vfftpk" describes a vectorized version of fftpack]
*/
#endif
......@@ -12,6 +12,10 @@
#define SQRT2 1.4142135623730950488016887242097
#define TSQRT2 2.8284271247461900976033774484194
void rffti(int n, FLOAT *wsave, int *ifac);
void rfftf(int n,FLOAT *r,FLOAT *wsave,int *ifac);
void rfftb(int n, FLOAT *r, FLOAT *wsave, int *ifac);
#ifdef HAVE_RFFTW
void fftw_realft(int N, double *in, double *out, int isign);
#else
void rffti(int n, FLOAT *wsave, int *ifac);
void rfftf(int n,FLOAT *r,FLOAT *wsave,int *ifac);
void rfftb(int n, FLOAT *r, FLOAT *wsave, int *ifac);
#endif
......@@ -61,37 +61,43 @@ double determinante(double *a,int n) {
/* Routinen f"uer Fouriertransformation */
void realft(double *data, int n,int isign) {
if(n<=0) return;
#ifdef HAVE_RFFTW
double *spek=malloc(n*sizeof(double));
int j;
if(isign) {
spek[0]=data[0];
spek[n/2]=data[n-1];
for(j=1;j<(n+1)/2;j++) {
spek[j]=data[2*j-1];
spek[n-j]=data[2*j];
}
fftw_realft(n,spek,data,isign);
} else {
fftw_realft(n,data,spek,isign);
data[0]=spek[0];
data[n-1]=spek[n/2];
for(j=1;j<(n+1)/2;j++) {
data[2*j-1]=spek[j];
data[2*j]=spek[n-j];
}
}
free(spek);
#else
static int oldn=-1;
static double *wsave=NULL;
static int ifac[15];
if(n<=0) return;
if(n>1 && n!=oldn) {
wsave=realloc(wsave,(2*n+15)*sizeof(double));
rffti(n,wsave,ifac);
#ifdef DEBUG
printf("Initialized for n=%d:\n",n);
if(n>1 && ifac[1]>0) {
int i;
for(i=2;i<2+ifac[1];i++) {
printf("ifac[%d]=%d\n",i,ifac[i]);
}
}
#endif
oldn=n;
}
if(isign) {
#ifdef DEBUG
int i;
for(i=0;i<10;i++) printf("data[%d]=%g\n",i,data[i]);
#endif
rfftf(n,data,wsave,ifac);
} else rfftb(n,data,wsave,ifac);
if(isign) rfftf(n,data,wsave,ifac);
else rfftb(n,data,wsave,ifac);
#endif
}
#if 0
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment