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.