aboutsummaryrefslogtreecommitdiff
path: root/fft.c
diff options
context:
space:
mode:
Diffstat (limited to 'fft.c')
-rw-r--r--fft.c1
1 files changed, 1 insertions, 0 deletions
diff --git a/fft.c b/fft.c
new file mode 100644
index 0000000..d6bb1bc
--- /dev/null
+++ b/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