/*  I scammed this from the CARL software  */

#include <stdio.h>
#include <math.h>

#define FORWARD 1
#define INVERSE 0
typedef struct { float re ; float im ; } complex ;

#define CABS(x) hypot( (x).re, (x).im )

complex cadd(), csub(), cmult(), smult(), cdiv(), conjg(), csqrt() ;

extern complex zero ;
extern complex one ;
extern float synt ;
float PI ;
float TWOPI ;


/*
 * 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.
 */
rfft( x, N, forward ) float x[] ; int N, 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 ;
    if ( first ) {
        PI = 4.*atan( 1. ) ;
        TWOPI = 8.*atan( 1. ) ;
        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.
 */
cfft( x, NC, forward ) float x[] ; int NC, forward ; {
 float wr, wi, wpr, wpi, theta, scale ;
 int mmax, ND, m, i, j, delta ;
    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
 */
bitreverse( x, N ) 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 ;
    }
}

/*************************************************************/
/*****  MY MAIN  *********************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

/********  NeXT Soundfile Header Struct   *******/
struct headerform {
    char pref[4];
    long hdr_length;
    long file_length;
    long mode;
    long samp_rate;
    long num_channels;
    char comment[16];
};

void main(int argc, char *argv[])   {
    float *result,temp;
    short *data;
    long i;
    int soundin;
    int length, power_of_two;
    FILE *input, *output;
    struct headerform hdr = {".snd",40,0,3,(long) 44100,1,"NBody File"};
    
    if (argc != 3)   {
        printf("useage: %s filein.fft<snd> fileout.fft<snd>\n", argv[0]);
	printf("        Converts .snd to .fft, and .fft to .snd\n");
	exit(0);
    }
    input = fopen(argv[1], "rb");
    output = fopen(argv[2], "wb");
    if (input && output)	{
	if (strchr(argv[1],'.')[1] == 's')	{
	    fread(&hdr.pref, 8, 1, input);
	    fread(&hdr.file_length, hdr.hdr_length-8, 1, input);
	    soundin = 1;
 	    length = hdr.file_length / 2;
 	}
	else {
	    length = 0;
	    while (fread(&temp,4,1,input)) length++;
	    fseek(input,0,0);
	    soundin = 0;
	    hdr.file_length = length * 2;
	}
	power_of_two = (log(length) / log(2));
	length = pow(2, power_of_two);
	printf("%i, %i, %i\n", hdr.file_length/2, power_of_two, length);
	
	data = (short *) malloc(length * 2);
	result = (float *) malloc(4 * length);
	if (soundin == 1)	{
	    fread(data, length, 2, input);
    	}
	else	{
	    fread(result,length,4,input);
	    fwrite(&hdr, hdr.hdr_length, 1, output);
	}

	if (soundin == 1)	{
	    for (i=0;i<length;i++)	
		result[i] = data[i];
	    rfft(result,length/2,FORWARD);
	    fwrite(result, length, 4, output);
	}
	else 	{
 	    rfft(result,length/2,INVERSE);
	    for (i=0;i<length;i++)	{
		if (fabs(result[i] > 32000))
			printf("%i %f Clipping here!!!\n",i,result[i]);
		data[i] = result[i];
	    }
	    fwrite(data, length, 2, output);
	}
	free(data);
	free(result);
    }
    else    {
	printf("Can't open input (or output) file!!\n");
    }
    
}