diff src/samer/units/FFT.java @ 0:bf79fb79ee13

Initial Mercurial check in.

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