#include<stdio.h>
#include<math.h>
#define N 256

/* FFT関数 */
void fft(double xRe[],double xIm[],int p)
{
 	int i,ip,j,k,m,me,me1,n,nv2;
 	double uRe,uIm,vRe,vIm,wRe,wIm,tRe,tIm;
 	double sin(),cos();

	n=1<<p;
	nv2=n/2;

	/* ビットリバース */
 	j=0;
 	for(i=0;i<n-1;i++){
 		if(j>i){
 			tRe=xRe[j];
 			tIm=xIm[j];
 			xRe[j]=xRe[i];
 			xIm[j]=xIm[i];
 			xRe[i]=tRe;
 			xIm[i]=tIm;
 		}
 		k=nv2;
 		while(j>=k){
 			j-=k;
 			k/=2;
 		}
 		j+=k;
 	}
 
 	/* 主計算 */
	for(m=1;m<=p;m++){
 		me=1<<m;
 		me1=me/2;
 		uRe=1.0;
 		uIm=0.0;
 		wRe=cos(M_PI/me1);
 		wIm=-sin(M_PI/me1);
 		
 		for(j=0;j<me1;j++){
 			for(i=j;i<n;i+=me){
 				ip=i+me1;
 				/* バタフライ演算 */
 				tRe=xRe[ip]*uRe-xIm[ip]*uIm;
 				tIm=xRe[ip]*uIm+xIm[ip]*uRe;
 				xRe[ip]=xRe[i]-tRe;
 				xIm[ip]=xIm[i]-tIm;
 				xRe[i]+=tRe;
 				xIm[i]+=tIm;
 			}
 			vRe=uRe*wRe-uIm*wIm;
 			vIm=uRe*wIm+uIm*wRe;
 			uRe=vRe;
 			uIm=vIm;
		}
	}
}

/* main関数 */
int main(){
	int d,i,p;
	double xRe[N],xIm[N];
	FILE *fpr,*fpw;

	p=log10(N) / log10(2);

	/* ファイル入力 */
	fpr=fopen("./CYMBAL.txt","r");
	for(i=0;i<N;i++){
		fscanf(fpr,"%d",d);
		xRe[i]=(double)d;
		xIm[i]=0.0;
		//printf("%lf\n",xRe[i]);
		}
	fclose(fpr);

	/* FFT関数呼び出し */
	fft(xRe,xIm,p);
	for(i=0;i<N;i++){
		printf("%d,%lf,%lf\n",i,xRe[i],xIm[i]);
	}

	/* ファイル出力 */
	fpw=fopen("./test.txt","w");
	for(i=0;i<N;i++){
		fprintf(fpw,"%lf,%lf\n",xRe[i],xIm[i]);
	}
	fclose(fpw);
}