#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define PI 3.1415926535898
#define ROW		128
#define COL		128
#define LEVEL	256

int* sinw;
int* cosw;

extern void MYFFT(int* a,int* b,int* si,int* co);
//int n_fft(unsigned char *inimage, int *prdata,int *pidata,int *ap,
//				int row, int col);
int n_fft(unsigned char *inimage, double *prdata,double *pidata,double *ap,
				int row, int col);
//int fft(int n, int *ar, int *ai);
void swap(double *f1, double *f2);
double sqr(double ff);
//fast fourier convert

int fft(int n, double *ar, double *ai)
{
	int n2, a, c, d, f, g, h, j,count;
	double b, e, k, i, n1, co, si;
	int m, p, q, r, s;
	n1=log10(n)/log10(2);
	n2=-1;
	a=n;
	b=2*PI/n;
	for(c=1;c<=n1;c++)
	{
		d=a;
		a=a/2;
		e=0;
		for(f=0;f<a;f++)
		{
			co=cos(e);
			si=sin(e)*n2;
			e=e+b;
			for(g=d;g<=n;g+=d)
			{
				h=g-d+f;
				j=h+a;
				k=*(ar+h)-*(ar+j);
				i=*(ai+h)-*(ai+j);
				*(ar+h)=*(ar+h)+*(ar+j);                                    //butterfly process
				*(ai+h)=*(ai+h)+*(ai+j);
				*(ar+j)=(co*k+si*i);
				*(ai+j)=(co*i-si*k);
			}
		}
		b=2*b;
	}
	m=0;
	p=n/2;
	q=n-1;
	for(r=0;r<q;r++)
	{
		if((r-m)<0)
		{
			swap((double *)ar+r, (double *)ar+m);                            //sort the sequence
 			swap((double *)ai+r, (double *)ai+m);
		}
		s=p;
		while((s-m)<1)
		{
			m=m-s;
			s=s/2;
		}
		m=m+s;
	}
	for(s=0;s<n;s++)
	{ 
		*(ar+s)=*(ar+s)/n;                                                   //calculate the results
		*(ai+s)=*(ai+s)/n;
	}
	//for(count=0;count < n; count++)
	//{
	//	printf("%lf %lfi ",*(ar+count),*(ai+count));
	//}
	return 1;
}

/* f1, f2 swap */
void swap(double *f1, double *f2)
{
	double ttf;
	ttf= *f1;
	*f1= *f2;
	*f2=ttf;
	return;
}

double sq (double ff)
{
	return(ff*ff);
}
//whole fft
int n_fft(unsigned char *inimage, double *prdata,double *pidata,double *ap,
				int row, int col)
{	
	FILE* destfile;
	int ii,jj;
	int *ar;
	int *ai;
	double* ard;
	double* aid;
	//int gmax;
	double gmax;
	int j;
	char* destname = "d:\\fft.raw";
	
	destfile = fopen(destname,"wb+");
	if(destfile == (FILE*)0)
		return -1;
	
	ar=malloc(2*col*sizeof(int));
	ai=malloc(2*col*sizeof(int));
	ard=malloc(2*col*sizeof(double));
	aid=malloc(2*col*sizeof(double));
	//printf("%ld\n",&col);
	printf("\n\n%d\n\n%d\n\n",row,col);
	for(ii=0;ii<row;ii++)
	{
		for(jj=0;jj<col;jj++)
		{
			
			//*(prdata+ii*col+jj) = (int)*(inimage+ii*col+jj);
			*(prdata+ii*col+jj) = (double)*(inimage+ii*col+jj);
			*(pidata+ii*col+jj) = 0;
			//printf("%d ",*(inimage+ii*col+jj));
			
		}
		//printf("\n");
		
	}
	for(j=0; j < 128; j++)
	{
		
		printf(" %d %d %d\n",j,*(sinw + j),*(cosw + j));
		
	}
	
	printf("\n\n%d\n\n%d\n\n",sinw,cosw);
	for(ii=0;ii<row;ii++)
	{
		for(jj=0;jj<col;jj++)
		{
			*(ar+jj)=*(prdata+ii*col+jj)*100;                                
			*(ai+jj)=*(pidata+ii*col+jj)*100;
			//*(ard+jj) = (double)*(prdata+ii*col+jj); 
			//*(aid+jj) = (double)*(pidata+ii*col+jj);
			
			
		}
		printf("\n");
		printf("%d \n",col);
		//fft(col,ar,ai);
		//fft(col,ard,aid);
		MYFFT(ar,ai,sinw,cosw);
	//	for(jj=0;jj<col;jj++)
	//	{
	//		printf("%d %d ",*(ar+jj),*(ai+jj));
	//	}
	
		
		for(jj=0;jj<col;jj++)
		{
			*(prdata+ii*col+jj)=*(ar+jj)/100;
			*(pidata+ii*col+jj)=*(ai+jj)/100;                                 //after fft give the values back
			
			
		}
	}
	
	
	for(ii=0;ii<col;ii++)                                                 
	{
		for(jj=0;jj<row;jj++)
		{
			*(ar+jj)=*(prdata+jj*col+ii)*100;
			*(ai+jj)=*(pidata+jj*col+ii)*100;
		}
		//fft(row,ar,ai);
		MYFFT(ar,ai,sinw,cosw);                                               //fft to each column
		for(jj=0;jj<row;jj++)
		{
			*(prdata+jj*col+ii)=*(ar+jj)/100;
			*(pidata+jj*col+ii)=*(ai+jj)/100;                                 //after fft give the values back
		}
	}
	
	
	/*for(ii=0;ii<col;ii++)                                                 
	{
		                                              
		for(jj=0;jj<row;jj++)
		{
			printf("%d %d ",*(prdata+jj*col+ii),*(pidata+jj*col+ii));
	                             
		}
		printf("\n");
	}*/

	
	printf("the end\n");
	gmax=0;
	for(ii=0;ii<row;ii++)
		for(jj=0;jj<col;jj++)
		{
			*(ap+col*ii+jj)=sqrt(sq(*(prdata+ii*col+jj))+sq(*(pidata+ii*col+jj)));      //choose the max one
			if(gmax<*(ap+col*ii+jj))
				gmax=*(ap+col*ii+jj);
		}
		printf("\n\n%d\n\n",gmax);
	for(ii=0;ii<row;ii++)
	{
		for(jj=0;jj<col;jj++)
		{
			*(ap+ii*col+jj)=(unsigned char)(*(ap+col*ii+jj)*LEVEL/gmax);   //calculate the result
			//printf("%d ",*(ap+ii*col+jj));
			fputc((unsigned char)(*(ap+ii*col+jj)),destfile);
			
		}
		printf("\n");
	}
	free(ar);
	free(ai);
	return 1;
	
}

int main()
{
	//int *ap;
	double* ap;
	double *pr;
	double *pi;
	double *si;
	double *co;
	int *prd;
	int *pid;
	
	FILE* file;
	char* filename = "d:\\lab2_d.raw";
	
	int j = 0;
	int ch;
	int i = 0;
	int ci = 0;
	int ri = 0;
	unsigned char* matrix;
	matrix = (unsigned char*)malloc(128*128*sizeof(char));
	pr = (double*)malloc(128*128*sizeof(double));
	pi = (double*)malloc(128*128*sizeof(double));
	//ap = (int*)malloc(128*128*sizeof(int));
	ap = (double*)malloc(128*128*sizeof(double));
	prd = (int*)malloc(128*128*sizeof(int));
	pid = (int*)malloc(128*128*sizeof(int));
	for(ri = 0; ri < 128; ri++)
	{
		for(ci = 0; ci < 128; ci++)
		{
				*(matrix+ri*128 + ci) = 0;
		}
	
	}
	/*open the file*/
	printf("%s\n",filename);
	file = fopen(filename,"rb+");
	if(file == (FILE*)0)
		return -1;
	ri = 0;
	ci = 0;
	while((ch = fgetc(file))!=EOF)
	{
		
		if(i!= 0 && i % 128 == 0)
		{
		ri++;
		printf("\n");
		ci = 0;
		}
		*(matrix+ri*128 + ci) = ch;
	//	printf("%d ",ch);
		ci++;
		i++;
	}
	printf("\n\n\n\n\n\n%d\n\n\n\n\n",i);
	
	/*for(ri = 0; ri < 128; ri++)
	{
		for(ci = 0; ci < 128; ci++)
		{
				printf("%d ",*(matrix+ri*128 + ci));
		}
		printf("\n");
	}*/
	for(ri = 0; ri < 128; ri++)
	{
		for(ci = 0; ci < 128; ci++)
		{
				*(ap+ri*128 + ci) = 0;
				*(pr+ri*128 + ci) = 0;
				*(pi+ri*128 + ci) = 0;
				*(prd+ri*128 + ci) = 0;
				*(pid+ri*128 + ci) = 0;
				
		}
		
	}
	/*n_fft(matrix,pr,pi,ap,128,128);
	for(ri = 0; ri < 128; ri++)
	{
		for(ci = 0; ci < 128; ci++)
		{
				printf("%lf ",*(ap+ri*128 + ci));
		}
		printf("\n");
	}*/
	si = (double*)malloc(128*sizeof(double));
	co = (double*)malloc(128*sizeof(double));
	sinw = (int*)malloc(128*sizeof(int));
	cosw = (int*)malloc(128*sizeof(int));
	for(j=0; j < 128; j++)
	{
		*(si + j) = sin(j*2*PI/128);
		*(co + j) = cos(j*2*PI/128);
		*(sinw + j) =  *(si+j) * 128 * (-1);
		*(cosw + j) =  *(co+j) * 128;
		printf(" %d %d %d\n",j,*(sinw + j),*(cosw + j));
		
	}
	n_fft(matrix,pr,pi,ap,128,128);
	//n_fft(matrix,prd,pid,ap,128,128);
	
	
	
	
	
	return 0;
}