1 files changed, 1 insertions, 0 deletions
diff --git a/lib/fft.c b/lib/fft.c
new file mode 100644
index 0000000..d6bb1bc
--- /dev/null
+++ b/lib/fft.c
@@ -0,0 +1 @@
+#include "fftease.h"
+/* If forward is true, rfft replaces 2*N real data points in x with
+   N complex values representing the positive frequency half of their
+   Fourier spectrum, with x[1] replaced with the real part of the Nyquist
+   frequency value.  If forward is false, rfft expects x to contain a
+   positive frequency spectrum arranged as before, and replaces it with
+   2*N real values.  N MUST be a power of 2. */
+void rfft( float *x, int N, int forward )
+{
+  float 	c1,c2,
+  		h1r,h1i,
+		h2r,h2i,
+		wr,wi,
+		wpr,wpi,
+  		temp,
+		theta;
+  float 	xr,xi;
+  int 		i,
+		i1,i2,i3,i4,
+		N2p1;
+  static int 	first = 1;
+/*float PI, TWOPI;*/
+void cfft();
+    if ( first ) {
+	first = 0;
+    }
+    theta = PI/N;
+    wr = 1.;
+    wi = 0.;
+    c1 = 0.5;
+    if ( forward ) {
+	c2 = -0.5;
+	cfft( x, N, forward );
+	xr = x[0];
+	xi = x[1];
+    } else {
+	c2 = 0.5;
+	theta = -theta;
+	xr = x[1];
+	xi = 0.;
+	x[1] = 0.;
+    }
+    wpr = -2.*pow( sin( 0.5*theta ), 2. );
+    wpi = sin( theta );
+    N2p1 = (N<<1) + 1;
+    for ( i = 0; i <= N>>1; i++ ) {
+	i1 = i<<1;
+	i2 = i1 + 1;
+	i3 = N2p1 - i2;
+	i4 = i3 + 1;
+	if ( i == 0 ) {
+	    h1r =  c1*(x[i1] + xr );
+	    h1i =  c1*(x[i2] - xi );
+	    h2r = -c2*(x[i2] + xi );
+	    h2i =  c2*(x[i1] - xr );
+	    x[i1] =  h1r + wr*h2r - wi*h2i;
+	    x[i2] =  h1i + wr*h2i + wi*h2r;
+	    xr =  h1r - wr*h2r + wi*h2i;
+	    xi = -h1i + wr*h2i + wi*h2r;
+	} else {
+	    h1r =  c1*(x[i1] + x[i3] );
+	    h1i =  c1*(x[i2] - x[i4] );
+	    h2r = -c2*(x[i2] + x[i4] );
+	    h2i =  c2*(x[i1] - x[i3] );
+	    x[i1] =  h1r + wr*h2r - wi*h2i;
+	    x[i2] =  h1i + wr*h2i + wi*h2r;
+	    x[i3] =  h1r - wr*h2r + wi*h2i;
+	    x[i4] = -h1i + wr*h2i + wi*h2r;
+	}
+	wr = (temp = wr)*wpr - wi*wpi + wr;
+	wi = wi*wpr + temp*wpi + wi;
+    }
+    if ( forward )
+	x[1] = xr;
+    else
+	cfft( x, N, forward );
+}
+/* cfft replaces float array x containing NC complex values
+   (2*NC float values alternating real, imagininary, etc.)
+   by its Fourier transform if forward is true, or by its
+   inverse Fourier transform if forward is false, using a
+   recursive Fast Fourier transform method due to Danielson
+   and Lanczos.  NC MUST be a power of 2. */
+void cfft( float *x, int NC, int forward )
+{
+  float 	wr,wi,
+		wpr,wpi,
+		theta,
+		scale;
+  int 		mmax,
+		ND,
+		m,
+		i,j,
+		delta;
+void bitreverse();
+    ND = NC<<1;
+    bitreverse( x, ND );
+    for ( mmax = 2; mmax < ND; mmax = delta ) {
+	delta = mmax<<1;
+	theta = TWOPI/( forward? mmax : -mmax );
+	wpr = -2.*pow( sin( 0.5*theta ), 2. );
+	wpi = sin( theta );
+	wr = 1.;
+	wi = 0.;
+	for ( m = 0; m < mmax; m += 2 ) {
+	 register float rtemp, itemp;
+	    for ( i = m; i < ND; i += delta ) {
+		j = i + mmax;
+		rtemp = wr*x[j] - wi*x[j+1];
+		itemp = wr*x[j+1] + wi*x[j];
+		x[j] = x[i] - rtemp;
+		x[j+1] = x[i+1] - itemp;
+		x[i] += rtemp;
+		x[i+1] += itemp;
+	    }
+	    wr = (rtemp = wr)*wpr - wi*wpi + wr;
+	    wi = wi*wpr + rtemp*wpi + wi;
+	}
+    }
+/* scale output */
+    scale = forward ? 1./ND : 2.;
+    { register float *xi=x, *xe=x+ND;
+	while ( xi < xe )
+	    *xi++ *= scale;
+    }
+}
+/* bitreverse places float array x containing N/2 complex values
+   into bit-reversed order */
+void bitreverse( float *x, int N )
+{
+  float 	rtemp,itemp;
+  int 		i,j,
+		m;
+    for ( i = j = 0; i < N; i += 2, j += m ) {
+	if ( j > i ) {
+	    rtemp = x[j]; itemp = x[j+1]; /* complex exchange */
+	    x[j] = x[i]; x[j+1] = x[i+1];
+	    x[i] = rtemp; x[i+1] = itemp;
+	}
+	for ( m = N>>1; m >= 2 && j >= m; m >>= 1 )
+	    j -= m;
+    }
+}
+\ No newline at end of file