WHY:
Recently I was asked to Fourier analyze recorded acoustical data.
I have programs that handle data which does not exceed two 64K
segments and one program for extremely long series that will not
fit in RAM (conventional plus extended). The data I had to analize
fell between this two extremes and I suspected it could be handled
easily using huge pointers.
DATA SIZE RESTRICTION:
The series has to fit in conventional memory (ie. far heap). Using
double precision it means somewhere around 70.000 real data points
(547 K bytes), depending on your system configuration.
If your data exceeds this limit my program TOGAHOCK.PAS might
be useful or consult the following reference:
Performing Fourier transforms on extremely long data streams
W. K. Hocking
Computers in Physics- JAN FEB 1989.
TRUNCATION ERRORS ARE A BIG PROBLEM:
First I used two FFT subroutines found in SYMTEL and modified them
using huge pointers to access the data. Everything worked fine until
I started transforming series 12K long and above. It took me a while
to figure out the problem was in the truncation errors when calculating
the sines and cosines using the library functions.
METHOD:
I translated to C, R. C. Singleton's mixed radix fast Fourier transform
algorithm (see reference in SING.C). His algorithm generates the
sines and cosines recursively and corrects for truncation errors.
DATA LENGTH
Does not have to be a power of 2 necessarily. The length can contain
even factors as 2 and 4, and also odd factors as 3,5, 7, 11, etc.
The algorithm is most efficient if the length is a power of four.
Data lengths with odd factors of 3 and 5 can be used without a great
loss in performance.
USAGE TO OBTAIN THE DIRECT TRANSFORM:
In the calling program:
1-) allocate memory in the far heap for the real and imaginary
parts (ie. using farcalloc)
2-) equate two huge pointers to the beginning of the real and
imaginary parts, respectively
3-) read data into real and imaginary parts
4-) call
sing(huge_points_to_real, huge_points_to_imaginary,
order, order, order, -1);
5-) divide output by order.
USAGE TO OBTAIN THE DIRECT TRANSFORM OF REAL DATA USING THE DOUBLING
ALGORITHM: (total_length = 2*half_length)
In the calling program:
1-) allocate memory in the far heap for the real and imaginary
parts (ie. using farcalloc), each of length (half_length +1)
2-) equate two huge pointers to the beginning of the real and
imaginary parts, respectively
3-) read real data alternatively into real and imaginary arrays:
real_array contains r(0), r(2), r(4), ....
and imaginary_array r(1), r(3), r(5).......
Zero real_array(half_length)and Imaginary_array(half_length)
4-) call
sing(huge_points_to_real, huge_points_to_imaginary,
half_order,half_order,half_order, -1);
realtr(huge_points_to_real, huge_points_to_imaginary,
half_order, -1);
5-) divide output by 4*half_order.
USAGE TO OBTAIN THE INVERSE TRANSFORM:
In the calling program:
1-) allocate memory in the far heap for the real and imaginary
parts (ie. using farcalloc)
2-) equate two huge pointers to the beginning of the real and
imaginary parts spectra
4-) call
sing(huge_points_to_real, huge_points_to_imaginary,
order, order, order, 1);
USAGE TO OBTAIN THE INVERSE TRANSFORM WHEN ONLY THE SIMMETRICAL
HALF IS AVAILABLE: (ie, transform of real data using the doubling
algorithm)
In the calling program:
1-) allocate memory in the far heap for the real and imaginary
parts (ie. using farcalloc)
2-) equate two huge pointers to the beginning of the real and
imaginary parts spectra
4-) call
realtr(huge_points_to_real, huge_points_to_imaginary,
half_order, +1);
sing(huge_points_to_real, huge_points_to_imaginary,
order, order, order, 1);
5-) divide output by 4*half_order.
6-) the real time series is alternatively in the real and
imaginary arrays:
real_array contains r(0), r(2), r(4), ....
and imaginary_array r(1), r(3), r(5).......
PERFORMANCE :
Using the doubling algorithm and an AT with 80287 math coprocessor
Total Length Half Length Approx time in seconds
2048 1024 2
6144 3072 8
8192 4096 11
10240 5120 15
16384 8192 24
20480 10240 33
24576 12288 37
32768 16384 48
40960 20480 64
51200 25600 84
61440 30720 112
65536 32768 109
70000 35000 144
VERIFICATION:
Versing.c uses the doubling algorithm. This program generates
a rectangular pulse, calculates its fourier transform and compares
it with the known analytic closed form expression. It then
calculates the inverse transform to get back the original
rectangular pulse. This is a practical way to asses the
propagation of truncation errors.
FINAL COMMENT:
Constructive criticisms and suggestions are welcomed. I am
willing to adapt this programs to your special needs as long as
I can find free time to do it.