Make your own free website on Tripod.com

To see the Dsp1.h Header File
To see the dft.h Header File
To see the makefile
To see the input
To see the output
To go back to index
To go back to homepage

/*----------------------------------------------------------------------*/
/*									*/
/*	Program to compute DFT coefficients using DIT FFT		*/
/*	3 other functions are used					*/
/*									*/
/*	Program written by: EC Ifeachor, July, 1992			*/
/*									*/
/*									*/
/*----------------------------------------------------------------------*/

#include "dsp1.h"
#include "dft.h"

main()
{
	extern long npt;
	extern int inv;
	printf("select type of transform \n");
	printf("\n");
	printf("0 for forward DFT \n");
	printf("1  for inverse DFT \n");
	scanf("%d",&inv);
	read_data();
	fft();
	save_data();
	exit();
}

void fft()
{
	int sign;
	long m,irem,l,le,le1,k,ip,i,j;
	double ur,ui,wr,wi,tr,ti,temp;
	extern long npt;
	extern int inv;
	extern complex x[size];

	/* in-place bit reverse shuffling of data */

	j=1;
	for (i=1;i< npt; ++i)
	{
		if(i < j)
		 {
			 tr = x[j].real;
			 ti = x[j].imag;
			 x[j].real = x[i].real;
   			 x[j].imag = x[i].imag;
			 x[i].real = tr;
			 x[i].imag = ti;
			 k = npt/2;
			 while (k < j)
			 {	j-=k;
				k/=2;
			 }
		}
		 else
		{
			k = npt/2;
			 while (k < j)
		         {
				j-=k;
				k/=2;
			 }
		}
		j+=k;
	}
	/* calculate the number of stages: m=log2(npt), and whether FFT or IFFT */
	m=0;
	irem=npt;
          irem = npt;
	  while( irem > 1)
	  {
		irem = irem/2;
		m+=1;
	  }
	  if (inv == 1)
		sign =1 ;
	  else
	       sign=-1;
/* perform the FFT computation for each stage */

	for (l=1; l<=m; l++)
	{
		le = pow(2,l);
		le1 = le/2;
		ur = 1.0;
		ui = 0;
		wr = cos(pi/le1);
		wi = sign*sin(pi/le1);
		for (j=1; j<= le1; ++j)
		{
			i=j;
			 while(i <= npt)
			 {
				ip = i +le1;
				tr = x[ip].real*ur - x[ip].imag*ui;
				ti = x[ip].imag*ur + x[ip].real * ui;
				x[ip].real = x[i].real - tr;
				x[ip].imag = x[i].imag - ti;
				x[i].real = x[i].real + tr;
				x[i].imag = x[i].imag + ti;
				i = i + le;
			}
			temp = ur * wr - ui * wi;
			ui = ui * wr + ui * wi;
			ur = temp; 
		}
	}
	/* If inverse FFT is desired divide each coefficient by npt */
	if (inv == 1)
	{
		for (i = 1; i <= npt; ++i)
		{
			x[i].real = x[i].real/npt;
			x[i].imag = x[i].imag/npt;
		}
	}	
}
/*--------------------------------------------------------------*/
/* 	Function to read the data from a file			*/
/*								*/
/*  	Program designed by E C Ifeachor July, 1992		*/
/*								*/
/*	Written by Ronald Bazillion 7/22/96			*/
/*								*/
/*--------------------------------------------------------------*/

void read_data()
{
	extern long npt;
	int n;
	extern complex x[size];

	for (n=0;n< size; ++n)
	{
		x[n].real=0;
		x[n].imag=0;
	}
	if ((in=fopen("coeff.dat","r"))==NULL)
	{
		printf("cannot open file coeff.dat\n");
		exit(1);
	}
	fscanf(in,"%ld",&npt);
	for(n=1; n<=npt; n++)
	{
		fscanf(in,"%lf %lf",&x[n].real,&x[n].imag);
	}
	fclose(in);
}
/*--------------------------------------------------------------*/
/* 	Function to save the data to a file 			*/
/*								*/
/*  	Program designed by E C Ifeachor July, 1992		*/
/*								*/
/*	Written by Ronald Bazillion 7/22/96			*/
/*								*/
/*--------------------------------------------------------------*/

void save_data()
{
	long k;
	int k1;
	extern  long npt;
	extern complex x[size];
	if((out=fopen("dftout.dat","w"))==NULL)
	{
		printf("Cannot open file dftout.dat\n");
		exit(1);
	}
	fprintf(out,"k \t\t XR(k) \t\t XI(k) \n");
	fprintf(out,"\n");
	for (k=1; k<=npt; ++k)
	{
		k1=k-1;
		fprintf(out, "%d \t %f \t %f \n", k1, x[k].real, x[k].imag);
	}
	fclose(out);
}