/* MAC 5749	Analise de Formas e Classificacao */
/* EP3		Convolucao utilizando FFT */
/* Entrega:	25/09/2003 */
/* Aluno:	Robson Ribeiro Correia */

/* Funcoes:
 *
 * void convlv (float data[], unsigned long n, float respns[], unsigned long n, int isign, float ans[])
 * void FFT(float data[], unsigned long nn, int isign)
 */


#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <math.h>
/* #include "nrutil.h" */

#define SWAP(a,b) tempr=(a); (a)=(b); (b)=tempr

static float sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)

#define NR_END 1
#define FREE_ARG char*

#define N 1024


main()
{
	/* Declaracao de variaveis */
	float data[2*N];	/* vetor de numeros complexos */



}


void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
	fprintf(stderr, "Numerical Recipes run-time error...\n");
	fprintf(stderr, "%s\n", error_text);
	fprintf(stderr, "...now exiting to system...\n");
	exit(1);
}


float *vector(long nl, long nh)
/* allocate a float vector with subscript range v[nl..nh] */
{
	float *v;

	v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
	if (!v) nrerror("allocation failure in vector()");
	return v-nl+NR_END;
}


void free_vector(float *v, long nl, long nh)
/* free a float vector allocated with vector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}


void convlv (float data[], unsigned long n, float respns[], unsigned long m, int isign, float ans[])

{
	void realfft(float data[], unsigned long n, int isign);
	void twofft(float data1[], float data2[], float fft1[], float fft2[], unsigned long n);
	unsigned long i, no2;
	float dum, mag2, *fft;

	fft= vector(1, n<<1);
	
	/* Put respns in array of length n */
	for(i=1;i<=(m-1)/2;i++)
		respns[n+1-1]=respns[m+1-i];

	/* Pad with zeros */
	for(i=(m+3)/2;i<=n-(m-1)/2;i++)
		respns[i]=0.0;
	
	/* FFT both at once */
	twofft(data, respns, fft, ans, n);
	no2=n>>1;
	for(i=2;i<=n+2;i+=2)
	{
		if(isign == 1)
		{
			/* Multiply FFTs to convolve */
			ans[i-1]=(fft[i-1]*(dum=ans[i-1])-fft[i]*ans[i])/no2;
			ans[i]=(fft[i]*dum+fft[i-1]*ans[i])/no2;
		}
		else if (isign == -1)
		{
			/*Divide FFTs to deconvolve */
			if((mag2=SQR(ans[i-1])+SQR(ans[i])) == 0.0)
				nrerror("Deconvolving at response zero in convlv");
			ans[i-1]=(fft[i-1]*(dum=ans[i-1])+fft[i]*ans[i])/mag2/no2;
			ans[i]=(fft[i]*dum-fft[i-1]*ans[i])/mag2/no2;
		}
		else
			nrerror("No meaning for isign in convlv");
	}
	
	/* Pack last element with first for realft */
	ans[2]=ans[n+1];
	
	/* Inverse transform back to time domain */
	realft(ans,n,-1);
	free_vector(fft,1,n<<1);
}


void FFT (float data[], unsigned long nn, int isign)

{
	unsigned long n, mmax, m, j, istep, i;
	double wtemp, wr, wpr, wpi, wi, theta;
	float tempr, tempi;
	
	/* sorts data in bit-reversed order */
	n = nn << 1;
	j = 1;
	for(i=1; i<n; i+=2)
	{
		if(j>i)
		{
			SWAP(data[j], data[i]);
			SWAP(data[j+1], data[i+1]);
		}
		m=nn;
		while(m>=2 && j>m)
		{
			j -= m;
			m >>= 1;
		}
		j += m;
	}

	/* calculates transforms of length 2, 4, 8, ..., N */
	mmax = 2;
	while(n>mmax)
	{
		istep = mmax << 1;				/* istep = mmax ^ 2 */
		theta = isign * (6.28318530717959/mmax);	/* theta = 2 PI i / N */
		wtemp = sin(0.5*theta);				
		wpr = -2.0*wtemp*wtemp;
		wpi = sin(theta);
		wr = 1.0;
		wi = 0.0;
		
		/* Two nested inner loops */
		for(m=1; m<mmax; m+=2)
		{
			for(i=m; i<=n; i+=istep)
			{
				j=i+mmax;
				
				/* This is the Danielson-Lanczos formula */
				tempr = wr*data[j] - wi*data[j+1];
				tempi = wr*data[j+1] + wi*data[j];
				data[j] = data[i] - tempr;
				data[j+1] = data [i+1] - tempi;
				data[i] += tempr;
				data[i+1] += tempi;
			}
			
			/* Trigonometric recurrence */
			wr = (wtemp = wr) * wpr - wi * wpi + wr;
			wi = wi * wpr + wtemp * wpi + wi;
		}
		mmax = istep;
	}
}