//this is a short C implementation to at Ex. 6.13 and 6.15
//since our input data are real, we will use realft() from
//	numerical recipe
//output is in freq instead of angular freq
#include <stdio.h>
#include <math.h>
#include "numcip.h"
#include "fft.h"

int maxN;//max length of data
double pi6;
double pi2;

main()
{
	int choice;
	int i,k;
	double *xdata;
	double funpara;//parameter for cos function
	int N;
	double delt,T;
	double freq,power;

	char filename[25];
	FILE *outfile;

	//system parameters
	maxN=64;
	pi2=8.0*atan(1.0);
	pi6=3.0*pi2;
	xdata=dvector(1,maxN);

	//user io
	printf("Choose the exercise:\n");
	printf("\t[0] Ex. 6.13\n");
	printf("\t[1] Ex. 6.15\n");
	scanf("%d",&choice);

	if(choice)
		{
		printf("\nEx 6.15: f(t)=cos 3 t\n");
		//period is pi/3.0;
		printf("Sampling rate is every second\n");
		delt=1.0;
		printf("Total length of sampled data:[8,16,32,64]\n");
		N=8;//starting sampling length
		funpara=3.0;
		
		}
	else
		{
		printf("\nEx 6.13: f(t)=cos 6pi t\n");
		//period is 1/3
		printf("The number of data point is fixed at N=32\n");
		N=32;
		printf("Sampling rate and sampling length will change accordingly\n");
		printf("\t[delt,totalT]\n");
		printf("\t[1,32],[0.5,16],[0.25,8]\n");
		delt=1.0;//starting sampling rate
		funpara=pi6;
		
		}

	//main part
	if(choice)
		{
		//Ex 6.15
		for(k=1;k<=4;k++)//outer loop over the diff cases
			{
			//generating data
			for(i=1;i<=N;i++)
				{
				xdata[i]=cos(i*delt*funpara);
				}
			realft(xdata,N,1);

			//output spectrum
			sprintf(filename,"ex15_%d.dat",k);
			outfile=fopen(filename,"w");
			
			fprintf(outfile,"%lg\t%lg\n",0.0,xdata[1]*xdata[1]/(N*N));//zero freq
			
			//other freqs
			freq=1.0/N;//freq increment 
			for(i=3;i<=N;i+=2)
				{
				//amplitude from real and imaginary parts
				power=xdata[i]*xdata[i]+xdata[i+1]*xdata[i+1];
				//normalize
				power/=N*N;
				fprintf(outfile,"%lg\t%lg\n",freq,power);
				freq+=1.0/N;
				}
			
			fprintf(outfile,"%lg\t%lg\n",0.5,xdata[2]*xdata[2]/(N*N));//critical freq
			fclose(outfile);

			N*=2;//increment sampling time
			}
		}
	else
		{
		//Ex 6.13
		for(k=1;k<=4;k++)//outer loop over diff cases
			{
			//generating data
			for(i=1;i<=N;i++)//number of data points is fixed at 32
				{
				xdata[i]=cos(i*delt*funpara);
				}
			
			realft(xdata,N,1);

			//output powere spec
			sprintf(filename,"ex13_%d.dat",k);
			outfile=fopen(filename,"w");
			
			fprintf(outfile,"%lg\t%lg\n",0.0,xdata[1]*xdata[1]/(N*N));//zero freq
			
			//other freqs
			T=N*delt;//total sampling time
			freq=1.0/T;//freq increment
			for(i=3;i<=N;i+=2)
				{
				//amplitude from real and imaginary parts
				power=xdata[i]*xdata[i]+xdata[i+1]*xdata[i+1];
				//normalize
				power/=N*N;
				fprintf(outfile,"%lg\t%lg\n",freq,power);
				freq+=1.0/T;//print out in units of pi
				}
			
			fprintf(outfile,"%lg\t%lg\n",0.5/delt,xdata[2]*xdata[2]/(N*N));//critical freq
			fclose(outfile);

			delt/=2.0;//halfing sampling rate
			}
		}

	return 0;
}