/*------------------------------------------------------------------------*/
/*									  */
/*			Program to compute DFT coefficients directly	  */
/*			3 other functions are used			  */
/*									  */
/*									  */
/*			Program design by: E C Ifeachor, July, 1992 	  */
/*									  */
/*			Program written  and modified 			  */
/*			by: Ronald Bazillion, 7/22/96	     		  */
/*									  */
/*------------------------------------------------------------------------*/

#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();
	dft();
	save_data();
	exit();
}
/*--------------------------------------------------------------*/
/* 	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);
}
/*------------------------------------------------------------------*/
/*								    */
/*	Function to compute the DFT of a discrete-time              */
/*	sequence directly					    */
/*								    */
/*								    */
/*      Program designed by: E C Ifeachor 31.10.91		    */
/*								    */
/*	Program written by: Ronald Bazillion 7/22/96		    */
/*								    */
/*------------------------------------------------------------------*/
void dft()
{
	extern	int inv;
	extern  long npt;
	long k,n;
	double WN, wk, c, s, XR[size], XI[size];
	extern complex x[size];

	WN=2*pi/npt;
	if(inv==1)
		WN=-WN;
	for (k=0;k<npt;++k)
	{
		XR[k]=0.0;
		XI[k]=0.0;
		wk=k*WN;
		for (n=0; n<npt;++n)
		{
		c=cos(n*wk);
		s=sin(n*wk);
		XR[k]=XR[k]+x[n+1].real*c+x[n+1].imag*s;
		XI[k]=XI[k]-x[n+1].real*s+x[n+1].imag*c;
	      	}
		if (inv==1)
		{	/* divide by N for IDFT */
			
			XR[k]=XR[k]/npt;
			XI[k]=XI[k]/npt;
		}
	}
	for (k=1;k<=npt; ++k)
	{
		/* store transformed data in x */
	
		x[k].real=XR[k-1];
		x[k].imag=XI[k-1];
	}
}