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



#define SAMPLE_RATE 16000
//#define STARTING_FREQ 300
//#define FREQ_STEP 3000
//#define ENDING_FREQ 8000
#define PI 3.141592653589793

//quick fix for scanf implimentation
int STARTING_FREQ;
int FREQ_STEP;
int ENDING_FREQ;

// precision transformer
int fixed_point_transform (int bits, int array_size, float* original , int* fixed_pt);



/*	wave generator and filter (do them both at the same time so that we dont have to have 
	a HUGE array laying around.  The wave can be modified by adjusting STARTING_FREQ, 
	FREQ_STEP and ENDING_FREQ.  This function writes out to file immediately with the format:

		<time (a double)> <original wave amp(float)>  <filtered wave amp (float)>

	filename written is waves.txt (I might make it a command line arg but not right now)
	This file will be fairly large, 2 seconds worth of waves will have 32000 lines.
*/
int wave_generator_and_filter (int taps_scale, int waves_scale, int array_size, int* fixed_pt);

int frequency_response(int taps_scale, int waves_scale, int array_size, int* fixed_pt);

/****** MAIN *****/
int main(int argc, char * argv[])
{
		
	int number_of_taps; // num of taps 
	int i;			// loop counter
	int wave_precision;  // # of bits to use for input wave
	int tap_precision; //# of bits to use for taps
	int wave_scale; // amount the input wave floats have been adjusted by (a power of 2)
	int tap_scale;  // amount that the floats have been adjusted by (a power of 2)
	
	int response;
	float largest_coef;  // largest coefficient (magnitude) found in file...could be used for good use of bits
	float* raw_coefficients; //array of floats read from the file
	int* rounded_coefficients;  // array of ints after binary bit precision is applied
	
	FILE* input_file;	

	response = '\0';
	largest_coef = 0.0;
	number_of_taps = 0; 
	wave_precision = 0;
	tap_precision = 0;

// for testing this is commented
	/*if (argc != 3){
		printf("Improper number of arguments\n");
		return -1;
	}
	
	
	precision = atoi (argv[2]);
	
	*/

	printf("\nEnter the number of bit precision for the input wave : ");
	scanf("%d",&wave_precision);
	printf("\nEnter the number of bit precision for the taps : ");
	scanf("%d:",&tap_precision);
	printf("\nInput Initial Frequency : ");
	scanf("%d",&STARTING_FREQ);
	printf("\nInput Frequency Step : ");
	scanf("%d",&FREQ_STEP);
	printf("\nInput Ending Frequency : ");
	scanf("%d",&ENDING_FREQ);

	input_file = fopen("taps.txt"/*filename*/, "r");

	if (input_file==NULL)
	{
		printf("File Not Found\n");
		return -1;
	}


	fscanf (input_file,"%i",&number_of_taps);

	raw_coefficients = malloc ((number_of_taps) * sizeof(float));

	for (i = 0; i < number_of_taps; i++)
	{
		if (!feof(input_file))
		{
			fscanf(input_file , "%f", &raw_coefficients[i]);
			if(raw_coefficients[i] > largest_coef || -1*raw_coefficients[i] > largest_coef || i == 0)
			{
				largest_coef = (raw_coefficients[i] > 0.0) ? raw_coefficients[i] : -1 * raw_coefficients[i];
			}
		}
		else
		{
			printf("improperly formatted file");
		}
	}

	fclose(input_file);


	rounded_coefficients = malloc (number_of_taps * sizeof(int));
	
	tap_scale = fixed_point_transform (tap_precision, number_of_taps, raw_coefficients , rounded_coefficients);

	wave_scale = (int)pow(2,wave_precision);

	wave_generator_and_filter (tap_scale, wave_scale, number_of_taps, rounded_coefficients);


	
	
	response = 0;
	printf("\n\n\nNOTE ** Calculating the Frequency Response of a filter takes some time \n");
	printf("(a few min) and can only changed if you have updated taps.txt.\n\n");
	printf("Enter 1 if you would like to calculate the Frequency Response of the filter : ");
	scanf("%d",&response);
	

	if (response == 1)
	{
		printf ("Calculating Frequency Response");
		frequency_response (tap_scale, wave_scale, number_of_taps, rounded_coefficients);
	}else
	{
		printf("Skipping Frequency Response");
	}
/*
	printf("Largest : %f  Scale : %d\n",largest_coef,scale);
	
	for (i = 0; i < number_of_taps; i++)
	{
		printf("%f (original) , %f (rounded)\n",raw_coefficients[i], (float)rounded_coefficients[i]/(float)scale);
	}
*/
	free (raw_coefficients);
	free (rounded_coefficients);
	return 0;

}





// this takes the float array and turns it into an int array that is a fixed point scalar 
// of the original.  It returns the scalar used to make the fixed pt.
int fixed_point_transform (int bits,int array_size, float* original , int* fixed_pt)
{

	int scalar; //2^bits .. nice to compute here if we want to ensure full bit usage we can modify our shift
	int i; //counter
	
	scalar= (int)pow(2, bits);
	
	for (i = 0; i < array_size; i++)
	{
		fixed_pt[i] = (int)(original[i] * (float)scalar);
	}

	return scalar;
}



int wave_generator_and_filter (int taps_scale, int waves_scale, int array_size, int* fixed_pt)
{	
	int index; //position in the circular storage bins
	int i; //for looping the bins
	int iterations;
	int number_of_freq;
	double d_tap_scale, d_wave_scale; //just limit the number of times I have to cast

	double current_time, time_step;
	double frequency, frequency_step;

	int * unfiltered_wave; //must hold wave for as many timesteps back as we have taps
	int filtered_wave;
	FILE* output_file;


		
//setup constants and resources
	time_step = (double)1.0 / (double)SAMPLE_RATE;
	number_of_freq = (ENDING_FREQ - STARTING_FREQ)/FREQ_STEP + 1;
	frequency_step = (double) FREQ_STEP;
	d_tap_scale = (double)taps_scale;
	d_wave_scale = (double) waves_scale;
	
	unfiltered_wave = malloc (array_size * sizeof(int));
	
	output_file = fopen("waves.txt","w");
	if (output_file == NULL)
	{
		printf("Cannot open output file\n");
		return -1;
	}

	

//initialization
	for(i = 0; i < array_size; i++)
		unfiltered_wave[i] = 0;
	index = 0;
	current_time = 0.0;

	
//running routine

	for(iterations = 0; iterations < 32000; iterations++)
	{
		//compute the unfiltered wave for time t
		frequency = (double)STARTING_FREQ;
		unfiltered_wave[index] = 0;
		for (i = 0; i < number_of_freq; i++)
		{
			unfiltered_wave[index] += (int)( (d_wave_scale*sin(frequency * 2.0 * PI * current_time))/
					((double)number_of_freq) );
			frequency += frequency_step;
		}
		
		
		// compute filtered wave for time t
		filtered_wave = 0;
		
		for (i = 0; i < array_size; i++)
		{
			filtered_wave += fixed_pt[i]*unfiltered_wave[(index + array_size - i)%array_size];
		}
		
		//send to file
		fprintf(output_file,"%lf\t%lf\t%lf\n",current_time, (double)(unfiltered_wave[index])/d_wave_scale, (double)(filtered_wave)/(d_wave_scale*d_tap_scale));

		//move to next time slice

		index = (index + 1)% array_size;
		current_time += time_step;
	}

// clean up
	fclose (output_file);
	free(unfiltered_wave);
	return 0;
}


//calculate the FR of the filter
int frequency_response(int taps_scale, int waves_scale, int array_size, int* fixed_pt)
{
	int index; //position in the circular storage bins
	int i; //for looping the bins
	int iterations;  //for looping over each freq
	int freq_max;

	int num_steps; //for log spacing

	double d_tap_scale; //just limit the number of times I have to cast
	double d_wave_scale;

	double current_time, time_step;
	int current_frequency, frequency_iteration;

	int * unfiltered_wave; //must hold wave for as many timesteps back as we have taps
	int filtered_wave;
	FILE* freq_file;


	//current_frequency [k] = (int)pow(8000,(k/num_steps) 	
//setup constants and resources
	time_step = (double)1.0 / (double)SAMPLE_RATE;
	
	current_frequency = 0;
	frequency_iteration = 0;
	num_steps = 8000/FREQ_STEP;

	d_tap_scale = (double)taps_scale;
	d_wave_scale = (double)waves_scale;
	
	unfiltered_wave = malloc (array_size * sizeof(int));
	
	freq_file = fopen("frequency.txt","w");
	if (freq_file == NULL)
	{
		printf("Cannot open output file\n");
		return -1;
	}

	while(current_frequency < 8000)
	{
		current_frequency = (int)pow( (double)8000, ( (double)frequency_iteration/(double)num_steps ) );
		
		//initialization for each frequency
		for(i = 0; i < array_size; i++)
			unfiltered_wave[i] = 0;
		index = 0;
		current_time = 0.0;
		freq_max = 0;
		
	//running routine

		for(iterations = 0; iterations < 4000; iterations++)
		{
			//compute the unfiltered wave for time t
			//frequency = (double)STARTING_FREQ;
				
			unfiltered_wave[index] =(int)(d_wave_scale*sin((double)current_frequency * 2.0 * PI * current_time));
			
			filtered_wave = 0;
			// compute filtered wave for time t
			for (i = 0; i < array_size; i++)
			{
				filtered_wave += fixed_pt[i]*unfiltered_wave[(index + array_size - i)%array_size];
			}
			
			if( abs(filtered_wave) > freq_max)
				freq_max = abs(filtered_wave);

			//move to next time slice

			index = (index + 1)% array_size;
			current_time += time_step;
		}
		//send max magnitude for current_freqency to file
		fprintf(freq_file,"%d\t%lf\n",current_frequency, (double)(freq_max)/(d_tap_scale*d_wave_scale));
		frequency_iteration++;
	}

// clean up
	fclose (freq_file);
	free(unfiltered_wave);
	return 0;
}