#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

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; unsigned long N, poz; int win_type; char wave_name[30];

void read_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 citeste_N_valori(FILE *, float *);
void window(float *, int);
void aplica_fereastra(float *, float *);
void fourier(float *, unsigned long, int);
void corectie_ifft(float *);
void salvare_date_fft(char *, float *);

void main(void)
{	clrscr();
	printf("Program de calcul SPECTRU INSTANTANEU intr-o fereastra centrata pe un esantion\n\n");
	listing_fisiere();
	introducere_date_wav();
	float *data = aloca_N_float(2*N+1);
	float *win = aloca_N_float(N+1);
	window(win,win_type);
	citeste_N_valori(pf, data);
	aplica_fereastra(data, win);
	fourier(data,N,-1); //fourier(data,N,1); corectie_ifft(data);
	salvare_date_fft(wave_name, data);
	free(win); free(data);
	fclose(pf);
}
//----------------------------------------------------------------------
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 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];
}
//----------------------------------------------------------------------
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 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;
	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 fereastra (%ld - %ld): ", i, j);
		scanf("%ld",&poz);
	} while( poz<i || poz>j);
	do{
		printf("Tip fereastra (0-No window, 1-Hamming, 2-Hanning 3-Blackman): ");
		scanf("%d",&win_type);
	}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);
}
//----------------------------------------------------------------------
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;
}
//----------------------------------------------------------------------
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);
}