/* this version puts things into files for future reference */
/* 06 apr 90 */





#include<\tc\include\stdlib.h>
#include<\tc\include\math.h>
#include<\tc\include\graphics.h>
#include<\tc\include\time.h>
#include<\tc\med\medlib\paramete.c>

/*	#include<\tc\med\medlib\inverse.c>
	#include<\tc\med\medlib\printa.c>*/

/*	#define L_Max 50   hey dumkopf in paramete.c L = 50 (max poss L) */

#define N 1
#define M 96

	int L;

	double data[N][M];
	double dataX[N][M+L_Max];

	double data_max = 0.0;



main()
{

#include <\tc\include\dos.h>

	int i,j;
	int gra = 1;
	db(0.0);


/*	get_data(M,&data[0][0]);*/

	get_data(M,&dataX[0][0]);

	if (gra == 1)
	{
		for (i = 0; i < N; i ++)
			for (j = 0; j < M; j ++)

				data[i][j] = dataX[i][j];
/*			printf("dataX[%i][%i] = %f\n",j,i,dataX[j][i]);
			getch();*/

	}


	if(gra == 1)
	{
	file_all(0,M,&data[0][0],"input");
	}

/*	if (gra == 1)
	{
		for (j = 0; j < M; j ++)
			dataX[0][j] = dataX[0][j] * .7;
	}
*/


	deconvolve1();


	if(gra == 1)
	{
	file_all(0,M,&dataX[0][0],"Output");
	}


	exit(0);
}

deconvolve1()
{

#define L_inc 1
#define L_begin 4
#define L_max 21		/*minus one */

	int ickx,icky;

	double C[N][L_Max];
	int i,j,k,s,u;
	double R[L_Max], W[L_Max][L_Max], F[L_Max], Y[N][M + L_Max];
	double expect_Y;
	double D_norm, D_norm_final;
	double medd_filter,medd_output;
	double Y_optimum[N][M + L_Max],E[L_Max];
	double optimum_D_norm, optimum_output;
	int optimum_L = L_begin;
	double f[N][M + L_Max];

	char filter_for_this_length[30];
	char output_for_this_length[30];

	filter_for_this_length[0] = "f";
	filter_for_this_length[1] = "i";
	filter_for_this_length[2] = "l";
	filter_for_this_length[3] = "_";

	output_for_this_length[0] = "o";
	output_for_this_length[1] = "u";
	output_for_this_length[2] = "t";
	output_for_this_length[3] = "_";



	for (i = 0; i < N; i ++)
	{
		for (j = 0; j < M + L_Max - 1; j ++)
		{
			Y_optimum[i][j] = 0;
		}

	}


	for (L = L_begin; L < L_max; L = L + L_inc)
	{

	D_norm_final = 0;

					/* calculate auto-corelation
					coefficients of input traces */


	for (k = 0; k < L; k ++)
	{

		for (i = 0; i < N; i ++)
		{
			C[i][k] = 0;
			for (j = 0; j < M; j ++)
			{
				if( j + k  >= 0 && j + k < M)
				{

/* printf("C: dataX[i][j]     = dataX[%i][%i] = %f\n",i,j,dataX[i][j]);*/
/* printf("C: dataX[i][j + k] = dataX[%i][%i] = %f\n",i,j+k,dataX[i][j+k]);*/

				    C[i][k] =
				    C[i][k] + dataX[i][j] * dataX[i][j+k];
				}
			}
			printf("C[%i][%i] = %f\n",i,k,C[i][k]);

		}

	}


	for (k = 0; k < L; k ++)
	{
		R[k] = 0;
		for (i = 0; i < N; i ++)
		{
			R[k] = R[k] + C[i] [k];

		 }
	printf("R[%i] = %f\n",k,R[k]);

	}







					/* inversion of autocorrelation
					matrix R[I-J] using
					Levinson Recursion:
					     -1
					W = R   (L x L) - matrix.	*/
	dbn(4.00);
	invert_matrix(L,L-1,&R[0],&R[1],&R[1],&W[0][0]);
	dbn(4.01);
/*	for (ickx = 0; ickx < L; ickx++)
	{
		for(icky = 0; icky < L; icky++)
		{

	printf("W[%i][%i] = %f\n",ickx,icky,W[ickx][icky]);
		}

	}

	getch();*/



	for (i = 0; i < N; i ++)
	{
		for (j = 0; j < M + L - 1; j ++)
		{

					/* calculation of [I] [J]
					filter, F	*/

			for (k = 0; k < L; k ++)
			{
				F[k] = 0;
				for (s = 0; s < L; s ++)
				{
					if ( j - s >= 0 && j - s < M)
					{
/* this assumes that j-s is terminative not cyclic*/


/*		printf("dataX[i][j-s] = dataX[%i][%i] = %#g\n",
		i,j-s,dataX[i][j-s]);*/

					F[k] = F[k] + dataX[i][j-s] * W[s][k];
					}

/*	printf("\t\ts,i,j-s,k,W[s][k],dataX[i][j-s],F[k]: %i %i %i %i %f, %f, %f\n",
	s,i,j-s,k,W[s][k],dataX[i][j-s],F[k]);*/

				}
/*				printf("\tFilter F[%i] = %#g\n",k,F[k]);*/
			}

/*getch();
for (icky = 0; icky < M; icky ++)
		printf("before s-k bussiness:	dataX[0][%i] = %#g\n",
		icky,dataX[0][icky]);
getch();*/
					/* calculation of the matrix Y
					corosponding to the output of the
					[I] [J] filter, F	*/
			for (k = 0; k < N; k ++)
			{
				for (s = 0; s < M+L-1; s ++)
				{
					Y[k] [s] = 0;
					for (u = 0; u < L; u ++)
					if (s - u + 1 >= 0 && s - u + 1 < M)
					{
/*printf("dataX[k][s-u+1] = dataX[%i][%i] = %f\n",k,s-u+1,dataX[k][s-u+1]);*/

						Y[k][s] =
					Y[k][s] + F[u] * dataX[k][s-u+1];
					}

/*				printf("\t\tY[%i][%i] = %#g\n",k,s,Y[k][s]);*/
				}
			}
						/* Calculation of ||Y|| */

/*
for (ickx = 0; ickx < M+L-1; ickx ++)
printf("\t\t\tY[%i][%i] = %#g\n",0,ickx,Y[0][ickx]);
*/

			expect_Y = 0;
			for (k = 0; k < N; k ++)
			{
				for (s = 0; s < M + L - 1; s ++)
					expect_Y =
					expect_Y + Y[k][s] * Y[k][s];

/*printf("\t\t\tk = %i, j = %i, M + L = %i, ||Y|| = %#g\n",
				k,j,M+L,expect_Y);*/
			 }
			expect_Y = sqrt(expect_Y);



						/* Calculation of the [I][J]
						coefficients of the [I][J]
						output normalized	*/

			if(expect_Y != 0.0)
				D_norm = (fabs(Y[i][j]))/(expect_Y);
			else
				D_norm = 0;

printf(".");
/*	printf("D_norm = %#g\n",D_norm);*/

			if (D_norm > D_norm_final)
			{
				D_norm_final = D_norm;
/*				medd_filter  = F[?];*/
				medd_output  = expect_Y;

			}


		}
	}
	printf("\n");

	E[L] = 0;
	for (i = 0; i < N; i ++)
	{
		for (j = 0; j < M + L - 1; j ++)
			E[L] = E[L] + Y[i][j] * Y[i][j];
	}

	printf("E[%i] = %#g\n",L,E[L]);


/*	for (ickx = 0; ickx < M+L-1; ickx ++)
printf("Y[%i][%i] = %f\n",0,ickx,Y[0][ickx]);*/

	for (j = 0; j < M; j ++)
	{
		f[0][j] = 0;
		if (j < L) f[0][j] = F[j];
	}


	file_all(0,M,&f,&filter_for_this_length);

	for (j = 0; j < M; j ++)
	{
		f[0][j] = 0;
		f[0][j] = Y[0][j];
	}
	file_all(0,M,&f,&output_for_this_length);


	printf("D_norm: %f\nmedd_output: %f\n",D_norm_final,medd_output);
	printf("number of filters: %i, error: %#g\n",L,E[L]);

/*	getch();*/

		 if (L == L_begin) goto label1;

		 if (E[L] <= E[optimum_L])
		 {
label1:		 optimum_L = L;
		 optimum_D_norm = D_norm_final;
		 optimum_output = medd_output;
			for (i = 0; i < N; i ++)
			{
				for (j = 0; j < M + L - 1; j ++)
				Y_optimum[i][j] = Y[i][j];
			}
		}
	}


	i = 0;
	for (L = L_begin; L < L_max; L = L + L_inc)
	{
		f[0][i] = log10(E[L]);
		i++;
	}
	file_all(0,i,&f,"E_vs_no");

	printf("FINAL RESULTS:\n");
	printf("D_norm: %f\nmedd_output: %f\n",optimum_D_norm,optimum_output);
	printf("Best error[# of filters] = E[%i] = %#g\n",
	optimum_L,E[optimum_L]);

/*	getch();*/

			for (i = 0; i < N; i ++)
			{
				for (j = 0; j < M + optimum_L - 1; j ++)
				{
					dataX[i][j] = 0;
					dataX[i][j] = Y_optimum[i][j];
				}
			}


}

xchange_data(double *in_data, double *out_data)
{
	int j;

	for (j = 0; j < M; j ++)

		*(out_data + j) = *(in_data + j);
}

get_data(int lngth, double *in_data)
{

	int wave_source = 0;			/* 0 => Cabrelli ex.1 */
						/* 1 => Gray ex. 1 */

	int i, j,k;
	double source_wavelet[M], spike_trace[M];

	double real_i, real_j, real_k, sum;

	double alpha,phi,time;
	double decay_rate;


	int reflect_1, reflect_2, reflect_3,reflect_4;

	double noise[M];

	double init_amp_0 =  1.0;
	double init_amp_1 =   .3880;
	double init_amp_2 =   .8274;
	double init_amp_3 =  -.6939;
	double init_amp_4 =   .9687;


	reflect_1 = lngth * .1059;
	reflect_2 = lngth * .1595;
	reflect_3 = lngth * .3050;
	reflect_4 = lngth * .6041;

	alpha = -18.778;
	alpha = -25;
	phi   = -.1960;
	decay_rate = 1.46;
	decay_rate = 5.5;


	for (j = 0; j < lngth; j ++)
	{
		real_j = j;
		time = real_j/lngth;
		printf("+");

		srand(j + 42);
		k = random(1000) - 500;
		real_k = k;
		noise[j] = (real_k/1000) * .02;

		source_wavelet[j] = 0;
		spike_trace[j] = 0;

		if(wave_source == 1) goto Gray1;

	/* create source wavelet (Cabrelli ex. 1) */

		source_wavelet[j] =
		(sin(alpha * time + phi) * exp(-decay_rate * time) *
		init_amp_0);


	/* create spike trace (as ex. 1) */

		if(j == reflect_1)
		spike_trace[j] = init_amp_1;

		if(j == reflect_2)
		spike_trace[j] = init_amp_2;

		if(j == reflect_3)
		spike_trace[j] = init_amp_3;

		if(j == reflect_4)
		spike_trace[j] = init_amp_4;

		if(wave_source == 0) goto Cab1;
		/* create wavelet as per Gray (ex. 1) */

Gray1:		source_wavelet[ 0] = 0.0;
		source_wavelet[ 1] = .033175;
		source_wavelet[ 2] = .076922;
		source_wavelet[ 3] = .165922;
		source_wavelet[ 4] = .52652;
		source_wavelet[ 5] = .80902;
		source_wavelet[ 6] = -.0039796;
		source_wavelet[ 7] = -.30637;
		source_wavelet[ 8] = -.30637;
		source_wavelet[ 9] = -.15650;
		source_wavelet[10] = .035809;
		source_wavelet[11] = .47878;
		source_wavelet[12] = -.17639;
		source_wavelet[13] = -.15385;
		source_wavelet[14] = -.013263;
		source_wavelet[15] = .31300;
		source_wavelet[16] = -.027852;
		source_wavelet[17] = -.14854;
		source_wavelet[18] = .014589;
		source_wavelet[19] = .13793;
		source_wavelet[20] = -.026526;
		source_wavelet[21] = -.12600;
		source_wavelet[22] = -.022547;
		source_wavelet[23] = .007957;


	/* create spike trace (as ex. 1) */

		spike_trace[ 1] = -1.00;
		spike_trace[ 8] =  -.582;
		spike_trace[22] =  -.900;
		spike_trace[28] =   .440;
		spike_trace[40] =  -.420;
		spike_trace[44] = -1.00;

	/* create convolution of source and spike */

Cab1:		sum = 0;
		for (i = 0; i <= j; i ++)
		{

			if (j-i >= 0 && j-i < lngth)
			sum =
			sum + source_wavelet[i] * (spike_trace[j-i]);
		}

		*(in_data + j) = sum + noise[j];
	}

/*	for (j = 0; j < lngth; j ++)
		*(in_data + j) = *(in_data + j) - sum;*/

	file_all(0,lngth,&source_wavelet,"Source");
	file_all(0,lngth,&spike_trace,"Reflect");
}

db(double idbp)
{
	printf("debug point %#g  \n",idbp);
/*	getch();*/
}
dbn(double i)
{
	printf("life is a bowl of fruits. \n");
	printf("no wait debug point %#g \n",i);
}	/* create source wavelet (Gray ex. 1) */