#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; char wave_name[30];
//unsigned long N, poz, poz_start, poz_end, step;

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);
void high_pass_70Hz_Butter_filter(double *,double *);
void aplica_filtru(double *,double *,int);

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);
	limite_energie(pf,win);
	window(win,win_type);
	FILE *pf_F0c;
	if((pf_F0c=fopen("F0cepstr.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)
			fprintf(pf_F0c,"%g\n",F0_cepstru_win(pf,data,win));

	free(win);
	free(data);
	fclose(pf_F0c);
*/
	fclose(pf);
}
//----------------------------------------------------------------------
//		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);
}
//----------------------------------------------------------------------
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;
}