view src/samer/units/FFT.java @ 8:5e3cbbf173aa tip

Reorganise some more

/*
 *	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));
	}
}