1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
|
/* 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;
}
|
