samer@0: /*
samer@0:  *	Copyright (c) 2000, Samer Abdallah, King's College London.
samer@0:  *	All rights reserved.
samer@0:  *
samer@0:  *	This software is provided AS iS and WITHOUT ANY WARRANTY; 
samer@0:  *	without even the implied warranty of MERCHANTABILITY or 
samer@0:  *	FITNESS FOR A PARTICULAR PURPOSE.
samer@0:  */
samer@0: 
samer@0: package samer.units;
samer@0: import  samer.maths.*; // uses ComplexVector and Vec
samer@0: 
samer@0: /**
samer@0:   *	<h2>Fast Fourier Transform</h2>
samer@0:   *	<h4>by Samer Abdallah</h4>
samer@0:   *	<h4>15 December 1999</h4>
samer@0:   *
samer@0:   *	<p>
samer@0:   *	Mostly from my old C++ version, though I should credit the following
samer@0:   *	FFT code for the idea of using a lookup table for the bit reversal
samer@0:   *   and also a complete lookup table of sines:
samer@0:   *
samer@0:   *	<pre>
samer@0:   *	>	From the Unix Version 2.4 by Steve Sampson, Public Domain,
samer@0:   *	>	September 1988.
samer@0:   *	>	Adapted for Java by Ben Stoltz <stoltz@sun.com>, September 1997
samer@0:   *	>
samer@0:   *	>	Refer to http://www.neato.org/~ben/Fft.html for updates and related
samer@0:   *	>	resources.
samer@0:   *	</pre>
samer@0:   *
samer@0:   *	<p>
samer@0:   *	My old version computed the bit-reversed indices on the fly - silly
samer@0:   *	really - don't know why I didn't think of a lookup table myself.
samer@0:   *	Also, it had a lookup table of sines and cosines of the form
samer@0:   *	sin( 2*PI/k) with k going 1, 2, 4, 8 .. N. This version has complete
samer@0:   *	tables for sin and cosine - which is a bit of a waste of memory,
samer@0:   *	but frankly, who cares?
samer@0:   *
samer@0:   *	<h3>Usage</h3>
samer@0:   *	<p>
samer@0:   *	It is the user's responsibility to fill the arrays real
samer@0:   *   and imag with data IN THE CORRECT BITREVERSED ORDER. This
samer@0:   *	is easy enough using the bitrev lookup table - something like this
samer@0:   *	would do the trick:
samer@0:   *	<pre>
samer@0:   *		for (int i=0; i<N; i++) {
samer@0:   *			int k=fft.bitrev[i]
samer@0:   *			fft.real[k] = theRealData[i];
samer@0:   *			fft.imag[k] = 0;
samer@0:   *		}				
samer@0:   *	</pre>
samer@0:   *	The arrays are used destructively, so it is also the user's 
samer@0:   *	responsibility, eg, to reset the imaginary parts to zero even
samer@0:   *	if they are always zero.
samer@0:   *
samer@0:   *	<p>
samer@0:   *	The results ends up in the arrays real and imag
samer@0:   */
samer@0: 
samer@0: 
samer@0: public class FFT 
samer@0: {
samer@0: 	public	int				bitrev[];	// lookup table
samer@0: 	public	double			H[];			// window - not yet used internally
samer@0: 	public	double			real[];			
samer@0: 	public	double			imag[];			
samer@0: 
samer@0: 	protected	ComplexVector	W;				// sines and cosines lookup table
samer@0: 	protected	int				N;
samer@0: 	protected	double			zeros[];
samer@0: 
samer@0: 	public FFT( int n)
samer@0: 	{
samer@0: 		N=n;
samer@0: 		real = new double[N];
samer@0: 		imag = new double[N];
samer@0: 		W = new ComplexVector(N);
samer@0: 		bitrev = new int[N];
samer@0: 		H = new double[N];
samer@0: 
samer@0: 		// calculate bit reversal lookup table
samer@0: 		int m=N/2, j=0;
samer@0: 		for (int i=0; i<n; i++) {
samer@0: 			bitrev[i]=j;
samer@0: 			int b;
samer@0: 			for (b=m; (b>1) && (j>=b); b>>=1) j-=b;
samer@0: 			j+=b; 
samer@0: 		}
samer@0: 
samer@0: 		// calculate sines and cosines lookup
samer@0: 		double th=2*Math.PI/N;
samer@0: 		for (int i=0; i<n; i++) {
samer@0: 			W.real[i] = Math.cos(i*th); 
samer@0: 			W.imag[i] = Math.sin(i*th); 
samer@0: 		}
samer@0: 
samer@0: 		// table of zeros for fast initialisation
samer@0: 		zeros = new double[N];
samer@0: 		for (int i=0; i<N; i++) zeros[i]=0.0;
samer@0: 	}
samer@0: 
samer@0: 	public int size() { return N; }
samer@0: 	public final void input( double Y[])
samer@0: 	{
samer@0: 		System.arraycopy(zeros,0,imag,0,N);
samer@0: 		for (int i=0; i<N; i++) {
samer@0: 			real[bitrev[i]] = H[i]*Y[i];
samer@0: 		}
samer@0: 	}				
samer@0: 
samer@0: 	public final void brevcopy( double src[], double dst[] )
samer@0: 	{
samer@0: 		for (int i=0; i<src.length; i++) {
samer@0: 			dst[bitrev[i]] = src[i];
samer@0: 		}
samer@0: 	}				
samer@0: 
samer@0: 	public final void input( Vec.Iterator it)
samer@0: 	{
samer@0: 		System.arraycopy(zeros,0,imag,0,N);
samer@0: 		for (int i=0; i<N; i++) {
samer@0: 			real[bitrev[i]] = H[i]*it.next();
samer@0: 		}
samer@0: 	}				
samer@0: 
samer@0: 	public final void outputPower( double S[])
samer@0: 	{
samer@0: 		// extract power
samer@0: 		for (int i=0; i<S.length; i++) {
samer@0: 			S[i] = real[i]*real[i] + imag[i]*imag[i];
samer@0: 		}
samer@0: 	}
samer@0: 
samer@0: 	public final void calculate()
samer@0: 	{
samer@0: 		double	a1,a2,b1,b2;
samer@0: 		int		step, jump, i, k;
samer@0: 		int		m=N/2, p;
samer@0: 		double	wReal, wImag;
samer@0: 
samer@0: 		// these are cached for performance -
samer@0: 		// the COMPILER ought to be doing this!
samer@0: 		double	Real[] = real;
samer@0: 		double	Imag[] = imag;
samer@0: 		double	wr[] = W.real;
samer@0: 		double	wi[] = W.imag;
samer@0: 
samer@0: 		step=1; jump=2;
samer@0: 		while (step<N) {
samer@0: 
samer@0: 			p = 0; 
samer@0: 			for (k=0; k<step; k++) {
samer@0: 
samer@0: 				// look up exp( (2*PI*i/N)*p );
samer@0: 				wReal = wr[p]; 
samer@0: 				wImag = wi[p]; 
samer@0: 				p+=m;
samer@0: 
samer@0: 				for (i=k; i<N; i+=jump) {
samer@0: 
samer@0: 					int j=i+step;
samer@0: 
samer@0: 					// complex butterfly between i and j
samer@0: 					a1 = Real[i];
samer@0: 					a2 = Imag[i];
samer@0: 
samer@0: 
samer@0: 					b1 = wReal*Real[j] - wImag*Imag[j];
samer@0: 					b2 = wImag*Real[j] + wReal*Imag[j];
samer@0: 
samer@0: 					Real[i] = a1 + b1;
samer@0: 					Imag[i] = a2 + b2;
samer@0: 					Real[j] = a1 - b1;
samer@0: 					Imag[j] = a2 - b2;
samer@0: 				}
samer@0: 			}
samer@0: 			step=jump; jump<<=1; m>>=1;
samer@0: 		}
samer@0: 	}
samer@0: 
samer@0: 	// make everything conjugate - does inverse fft
samer@0: 	public void invert()	{
samer@0: 		for (int i=0; i<W.length; i++) W.imag[i] = -W.imag[i];
samer@0: 	}
samer@0: 
samer@0: 	// useful window functions
samer@0: 	public void setWindow(Function fn) {
samer@0: 		for (int i=0; i<N; i++) H[i] = (double)i/N;
samer@0: 		fn.apply(H); Mathx.mul(H,1/Math.sqrt(N));
samer@0: 	}
samer@0: }