#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
#include <math.h>
#include <dir.h>
#define    SWAP(a,b)    tempr=(a); (a)=(b); (b)=tempr
#define    PI			3.1415926
#define	   FLT_ORD		5

struct wav_structure
{   char id[5], riff_type[5], fmt_id[5], data_id[5];
	unsigned long size, fmt_size, data_size, Fs, AvgBytesPerSecond;
	short format_tag, nr_canale, BlockAlign, BytesPerSamples;
}wav;
FILE *pf; int win_type; char wave_name[30];
unsigned long N, poz, poz_start, poz_end, step;
float energie_min, energie_max, F0a, F0c;

void read_wav_header(FILE *);
void write_wav_header(FILE *);
void eroare_wav(int);
long filesize(FILE *);
int str_compare(char *, char *);
void listing_fisiere(void);
void introducere_date_wav(void);
float* aloca_N_float(unsigned long);
void limite_energie(FILE *, float*);
void citeste_N_valori(FILE *, float *);
void window(float *, int);
void aplica_fereastra(float *, float *);
void F0_corel_cepstru_win(FILE *, float*, float*, float*);
void autocorelatie(float *, float *, unsigned int);
void high_pass_70Hz_Butter_filter(double *,double *);
void aplica_filtru(double *,double *,int);
void fourier(float *, unsigned long, int);
void corectie_ifft(float *);
void salvare_date_fft(char *, float *);

void main(void)
{	clrscr();
	printf("Program de calcul PITCH prin metoda cepstrala\n");
	listing_fisiere();
	introducere_date_wav();
	double A[FLT_ORD+1],B[FLT_ORD+1]; //coeficientii filtrului Butter de ordinul 5
	high_pass_70Hz_Butter_filter(B,A);
	aplica_filtru(B,A,FLT_ORD);

	float *data = aloca_N_float(2*N+1);
	float *win = aloca_N_float(N+1);
	float *corr = aloca_N_float(N);

	limite_energie(pf,win);
	window(win,win_type);
	FILE *pf_F0a;
	if((pf_F0a=fopen("F0autocor.txt","w+"))==NULL)
	{	printf("Eroare la creare fisier ce pastreaza valori F0 determinate prin metoda autocorelatiei"); getch(); free(win);	free(data); fclose(pf); exit(1);	}
	FILE *pf_F0c;
	if((pf_F0c=fopen("F0cepstru.txt","w+"))==NULL)
	{	printf("Eroare la creare fisier ce pastreaza valori F0 determinate prin metoda cepstrala"); getch(); free(win);	free(data); fclose(pf); exit(1);	}

	for(poz=poz_start; poz<=poz_end; poz+=step)
	{		F0_corel_cepstru_win(pf,data,corr,win);
			fprintf(pf_F0a,"%g\n",F0a);
			fprintf(pf_F0c,"%g\n",F0c);
	}

	free(corr);
	free(win);
	free(data);
	fclose(pf_F0a);	fclose(pf_F0c);
	fclose(pf);
}
//----------------------------------------------------------------------
void F0_corel_cepstru_win(FILE *pf, float* data, float* corr, float* win)
{   unsigned int i, f_min, f_max;
	float x, energie, max_corel, max_cepstru, energie_banda, energie_spectru;
	F0a=F0c=max_corel=max_cepstru=energie_banda=energie_spectru=energie=0.0;
	citeste_N_valori(pf, data);
	for(i=1; i<=2*N; i+=2) energie+=fabs(data[i]);//calculeaza energia pe fereastra
//conditii de existenta Pitch (F0)
	if (energie<energie_min*10 || energie<energie_max/20)
	{	F0a=F0c=0; return;	}

	autocorelatie(data+1, corr, N); // data[0] neutilizat
	for(i=wav.Fs/500; i<=wav.Fs/80; i++)
		if (corr[i]>max_corel)
		{	max_corel=corr[i]; F0a=i+1;	}

	aplica_fereastra(data, win);
	fourier(data,N,-1);				// FFT de fereastra de analiza ponderata
// calcul energie banda de frecvente [80 2500] pe jumatate din spectru
	f_min = ceil(80*N/wav.Fs)+1; f_max = floor(2500*N/wav.Fs)+1;
	for(i=1; i<=2*N; i+=2)
	{   data[i]=sqrt(data[i]*data[i]+data[i+1]*data[i+1]);	//abs spectru
		data[i+1]=0;
		if(i<=N) energie_spectru+=data[i];
	}
	for(i=f_min; i<=f_max; i++)
		energie_banda+=data[2*i-1];
//conditii de existenta Pitch (F0)
	if (energie_banda<energie_spectru*0.7)
	{	F0a=F0c=0; return;	}

	for(i=1; i<=2*N; i+=2)
		data[i]=log(data[i]+0.00000001);
	fourier(data,N,1); corectie_ifft(data);		// IFFT(log(abs(FFT)))
// jumatate din valorile cepstrului se regasesc oglindite in cealalta jumatate(ca la FFT)
	for(i=wav.Fs/500; i<=wav.Fs/80 &&i<=N/2; i++)	//parte reala cepstru
		if (data[2*i-1]>max_cepstru)
		{	max_cepstru=data[2*i-1]; F0c=i;	}
	F0a=(float)(wav.Fs/F0a); F0c=(float)(wav.Fs/F0c);
}
//----------------------------------------------------------------------
// Functia de autocorelatie lucreaza cu partea reala a unui vector complex
// data[0], data[2],..., data[2*(N-1)]
//----------------------------------------------------------------------
void autocorelatie(float data[], float corr[], unsigned int dim)
{	unsigned int i,j;
	for(i=0;i<dim;i++)
	{	corr[i]=0;
		for(j=i;j<dim;j++)
//			corr[i]+=data[j-i]*data[j];
			corr[i]+=data[2*(j-i)]*data[2*j];
	}
}
//----------------------------------------------------------------------
float* aloca_N_float(unsigned long N)
{	float *data;
	if ((data = (float *) malloc(N*sizeof(float))) == NULL)
	{	printf("\nEroare la alocare memorie pentru vectorul asociat ferestrei de analiza\n"); getch(); fclose(pf); exit(1);	}
	return data;
}
//----------------------------------------------------------------------
void limite_energie(FILE *pf, float* data)
{   long x;
	fseek(pf, 44, SEEK_SET);
	unsigned long masca_semn, masca=1, i, j, nwin;
	float suma=0;
	masca<<=wav.BytesPerSamples;
	masca_semn = masca>>1; masca--;
	for(j=0; j<N; j++)
	{	fread(&x,wav.BytesPerSamples/8,1, pf);
		if(x & masca_semn)
			x=(x & masca)-masca-1;	//numar negativ in C2-Complement fata de 2
		else
			x&=masca;				//numar pozitiv
		data[j] = (float)x/masca_semn;
		suma += fabs(data[j]);
	}
	energie_min=energie_max=suma;
	nwin=(poz_end-poz_start)+1; //nwin=(poz_end-poz_start)/step+1;
	j=j%N;
	for(i=1; i<=nwin; i++)
	{	fread(&x,wav.BytesPerSamples/8,1, pf);
		if(x & masca_semn)
			x=(x & masca)-masca-1;	//numar negativ in C2-Complement fata de 2
		else
			x&=masca;				//numar pozitiv
		suma -= fabs(data[j]);
		data[j] = (float)x/masca_semn;
		suma += fabs(data[j]);
		j++; j=j%N;
		if (energie_min>suma) energie_min=suma;
		if (energie_max<suma) energie_max=suma;
	}
}
//----------------------------------------------------------------------
void citeste_N_valori(FILE *pf, float* data)
{	long x;
//	fseek(pf,44,SEEK_SET); fseek(pf, (poz - N/2)*wav.BytesPerSamples/8, SEEK_CUR);
	fseek(pf, 44+(poz - N/2)*wav.BytesPerSamples/8, SEEK_SET);
	unsigned long masca_semn, masca=1; masca<<=wav.BytesPerSamples;
	masca_semn = masca>>1; masca--;
	for(int i=1; i<=2*N; i+=2)
	{	fread(&x,wav.BytesPerSamples/8,1, pf);
		if(x & masca_semn)
			x=(x & masca)-masca-1;	//numar negativ in C2-Complement fata de 2
		else
			x&=masca;				//numar pozitiv
		data[i] = (float)x/masca_semn;
		data[i+1] = 0;
	}
}
//----------------------------------------------------------------------
void window(float* win, int type)
{	for(int i=1; i<=N; i++)
	switch(type)
	{   case 0:win[i]=1; break;	//No window
		case 1:		//Hamming window
		win[i]=0.54-0.46*cos(2*PI*(i-1)/(N-1)); break;
		case 2:		//Hann window
		win[i]=0.5*(1-cos(2*PI*(i-1)/(N-1))); break;
		case 3:		//Blackman window
		win[i]=0.42-0.5*cos(2*PI*(i-1)/(N-1))+0.08*cos(4*PI*(i-1)/(N-1)); break;
	}
}
void aplica_fereastra(float* data,float* win)
{
	for(int i=1, j=1; j<=N; i+=2,j++) data[i]=data[i]*win[j];
}
//----------------------------------------------------------------------
//		Listing, citire, verficare corectitudine fisiere WAV
//----------------------------------------------------------------------
void listing_fisiere(void)
{	struct ffblk ffblk;
	int done;
	printf("Listing fisiere cu extensia *.wav din directorul curent\n");
	done = findfirst("*.wav",&ffblk,0);
	while (!done)
	{
		printf("  %s", ffblk.ff_name);
		done = findnext(&ffblk);
	}
}
//----------------------------------------------------------------------
void introducere_date_wav(void)
{	unsigned long nr_esantioane,i,j;
	int n;
	printf("\n\nNume fisier: "); gets(wave_name);
	if((pf=fopen(wave_name,"rb"))==NULL)
	{	printf("Eroare la deschidere fisier \"%s\"",wave_name); getch(); exit(1);	}
	read_wav_header(pf);
	printf("Frecventa esantionare: %ld, Nr.canale %d, Dimensiune esantion %d-biti", wav.Fs, wav.nr_canale, wav.BytesPerSamples);
	nr_esantioane = wav.data_size/(wav.BytesPerSamples/8);
	printf("\nNumar esantioane: %ld\n", nr_esantioane);
	while (N>nr_esantioane || N<4)
	{	printf("Introduceti dimensiunea ferestrei de analiza \n(putere a lui 2: 256, 512, 1024,...):");
		scanf("%ld",&N);
		for(i=1,j=N; j>0 && i!=N; i<<=1) j>>=1;
		N=i; //daca N nu este putere a lui 2, atunci se determina aceasta valoare
	}
	i=N/2; j=nr_esantioane-N/2;
/*	do{
		printf("Introduceti esantionul pe care va fi centrata prima fereastra (%ld - %ld): ", i, j);
		if((n=scanf("%ld",&poz_start))!=1) {printf("\nValoare implicita = %ld\n", i); poz_start=i;}
		printf("Introduceti esantionul pe care va fi centrata ultima fereastra (%ld - %ld): ", i, j);
		if((n=scanf("%ld",&poz_end))!=1) {printf("\nValoare implicita = %ld\n", j); poz_end=j;}
	}while(poz_start<i || poz_start>j || poz_end<i || poz_end>j || poz_end<poz_start);
*/ poz_end=j;poz_start=i;
	do{
		printf("Introduceti pasul de deplasare a fereastrei (1 - %ld): ", j-i);
		if((n=scanf("%ld",&step))!=1) {puts("\nValoare implicita = 1"); step=1;}
	}while( step<1 || step>j-i);
	do{
		printf("Tip fereastra (0-No window, 1-Hamming, 2-Hanning 3-Blackman): ");
		if((n=scanf("%d",&win_type))!=1) {puts("\nValoare implicita = 0"); win_type=0;}
	}while( win_type<0 || win_type>3);
}
//----------------------------------------------------------------------
void read_wav_header(FILE *pf)
{
	if(filesize(pf)<44) eroare_wav(0);
	fread(&wav.id, 4, 1,pf); wav.id[4]='\0';
	if(str_compare(wav.id, "RIFF")) eroare_wav(1);
	fread(&wav.size,sizeof(long),1, pf);
	if (wav.size != filesize(pf)) eroare_wav(2);
	fread(&wav.riff_type, 4, 1,pf); wav.riff_type[4]='\0';
	if(str_compare(wav.riff_type, "WAVE")) eroare_wav(3);
	fread(&wav.fmt_id, 4, 1,pf); wav.fmt_id[4]='\0';
	if(str_compare(wav.fmt_id, "fmt ")) eroare_wav(4);
	fread(&wav.fmt_size,sizeof(long),1, pf);
	if (wav.fmt_size != 0x10) eroare_wav(5);
	fread(&wav.format_tag,sizeof(short),1, pf);
	if (wav.format_tag != 1) eroare_wav(6);
	fread(&wav.nr_canale,sizeof(short),1, pf);
	if (wav.nr_canale != 1) eroare_wav(7);
	fread(&wav.Fs,sizeof(long),1, pf);
	fread(&wav.AvgBytesPerSecond,sizeof(long),1, pf);
	//AvgBytesPerSecond = Fs*BlockAlign
	fread(&wav.BlockAlign,sizeof(short),1, pf);
	fread(&wav.BytesPerSamples,sizeof(short),1, pf);
	//BlockAlign = BytesPerSamples/8*nr_canale
	fread(&wav.data_id, 4, 1,pf); wav.data_id[4]='\0';
	if(str_compare(wav.data_id, "data")) eroare_wav(8);
	fread(&wav.data_size,sizeof(long),1, pf);
	if (wav.data_size != filesize(pf)) eroare_wav(9);
	//..............urmeaza datele propriu-zise
}
void eroare_wav(int er)
{	switch(er)
	{   case 0:	puts("Dimensiune fisier prea mica (numai headerul are 44 octeti"); break;
		case 1:	puts("Semnatura wav incorecta (nu este \"RIFF\")"); break;
		case 2:	puts("Dimensiune wav din header este incorecta (nu este filesize-8)"); break;
		case 3:	puts("Identificator fisier wav incorect (nu este \"WAVE\")"); break;
		case 4:	puts("Identificator header fmt incorect (nu este \"fmt \")"); break;
		case 5:	puts("Dimensiune header fmt incorecta (nu este 16) "); break;
		case 6:	puts("Fisierul de sunet nu este de tip PCM/necompresat"); break;
		case 7:	puts("Fisierul de sunet nu este de pe un singur canal"); break;
		case 8:	puts("Nu este gasit identificator sectiune de date \"data\""); break;
		case 9:	puts("Dimensiunea sectiunii de date incorecta"); break;
	}
	getch(); exit(1);
}
void write_wav_header(FILE *pf)
{
	fwrite(&wav.id, 4, 1,pf);
	fwrite(&wav.size,sizeof(long),1, pf);
	fwrite(&wav.riff_type, 4, 1,pf);
	fwrite(&wav.fmt_id, 4, 1,pf);
	fwrite(&wav.fmt_size,sizeof(long),1, pf);
	fwrite(&wav.format_tag,sizeof(short),1, pf);
	fwrite(&wav.nr_canale,sizeof(short),1, pf);
	fwrite(&wav.Fs,sizeof(long),1, pf);
	fwrite(&wav.AvgBytesPerSecond,sizeof(long),1, pf);
	fwrite(&wav.BlockAlign,sizeof(short),1, pf);
	fwrite(&wav.BytesPerSamples,sizeof(short),1, pf);
	fwrite(&wav.data_id, 4, 1,pf);
	fwrite(&wav.data_size,sizeof(long),1, pf);
	//..............urmeaza datele propriu-zise
}

//----------------------------------------------------------------------
long filesize(FILE *stream)
{
   long curpos, length;

   curpos = ftell(stream);
   fseek(stream, 0L, SEEK_END);
   length = ftell(stream);
   fseek(stream, curpos, SEEK_SET);
   return length-curpos;
}
int str_compare(char *s1, char *s2)
{
	for(int i=0; s1[i]!=NULL || s2[i]!=NULL; i++)
	{	if(s1[i]<s2[i]) return -1;
		if(s1[i]>s2[i]) return 1;
	}
	return 0;
}
//----------------------------------------------------------------------
// Coeficientii sunt furnizati de functia <<butter>> din MATLAB
// ord=5; fe=22050; cutoff=2*70/fe; [B,A]=butter(N, cutoff, 'high');
//----------------------------------------------------------------------
void high_pass_70Hz_Butter_filter(double B[],double A[]) //Filtrul are ordinul 5
{
	switch(wav.Fs)
	{	case 8000:
			B[0]=0.914873096279762; B[1]=-4.57436548139881; B[2]=9.14873096279762; B[3]=-B[2]; 			   B[4]=-B[1]; 			  B[5]=-B[0];
			A[0]=1; 				A[1]=-4.82209440271224; A[2]=9.30410015818109; A[3]=-8.97888212541994; A[4]=4.3338696123426;  A[5]=-0.836992782296519; break;
		case 12000:
			B[0]=0.942416322024886; B[1]=-4.71208161012443; B[2]=9.42416322024886; B[3]=-B[2]; 			   B[4]=-B[1]; 			  B[5]=-B[0];
			A[0]=1; 				A[1]=-4.88139385024464; A[2]=9.53257868664335; A[3]=-9.3091187469501;  A[4]=4.54608249693936; A[5]=-0.888148524018915; break;
		case 16000:
			B[0]=0.956494956106162; B[1]=-4.78247478053081; B[2]=9.56494956106162; B[3]=-B[2]; 			   B[4]=-B[1]; 			  B[5]=-B[0];
			A[0]=1; 				A[1]=-4.91104475245164; A[2]=9.64812269167023; A[3]=-9.47799193668147; A[4]=4.65579661353731; A[5]=-0.914882601056528; break;
		case 24000:
			B[0]=0.970782788541803; B[1]=-4.85391394270902; B[2]=9.70782788541803; B[3]=-B[2]; 			   B[4]=-B[1]; 			  B[5]=-B[0];
			A[0]=1; 				A[1]=-4.9406961990485;  A[2]=9.76453946304156; A[3]=-9.64940919408792; A[4]=4.76798515463072; A[5]=-0.942419222529;    break;
		case 11025:
			B[0]=0.937485573285627; B[1]=-4.68742786642814; B[2]=9.37485573285627; B[3]=-B[2]; 			   B[4]=-B[1]; 			  B[5]=-B[0];
			A[0]=1; 				A[1]=-4.87090526008662; A[2]=9.49191453761557; A[3]=-9.24998466497119; A[4]=4.507854682348;   A[5]=-0.878879200118681; break;
		case 22050:
			B[0]=0.9682403079291;   B[1]=-4.8412015396455;  B[2]=9.682403079291;   B[3]=-B[2]; 			   B[4]=-B[1]; 			  B[5]=-B[0];
			A[0]=1; 				A[1]=-4.93545169389009; A[2]=9.74388511298395; A[3]=-9.618903937548;   A[4]=4.74795981541053; A[5]=-0.937489293898639; break;
		case 44100:
			B[0]=0.983992270398004; B[1]=-4.91996135199002; B[2]=9.83992270398004; B[3]=-B[2]; 			   B[4]=-B[1]; 			  B[5]=-B[0];
			A[0]=1; 				A[1]=-4.96772572987171; A[2]=9.87142312065382; A[3]=-9.80790980824479; A[4]=4.87245320576278; A[5]=-0.968240788203018; break;
		default:
			printf("\nFiltrul Butter nu suporta frecventa de esantionare FS:%d\n", wav.Fs); getch(); fclose(pf); exit(1);
	}
}

//----------------------------------------------------------------------
// H[z]=B[z]/A[z] => y[n]*A[0]=x[n]*B[0]+x[n-1]*B[1]+...-y[n-1]*A[1]-y[n-2]*A[2]-...
// Functia e optimizata in sensul ca pastreaza valorile intr-o lista circulara
// cu dimensiunea specificata de ordinul filtrului
//----------------------------------------------------------------------
void aplica_filtru(double B[],double A[], int ord)
{	char new_wave_name[30];
	long x;
	unsigned long masca_semn, masca=1, i;
	double in[FLT_ORD+1], out[FLT_ORD+1], suma;
	int first, j, k;

	FILE *pf_flt;
	for(i=0; wave_name[i-1]!='.' && wave_name[i]!=NULL;i++)
		new_wave_name[i]=wave_name[i];
	// se adauga extensia 'flt'
	new_wave_name[i++]='f';new_wave_name[i++]='l';new_wave_name[i++]='t';new_wave_name[i]='\0';
	if((pf_flt=fopen(new_wave_name,"w+b"))==NULL)
	{	printf("Eroare la creare fisier ce pastreaza fisierul wav filtrat"); getch(); exit(1);	}

	write_wav_header(pf_flt);

	fseek(pf, 44, SEEK_SET);	//pozitionare pe primul esantion din fisierul wav
	masca<<=wav.BytesPerSamples; masca_semn = masca>>1; masca--;
	first=FLT_ORD;
	for(i=0; i<=FLT_ORD; i++) {	in[i]=0; out[i]=0;	}
	for(i=0; i<wav.data_size/(wav.BytesPerSamples/8); i++)
	{	fread(&x,wav.BytesPerSamples/8,1, pf);
		if(x & masca_semn)
			x=(x & masca)-masca-1;	//numar negativ in C2-Complement fata de 2
		else
			x&=masca;				//numar pozitiv
		first++; first%=(FLT_ORD+1);
		in[first] = (double)x/masca_semn;
		suma=in[first]*B[0];
		for(k=1,j=(first-1)<0?FLT_ORD:first-1; k<=FLT_ORD; k++, j=(j-1)<0?FLT_ORD:j-1)
			suma+=in[j]*B[k]-out[j]*A[k];
		out[first]=suma/A[0];
		x=(long)(out[first]*masca_semn);
		if(x<0) x=(x+masca+1)&masca;
		fwrite(&x,wav.BytesPerSamples/8,1, pf_flt);
	}
//	fclose(pf_flt);
	fclose(pf); pf=pf_flt;
}
//----------------------------------------------------------------------
//			    Instrumente calculare, salvare FFT, IFFT
// Preluate din "Numerical Recipes in C. The Art of Scientific Computing"
// Second Edition, William H. Press, Saul A. Teukolsky, William T. Vetterling,
// Brian P. Flannery. Cambridge Uniersity Press, 1992. Chapte 12, pp. 496-536
//----------------------------------------------------------------------
void fourier(float data[], unsigned long nn, int isign)
{
	unsigned long n, mmax, m, j, istep, i;
	double wtemp, wr, wpr, wpi, wi, theta;
	float tempr, tempi;

	n = nn<<1;
	j = 1;
	for (i=1; i<n; i+=2)
	{   if (j>i)
		{    SWAP(data[j], data[i]);
			 SWAP(data[j+1], data[i+1]);
//			 printf("%ld <-> %ld\n", (i+1)/2, (j+1)/2);
		}
		m = n >> 1;
		while (m>=2 && j>m)
		{   j -= m;
			m >>=1;
		}
		j += m;
	}
	mmax=2;
	while (n>mmax)
	{   istep = mmax << 1;
		theta = isign*(6.28318530717959/mmax);
		wtemp = sin(0.5*theta);
		wpr = -2.0*wtemp*wtemp;
		wpi = sin(theta);
		wr = 1.0;
		wi = 0.0;
		for (m=1; m<mmax; m+=2)
		{   for (i=m; i<=n; i+=istep)
			{	j = i+mmax;
				tempr = wr*data[j]-wi*data[j+1];
				tempi = wr*data[j+1]+wi*data[j];
				data[j] = data[i]-tempr;
				data[j+1] = data[i+1]-tempi;
				data[i] +=tempr;
				data[i+1] +=tempi;
			}
			wr = (wtemp=wr)*wpr-wi*wpi+wr;
			wi = wi*wpr+wtemp*wpi+wi;
		}
		mmax=istep;
	}
}
void corectie_ifft(float data[])
{	for(int i=1; i<=2*N; i+=2)
		data[i]=data[i]/N;
}
//----------------------------------------------------------------------
void salvare_date_fft(char *wave_name, float data[])
{   int i, answer;
	FILE *pf_fft;
	for(i=0; wave_name[i]!='.' && wave_name[i]!='\0' ; i++);
	if (wave_name[i]=='\0') wave_name[i]='.';
	wave_name[++i]='f';wave_name[++i]='f';wave_name[++i]='t'; wave_name[++i]='\0';
	if((pf_fft=fopen(wave_name,"w+"))==NULL)
	{	printf("Eroare la creare fisier \"%s\"",wave_name); getch(); exit(1);	}
	printf("Salvare datele fft ca (0) numere complexe sau (1) modul: ");
	scanf("%d",&answer);
	for(i=1; i<=2*N; i+=2)
		if (answer)	fprintf(pf_fft,"%g\n", sqrt((double)(data[i]*data[i]+data[i+1]*data[i+1])));
		else		fprintf(pf_fft,"%g\t%g\n", data[i], data[i+1]);
	fclose(pf_fft);
}