Update of
/cvsroot/audacity/lib-src/libnyquist/nyquist/ffts/Numerical-Recipes-testing
In directory
23jxhf1.ch3.sourceforge.com:/tmp/cvs-serv8506/ffts/Numerical-Recipes-testing
Added Files:
fftTest.c fftTest2d.c fftTest3d.c rfftTest.c rfftTest2d.c
Log Message:
Updating to Nyquist v3.03.
--- NEW FILE: rfftTest2d.c ---
/* A program to test real 2d forward and inverse fast fourier transform
routines */
#include <NR.H> /* uses rlft3 from numerical recipes in C to verify rifft2d */
/*change fmin in numerical recipes to fminnr to avoid
conflict with fp.h */
#include <stdio.h>
#include <stdlib.h>
#include <fp.h>
#include <math.h>
#include "fftlib.h"
#include "fftext.h"
#include "fft2d.h"
#include <NRUTIL.H> // uses ugly tensors from numerical recipes; so can
call rlft3
#if macintosh
#include <timer.h>
#endif
#define NSIZES 24 /* the number of different ffts sizes to test */
#define BIPRAND(a) (2.0/(RAND_MAX+1.0)*a-1.0)
void main(){
long fftSize[NSIZES] = /* size of FFTs, must be powers of 2 */
{2, 4, 8, 16,
32, 64, 128, 256,
512, 1024, 2048, 4096, 8192, 16384, 32768, 65536,
131072, 262144, 524288, 1048576,
2097152, 4194304, 8388608, 16777216};
float *a;
long N2 = 64; /* the number of rows in 2d ffts, must be power of 2 */
long isize;
long i1;
long i2;
long TheErr;
long N;
long M;
long M2;
float maxerrifft;
float maxerrfft;
float ***NRtensdata; /* needed for rlft3 */
float **NRmatdata; /* needed for rlft3 */
float *specdata; /* needed for rlft3 */
float t1,t2;
unsigned int randseed = 777;
int rannum;
#if macintosh
UnsignedWide TheTime1;
Microseconds(&TheTime1);
randseed = TheTime1.lo;
#endif
printf(" %6d Byte Floats \n", sizeof(a[0]));
printf(" randseed = %10u\n", randseed);
for (isize = 0; isize < NSIZES; isize++){
srand(randseed);
N = fftSize[isize];
M = roundtol(LOG2(N));
N = POW2(M);
M2 = roundtol(LOG2(N2));
N2 = POW2(M2);
printf("rffts size = %6d X%6d, ", N2, N);
TheErr = 0;
TheErr = fft2dInit(M2, M);
if(!TheErr){
NRmatdata=matrix(1,1,1,2*N2);
specdata = &NRmatdata[1][1];
NRtensdata=f3tensor(1,1,1,N2,1,N); // uses ugly tensors
from NRUTIL; so can call rlft3
a = &NRtensdata[1][1][1];
if ((a == 0)||(specdata == 0)) TheErr = 2;
}
if(!TheErr){
/* set up a simple test case */
for (i1=0; i1<N2*N; i1++){
rannum = rand();
a[i1] = BIPRAND(rannum);
}
/* first use rlft3 from numerical recipes in C to
verify rifft2d */
/* unfortunately numerical recipes in C uses
backwards time and space */
/* forward fft, so our answer comes out time and space
reversed */
rlft3(NRtensdata, NRmatdata, 1, N2, N, 1);
/* move data to my in-place order */
a[1] = a[N2/2*N]; // pack in nyquest
wavenumber point for DC freq transform
for (i2=0; i2<N2/2; i2++){ // move transform of nyquest
frequency
a[(N2/2+i2)*N] = specdata[i2*2];
a[(N2/2+i2)*N+1] = specdata[i2*2+1];
}
a[N2/2*N+1] = specdata[N2]; // pack in nyquest wavenumber
point for nyquest freq transform
rifft2d(a, M2, M);
srand(randseed);
rannum = rand();
maxerrifft = fabs(BIPRAND(rannum)-a[0]);
for (i1=1; i1<N; i1++){
rannum = rand();
t1=BIPRAND(rannum);
t2=a[N-i1];
maxerrifft = fmax(maxerrifft,
fabs(BIPRAND(rannum)-a[N-i1]));
}
for (i2=1; i2<N2; i2++){
rannum = rand();
t1=BIPRAND(rannum);
t2=a[(N2-i2)*N];
maxerrifft = fmax(maxerrifft,
fabs(BIPRAND(rannum)-a[(N2-i2)*N]));
for (i1=1; i1<N; i1++){
rannum = rand();
t1=BIPRAND(rannum);
t2=a[(N2-i2)*N+N-i1];
maxerrifft = fmax(maxerrifft,
fabs(BIPRAND(rannum)-a[(N2-i2)*N+N-i1]));
}
}
printf("maxerrifft = %6.4e, ", maxerrifft);
/* now use rifft2d to verify rfft2d */
srand(randseed);
for (i1=0; i1<N2*N; i1++){
rannum = rand();
a[i1] = BIPRAND(rannum);
}
rifft2d(a, M2, M);
rfft2d(a, M2, M);
maxerrfft = 0;
srand(randseed);
for (i1=0; i1<N2*N; i1++){
rannum = rand();
maxerrfft = fmax(maxerrfft,
fabs(BIPRAND(rannum)-a[i1]));
}
printf("maxerrfft = %6.4e\n", maxerrfft);
fft2dFree();
free_f3tensor(NRtensdata,1,1,1,N2,1,N);
free_matrix(NRmatdata,1,1,1,2*N2);
a = 0;
}
else{
if(TheErr==2) printf(" out of memory \n");
else printf(" error \n");
fft2dFree();
}
}
printf(" Done. \n");
return;
}
--- NEW FILE: fftTest.c ---
/* A program to test complex forward and inverse fast fourier transform
routines */
#include <NR.H> /* uses four1 from numerical recipes in C to verify iffts */
/*change fmin in numerical recipes to fminnr to avoid
conflict with fp.h */
#include <stdio.h>
#include <stdlib.h>
#include <fp.h>
#include <math.h>
#include "fftlib.h"
#include "fftext.h"
#if macintosh
#include <timer.h>
#endif
#define NSIZES 24 /* the number of different fft sizes to test */
#define BIPRAND(a) (2.0/(RAND_MAX+1.0)*a-1.0)
typedef struct{
float Re;
float Im;
} Complex;
void main(){
long fftSize[NSIZES] = /* size of FFTs, must be powers of 2 */
{2, 4, 8, 16,
32, 64, 128, 256,
512, 1024, 2048, 4096, 8192, 16384, 32768, 65536,
131072, 262144, 524288, 1048576,
2097152, 4194304, 8388608, 16777216};
Complex *a1;
const long N2 = 2; /* the number ffts to test at each size */
long isize;
long i1;
long TheErr;
long N;
long M;
float maxerrifft;
float maxerrfft;
unsigned int randseed = 777;
int rannum;
#if macintosh
UnsignedWide TheTime1;
Microseconds(&TheTime1);
randseed = TheTime1.lo;
#endif
printf(" %6d Byte Floats \n", sizeof(a1[0].Re));
printf(" randseed = %10u\n", randseed);
for (isize = 0; isize < NSIZES; isize++){
srand(randseed);
N = fftSize[isize];
printf("ffts size = %8d, ", N);
M = roundtol(LOG2(N));
TheErr = fftInit(M);
if(!TheErr){
a1 = (Complex *) malloc( N2*N*sizeof(Complex) );
if (a1 == 0) TheErr = 2;
}
if(!TheErr){
/* set up a1 simple test case */
for (i1=0; i1<N2*N; i1++){
rannum = rand();
a1[i1].Re = BIPRAND(rannum);
rannum = rand();
a1[i1].Im = BIPRAND(rannum);
}
/* first use four1 from numerical recipes in C to
verify iffts */
/* Note their inverse fft is really the conventional
forward fft */
for (i1=0; i1<N2; i1++){
four1((float *)(a1+i1*N)-1, N, -1);
}
iffts((float *)a1, M, N2);
maxerrifft = 0;
srand(randseed);
for (i1=0; i1<N2*N; i1++){
rannum = rand();
maxerrifft = fmax(maxerrifft,
fabs(BIPRAND(rannum)-a1[i1].Re));
a1[i1].Re = BIPRAND(rannum);
rannum = rand();
maxerrifft = fmax(maxerrifft,
fabs(BIPRAND(rannum)-a1[i1].Im));
a1[i1].Im = BIPRAND(rannum);
}
printf("maxerrifft = %6.4e, ", maxerrifft);
/* now use iffts to verify ffts */
iffts((float *)a1, M, N2);
ffts((float *)a1, M, N2);
maxerrfft = 0;
srand(randseed);
for (i1=0; i1<N2*N; i1++){
rannum = rand();
maxerrfft = fmax(maxerrfft,
fabs(BIPRAND(rannum)-a1[i1].Re));
rannum = rand();
maxerrfft = fmax(maxerrfft,
fabs(BIPRAND(rannum)-a1[i1].Im));
}
printf("maxerrfft = %6.4e\n", maxerrfft);
free(a1);
fftFree();
}
else{
if(TheErr==2) printf(" out of memory \n");
else printf(" error \n");
fftFree();
}
}
printf(" Done. \n");
return;
}
--- NEW FILE: fftTest2d.c ---
/* A program to test 2d complex forward and inverse fast fourier transform
routines */
#include <NR.H> /* uses fourn from numerical recipes in C to verify ifft2d */
/*change fmin in numerical recipes to fminnr to avoid
conflict with fp.h */
#include <stdio.h>
#include <stdlib.h>
#include <fp.h>
#include <math.h>
#include "fftlib.h"
#include "fftext.h"
#include "fft2d.h"
#if macintosh
#include <timer.h>
#endif
#define NSIZES 24 /* the number of different ffts col sizes to
test */
#define BIPRAND(a) (2.0/(RAND_MAX+1.0)*a-1.0)
//#define BIPRAND(a) round(100*(2.0/(RAND_MAX+1.0)*a-1.0))
typedef struct{
float Re;
float Im;
} Complex;
void main(){
long fftSize[NSIZES] = /* size of FFTs cols, must be powers of 2 */
{2, 4, 8, 16,
32, 64, 128, 256,
512, 1024, 2048, 4096, 8192, 16384, 32768, 65536,
131072, 262144, 524288, 1048576,
2097152, 4194304, 8388608, 16777216};
Complex *a1;
long N2 = 64; /* the number of rows in the 2d fft */
long isize;
long i1;
long TheErr;
long N;
long M;
long M2;
float maxerrifft;
float maxerrfft;
unsigned long nn[2];
unsigned int randseed = 777;
int rannum;
#if macintosh
UnsignedWide TheTime1;
Microseconds(&TheTime1);
randseed = TheTime1.lo;
#endif
printf(" %6d Byte Floats \n", sizeof(a1[0].Re));
printf(" randseed = %10u\n", randseed);
for (isize = 0; isize < NSIZES; isize++){
srand(randseed);
N = fftSize[isize];
M = roundtol(LOG2(N));
N = POW2(M);
M2 = roundtol(LOG2(N2));
N2 = POW2(M2);
printf("ffts size = %6d X%6d, ", N2, N);
nn[0] = N2;
nn[1] = N;
TheErr = fft2dInit(M2, M);
if(!TheErr){
a1 = (Complex *) malloc(N2*N*sizeof(Complex) );
if (a1 == 0) TheErr = 2;
}
if(!TheErr){
/* set up a simple test case */
for (i1=0; i1<N2*N; i1++){
rannum = rand();
a1[i1].Re = BIPRAND(rannum);
rannum = rand();
a1[i1].Im = BIPRAND(rannum);
}
/* first use fourn from numerical recipes in C to
verify ifft2d */
/* Note their inverse fft is really the conventional
forward fft */
fourn((float *)a1-1, nn-1, 2, -1);
ifft2d((float *)a1, M2, M);
maxerrifft = 0;
srand(randseed);
for (i1=0; i1<N2*N; i1++){
rannum = rand();
maxerrifft = fmax(maxerrifft,
fabs(BIPRAND(rannum)-a1[i1].Re));
a1[i1].Re = BIPRAND(rannum);
rannum = rand();
maxerrifft = fmax(maxerrifft,
fabs(BIPRAND(rannum)-a1[i1].Im));
a1[i1].Im = BIPRAND(rannum);
}
printf("maxerrifft = %6.4e, ", maxerrifft);
/* now use iffts to verify ffts */
ifft2d((float *)a1, M2, M);
fft2d((float *)a1, M2, M);
maxerrfft = 0;
srand(randseed);
for (i1=0; i1<N2*N; i1++){
rannum = rand();
maxerrfft = fmax(maxerrfft,
fabs(BIPRAND(rannum)-a1[i1].Re));
rannum = rand();
maxerrfft = fmax(maxerrfft,
fabs(BIPRAND(rannum)-a1[i1].Im));
}
printf("maxerrfft = %6.4e\n", maxerrfft);
free(a1);
fft2dFree();
}
else{
if(TheErr==2) printf(" out of memory \n");
else printf(" error \n");
fft2dFree();
}
}
printf(" Done. \n");
return;
}
--- NEW FILE: rfftTest.c ---
/* A program to test real forward and inverse fast fourier transform routines
*/
#include <NR.H> /* uses realft from numerical recipes in C to verify riffts */
/*change fmin in numerical recipes to fminnr to avoid
conflict with fp.h */
#include <stdio.h>
#include <stdlib.h>
#include <fp.h>
#include <math.h>
#include "fftlib.h"
#include "fftext.h"
#if macintosh
#include <timer.h>
#endif
#define NSIZES 24 /* the number of different fft sizes to test */
#define BIPRAND(a) (2.0/(RAND_MAX+1.0)*a-1.0)
void main(){
long fftSize[NSIZES] = /* size of FFTs, must be powers of 2 */
{2, 4, 8, 16,
32, 64, 128, 256,
512, 1024, 2048, 4096, 8192, 16384, 32768, 65536,
131072, 262144, 524288, 1048576,
2097152, 4194304, 8388608, 16777216};
float *a;
const long N2 = 2; /* the number ffts to test at each size */
long isize;
long i1;
long i2;
long TheErr;
long N;
long M;
float maxerrifft;
float maxerrfft;
unsigned int randseed = 777;
int rannum;
#if macintosh
UnsignedWide TheTime1;
Microseconds(&TheTime1);
randseed = TheTime1.lo;
#endif
printf(" %6d Byte Floats \n", sizeof(a[0]));
printf(" randseed = %10u\n", randseed);
for (isize = 0; isize < NSIZES; isize++){
srand(randseed);
N = fftSize[isize];
printf("rffts size = %8d, ", N);
M = roundtol(LOG2(N));
TheErr = 0;
TheErr = fftInit(M);
if(!TheErr){
a = (float *) malloc(N2*N*sizeof(float) );
if (a == 0) TheErr = 2;
}
if(!TheErr){
/* set up a simple test case */
for (i1=0; i1<N2*N; i1++){
rannum = rand();
a[i1] = BIPRAND(rannum);
}
/* first use realft from numerical recipes in C to
verify riffts */
/* unfortunately numerical recipes in C uses
backwards time */
/* forward fft, so our answer comes out time reversed
*/
for (i2=0; i2<N2; i2++){
realft((a+i2*N)-1, N, 1);
}
riffts(a, M, N2);
srand(randseed);
for (i2=0; i2<N2; i2++){
rannum = rand();
maxerrifft = fabs(BIPRAND(rannum)-a[i2*N]);
for (i1=1; i1<N; i1++){
rannum = rand();
maxerrifft = fmax(maxerrifft,
fabs(BIPRAND(rannum)-a[i2*N+N-i1]));
}
}
printf("maxerrifft = %6.4e, ", maxerrifft);
/* now use iffts to verify ffts */
srand(randseed);
for (i1=0; i1<N2*N; i1++){
rannum = rand();
a[i1] = BIPRAND(rannum);
}
riffts(a, M, N2);
rffts(a, M, N2);
maxerrfft = 0;
srand(randseed);
for (i1=0; i1<N2*N; i1++){
rannum = rand();
maxerrfft = fmax(maxerrfft,
fabs(BIPRAND(rannum)-a[i1]));
}
printf("maxerrfft = %6.4e\n", maxerrfft);
fftFree();
free(a);
a = 0;
}
else{
if(TheErr==2) printf(" out of memory \n");
else printf(" error \n");
fftFree();
}
}
printf(" Done. \n");
return;
}
--- NEW FILE: fftTest3d.c ---
/* A program to test 3d complex forward and inverse fast fourier transform
routines */
#include <NR.H> /* uses fourn from numerical recipes in C to verify ifft3d */
/*change fmin in numerical recipes to fminnr to avoid
conflict with fp.h */
#include <stdio.h>
#include <stdlib.h>
#include <fp.h>
#include <math.h>
#include "fftlib.h"
#include "fftext.h"
#include "fft2d.h"
#if macintosh
#include <timer.h>
#endif
#define NSIZES 24 /* the number of different ffts col sizes to
test */
#define BIPRAND(a) (2.0/(RAND_MAX+1.0)*a-1.0)
typedef struct{
float Re;
float Im;
} Complex;
void main(){
long fftSize[NSIZES] = /* size of FFTs cols, must be powers of 2 */
{2, 4, 8, 16,
32, 64, 128, 256,
512, 1024, 2048, 4096, 8192, 16384, 32768, 65536,
131072, 262144, 524288, 1048576,
2097152, 4194304, 8388608, 16777216};
Complex *a1;
long isize;
long i1;
long TheErr;
long N;
long M;
long N2 = 16; /* the number of rows in the 3d fft */
long M2;
long N3 = 32; /* the number of pages in the 3d fft */
long M3;
float maxerrifft;
float maxerrfft;
unsigned long nn[3];
unsigned int randseed = 777;
int rannum;
#if macintosh
UnsignedWide TheTime1;
Microseconds(&TheTime1);
randseed = TheTime1.lo;
#endif
printf(" %6d Byte Floats \n", sizeof(a1[0].Re));
printf(" randseed = %10u\n", randseed);
for (isize = 0; isize < NSIZES; isize++){
srand(randseed);
N = fftSize[isize];
M = roundtol(LOG2(N));
N = POW2(M);
M2 = roundtol(LOG2(N2));
N2 = POW2(M2);
M3 = roundtol(LOG2(N3));
N3 = POW2(M3);
printf("ffts size = %5d X%5d X%6d, ", N3, N2, N);
nn[0] = N3;
nn[1] = N2;
nn[2] = N;
TheErr = fft3dInit(M3, M2, M);
if(!TheErr){
a1 = (Complex *) malloc(N3*N2*N*sizeof(Complex) );
if (a1 == 0) TheErr = 2;
}
if(!TheErr){
/* set up a simple test case */
for (i1=0; i1<N3*N2*N; i1++){
rannum = rand();
a1[i1].Re = BIPRAND(rannum);
rannum = rand();
a1[i1].Im = BIPRAND(rannum);
}
/* first use fourn from numerical recipes in C to
verify ifft3d */
/* Note their inverse fft is really the conventional
forward fft */
fourn((float *)a1-1, nn-1, 3, -1);
ifft3d((float *)a1, M3, M2, M);
maxerrifft = 0;
srand(randseed);
for (i1=0; i1<N3*N2*N; i1++){
rannum = rand();
maxerrifft = fmax(maxerrifft,
fabs(BIPRAND(rannum)-a1[i1].Re));
a1[i1].Re = BIPRAND(rannum);
rannum = rand();
maxerrifft = fmax(maxerrifft,
fabs(BIPRAND(rannum)-a1[i1].Im));
a1[i1].Im = BIPRAND(rannum);
}
printf("errifft = %4.3e, ", maxerrifft);
/* now use iffts to verify ffts */
ifft3d((float *)a1, M3, M2, M);
fft3d((float *)a1, M3, M2, M);
maxerrfft = 0;
srand(randseed);
for (i1=0; i1<N3*N2*N; i1++){
rannum = rand();
maxerrfft = fmax(maxerrfft,
fabs(BIPRAND(rannum)-a1[i1].Re));
rannum = rand();
maxerrfft = fmax(maxerrfft,
fabs(BIPRAND(rannum)-a1[i1].Im));
}
printf("errfft = %4.3e\n", maxerrfft);
free(a1);
fft3dFree();
}
else{
if(TheErr==2) printf(" out of memory \n");
else printf(" error \n");
fft3dFree();
}
}
printf(" Done. \n");
return;
}
------------------------------------------------------------------------------
Open Source Business Conference (OSBC), March 24-25, 2009, San Francisco, CA
-OSBC tackles the biggest issue in open source: Open Sourcing the Enterprise
-Strategies to boost innovation and cut costs with open source participation
-Receive a $600 discount off the registration fee with the source code: SFAD
http://p.sf.net/sfu/XcvMzF8H
_______________________________________________
Audacity-cvs mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/audacity-cvs
