#include<stdlib.h>
	#include<math.h>
       	#include<graphics.h>
	#include<\bsm\medlib\paramete.c>
	#include<\bsm\medlib\inverse.c>
	#include<\bsm\medlib\printa.c>

	#define N 1
	#define M 96


	double data[N+1][M+1];
	double dataX[N+1][M+L+1+1];

	double data_max = 0.0;


main()
{
	#include <dos.h>
	int i,j;
	int gra = 1;

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



	get_data(M,&data[1][1]);

/*		for (j = 1; j < N+1; j ++)
			for( i = 1; i < M+1; i ++)
			printf("dataX[%i][%i] = %f\n",j,i,dataX[j][i]);

	getch();
*/	}


	if(gra == 1)
	{
	plot_all(M,&data[1][1],
	"Source wavelet convolved with reflection signal, x(i) = w(i) * q(i)");
	}

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

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

	deconvolve1();


	if(gra == 1)
	{
	plot_all(M,&dataX[1][1],
	"MEDD Output, Y(i) = f * x(i)");
	}


	exit(0);
}

plot_all(int lngth, double *data_pointer, char *id_string[30])

{
	graphics_on();
	plot_data(lngth, data_pointer);
	outtextxy(50,3,"Cabrelli example 1a, M = 96, L = 5");
	outtextxy(50,15,id_string);
	outtextxy(230,320,"Press any key to continue");
	getch();
	graphics_off();
}

graphics_on()
{
	int g_driver, g_mode, g_error;

	data_max = 0.0;
	detectgraph(&g_driver, &g_mode);
	initgraph(&g_driver, &g_mode, "");
}

graphics_off()
{
	closegraph();
}

plot_data(int lngth, double *in_data)
{
	int j;
	int int_x,int_data;
	int max_x, max_y;
	double real_max_x;

	max_x = getmaxx();
	max_y = getmaxy();

	real_max_x = max_x;

	if(data_max == 0.0)
		make_axis(max_x, max_y);

	scale_data(lngth, in_data, max_x, max_y);

/*           	*** pixel method of doiung things ***
	for (j = 0; j < lngth; j ++)
	{
		int_x = real_max_x/(lngth-1)*j;
		int_data = *(in_data + j);

		putpixel(int_x,int_data,5);
	}
		*** end of pixel method ***
*/

	setcolor(10);

	for (j = 1; j < lngth; j ++)

		line((real_max_x/(lngth-1)*(j-1)),*(in_data + j - 1),
		     (real_max_x/(lngth-1)*j),*(in_data + j));

	setcolor(15);

}

make_axis(int max_x,int max_y)
{
	setcolor(15);
	rectangle(0,0,max_x,max_y);
	line(0,max_y/2,max_x,max_y/2);

}

scale_data(int lngth, double *in_data,int max_x,int max_y)
{
	 int j;
	 double max_y_over_2 = max_y/2;
	 double real_data;

	 if(data_max == 0.0)
	 {

		for (j = 0; j < lngth; j ++)
		{
			if(*(in_data + j) < 0)
				real_data = - *(in_data + j);
			else
				real_data = *(in_data + j);

			data_max = max(data_max,real_data);

		 }
	 }

	 for (j = 0; j < lngth; j ++)
	 {
		real_data = max_y_over_2 - *(in_data + j)*max_y_over_2/data_max*.9;
		*(in_data + j) = real_data;
	 }

	 j = max_x;

}

deconvolve1()
{
	int ickx,icky;

	double C[N+1][L+1];
	int i,j,k,s,u;
	double R[L+1], W[L+1][L+1], F[L+1], Y[N+1][M+1 + L+1];
	double expect_Y;
	double D_norm, D_norm_final;
	double medd_filter[L+1],medd_output[N+1][M+1 + L+1];


	D_norm_final = 0;

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


	for (k = 0; k <= L-1; k ++)
	{
		for (i = 1; i <= N; i ++)
		{
			C[i][k] = 0;
			for (j = 1; j <= M; j ++)
			{
				if (j + k >= 1 && 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-1; k ++)
	{
		R[k] = 0;
		for (i = 1; 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.	*/


	invert_matrix(L-1,&R[0],&R[1],&R[1],&W[0]);


/*	for (ickx = 0; ickx < L+1; ickx++)
	{
		for(icky = 0; icky < L+1; icky++)
		{

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

	}

	getch();
*/


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

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

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


/*		printf("dataX[i][j-s] = dataX[%i][%i] = %f\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] = %f\n",k,F[k]);*/
			}

/*getch();
for (icky = 1; icky < M+1; icky ++)
		printf("before s-k bussiness:	dataX[0][%i] = %f\n",
		icky,dataX[0][icky]);
getch();*/

					/* calculation of the matrix Y
					corosponding to the output of the
					[I] [J] filter, F	*/
			for (k = 1; k <= N; k ++)
			{
				for (s = 1; s <= M+L-1; s ++)
				{
					Y[k] [s] = 0;
					for (u = 1; u <= L; u ++)
					if (s - u + 1 >= 1 &&
					    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] = %f\n",k,s,Y[k][s]);*/
				}
			}
						/* Calculation of ||Y|| */

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

			expect_Y = 0;
			for (k = 1; k <= N; k ++)
			{
				for (s = 1; 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|| = %f\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 = %f\n",D_norm);*/

			printf("d_norm = %f\n",D_norm);

			if (D_norm > D_norm_final)
			{
				D_norm_final = D_norm;
				for (u = 1; u <= N; u++)
				{
					for (s = 1; s <= M + L - 1; s ++)
					{
					medd_output[u][s] = Y[u][s];
					}
				}
				for (u = 1; u <= L; u++)
					medd_filter[u] = F[u];
			}
		}
	}
	printf("\n");
	for (i = 1; i <= N; i ++)
	{
		for (j = 1; j <= M + L - 1; j ++)
			dataX[i][j] = Y[i][j];
	}

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

	getch();
	for (j = 1; j <= M; j ++)
	{
		Y[1][j] = 0;
		if (j <= L) Y[1][j] = medd_filter[j];
	}



	plot_all(M,&Y[1][1],"MEDD filter values, f");

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



	plot_all(M,&Y[1][1],"MEDD filter last values, f");

	plot_all(M+L,&medd_output[1][1],"medd output");

}

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

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

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

get_data(int lngth, double *in_data)
{


	int i, j;
	double source_wavelet[M+1], spike_trace[M+1];

	double real_i, real_j, sum;

	double alpha,phi,time;
	double decay_rate;

	int reflect_1, reflect_2, reflect_3,reflect_4;

	double noise = 0;

	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;
	phi   = -.1960;
	decay_rate = 1.46;


	printf("reflect_2 = %i\n",reflect_2);

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


	/* 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) */

		spike_trace[j] = .0000000;

		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;

	/* create convolution of source and spike */

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

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

		*(in_data + j) = sum;

	}

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

	plot_all(lngth,&source_wavelet[1],"Source wavelet, w(i)");
	plot_all(lngth,&spike_trace[1],"Reflection signals, q(i)");
}