#include<\c\include\stdlib.h>
	#include<\c\include\math.h>
	#include<\c\include\graphics.h>
	#include<\c\eellib\pllib.h>


	#define N 1
	#define M 128
	#define nrows N
	#define row_lngth M*2



typedef struct
{
	double Re;
	double Im;

} complex;


	complex data[N][M];
	complex output[N][M];

	int plbuf[3072];
	int n2p;

double exponent(double a, double exponnt)
{
	double z;

	if (a != 0.0)
		z = exp(exponnt * log(a));
	else
		z = 0;

	printf("expontent: a = %e, e = %e, z = %e\n",a,exponnt,z);

	return(z);
}

double sign(double a)
{
	double c;

	if (a >= 0.0)
		c = 1;
	else
		c = -1;

/*	printf("sign: a = %f, c = %f\n",a,c);*/

	return(c);
}

main()
{
	int i,j;


	init_pl();
	getch();
/*	nrows = 1;
*/	n2p = log(row_lngth)/log(2);



	get_data(M,&data[0][0]);
	plot_complex_data(1,M,&data,"The synthetic wavelet");
	fft(1,row_lngth,nrows,n2p,&data,&data);  /* First do forward real fft's on rows. */
/*	print_complex(M,data);*/
	plot_complex_data(2,M/2,&data,
	"FFT of Source wavelet convolved with reflection signal, x(i) = w(i) * q(i)");
	fft(-1,row_lngth,nrows,n2p,&data,&data);      /* First do forward real fft's on rows. */
	plot_complex_data(1,M,&data,
	"IFFT of Source wavelet convolved with reflection signal, x(i) = w(i) * q(i)");

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

	plot_complex_data(1,M,&data,
	"Source wavelet convolved with reflection signal, x(i) = w(i) * q(i)");
skip:

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


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




	deconvolve(&data[0][0],&output[0][0]);


	plot_complex_data(1,M,&data,
	"MEDD Output, Y(i) = f * x(i)");


    return;
}

deconvolve(complex *y, complex *x)
{

/*	Gray's MED algorithm:

	Wiggin's MED algorithm with an exponentail transform
	in Fourier (frequency) space.

	conventions:

		all varibles that are lower case => time domain
		while their frequency domain counterparts are represented
		by the upper case letters.

	inputs:
		input complex matrix (address given by *y)

	outputs:
		MED processes output complex matrix (address given by x)

*/
	#define convergence_criteria .001

	int ir,it,iw;                   /* not used */

	int t,w,r,l;			/* element counters for main
					   matrices, t = tempus, w =
					   frequence, r = channel #,
					   l for padding elements */

	int t0;				/* tempus offset for filter
					   value is totally arbitary */

	int os;				/* temporary offset for pointer
					   arithematic		*/

					/* the filter vectors and their
					   respective pointers */
	complex f[N][M];
	complex F[N][M];

					/* y (input matrix) in frequency
					   domain and its pointer */
	complex Y[N][M];

					/* x (output matrix) in frequency
					   domain and its pointer */
	complex X[N][M];


	complex R[N][M];	/* frequency domain of autocorrelations */
	complex G[N][M];	/* "                 " crosscorrelations cubed */


	double a = 4.2;          /* Gray's alpha */
	complex a_minus_1;


	complex u1[N];
	complex u2[N];

	double   sum_of_xs;
	complex sum_of_Rs;
	complex sum_of_Fs;

	complex k = {0.000001,0.0};	/* Gray's lambda constant */

	double M_over_a = M/a;
	complex one_over_M = {1.0/M,0.0};

	complex Tmp[5], *tmp = &Tmp[0];

	complex x_tmp;

/*	complex norm_numberator;
	complex norm_denominator;
*/
	double norm_numberator;
	double norm_denominator;

	double V_norm;
	double last_V_norm = 0.0;

	int iteration = 0;

	complex One = {1.0, 0.0};
	complex Two = {2.0, 0.0};
	double minus_one = -1.0;


	t0 = 0;    		/* set up initial filter f0 as a simple
				   delta function, location of the spike is
				   arbatarary */


printf("here begins deconvolve...\n");
getch();

	a_minus_1.Re  = a-1.0, a_minus_1.Im  = 0.0;
/*	one_over_M.Re = 1.0/M, one_over_M.Im = 0.0;*/

			/* Initialize filter estimate f^0 */


	for (r = 0; r < N; r ++)
	{
		for (t = 0; t < M; t ++)
		{
			f[r][t].Re = 0.0, f[r][t].Im = 0.0;
		}

	f[r][t0].Re = 1.0, f[r][t0].Im = 0.0;

	}

/*printf("t0 = %i\n",t0);*/

	plot_complex_data(1,M,&f[0][0],
	"first approx. to filter f");


	fft(1,row_lngth,nrows,n2p,&f[0][0],&F[0][0]);

	plot_complex_data(2,M/2,&F,
	"first approx. to filter f, the FFT, F");

				/* input each trace (Gray calls them channels)
				   and process it with the MEDD algortihm */

iterate:

	for (r = 0; r < N; r++)
	{

				/* zero out R and W at the beginning of
				   each iteration */

		for (w = 0; w < M; w ++)
		{
			R[r][w].Re = 0, R[r][w].Im = 0;
			G[r][w].Re = 0, G[r][w].Im = 0;
		}

				/* take FFT of input matrix */

/* note: this is wrong, must pad y then take the FFT, not the way
	 shown here.  I don't knwo how to fix this yet */


	plot_complex_data(1,M,y,
	"the input matrix, y");

	fft(1,row_lngth,nrows,n2p,y,&Y[0][0]);

	plot_complex_data(2,M/2,&Y,
	"the FFT of the input matrix, Y");

				/* provide zero padding on Y */

/*		for (l = 0; l < M; l ++)
		{

			Y[r][M + l].Re = 0, Y[r][M + l].Im = 0;
		}
*/
				/* calculate the output matrix X
				   using the approximation F */

		for (w = 0; w < M; w ++)
		{
			complex_equate(&X[r][w],
				complex_times(tmp,&F[r][w], &Y[r][w]));

/*printf("Y[%i][%i] = (%f,%f), X = (%f,%f), F = (%f,%f)\n",
r,w,Y[r][w].Re,Y[r][w].Im,X[r][w].Re,X[r][w].Im,F[r][w].Re,F[r][w].Im);
getch();*/
		}

/*printf("X=F*Y done\n");
  */
				/* take inverse fft of X, use x to
				   get the scaling factors u1 and u2 */


	plot_complex_data(2,M/2,&X,
	"the FFT of the output matrix, X");

		fft(-1,row_lngth,nrows,n2p,&X[0][0],x);

	plot_complex_data(1,M,x,
	"the output matrix, x");

		u1[r].Re = 0.0, u1[r].Im = 0;
		u2[r].Re = 0.0, u2[r].Im = 0;

		for (t = 0; t < M; t ++)
		{
		     os = r*M + t;

/*printf("x[r][t] = (%e,%e) os = %i\n",((x+os)->Re),(x+os)->Im,os);*/



		     complex_equate(&u1[r],
			   complex_add(tmp+2, &u1[r],
				complex_exponent(tmp+1,
					complex_abs(tmp,(x + os)),&a)));


		     complex_equate(&u2[r],
			   complex_add(tmp+1,&u2[r],
				complex_exponent(tmp,(x + os),&Two)));

/*printf("u1 = (%e,%e)  u2 = (%e,%e)\n",u1[r].Re,u1[r].Im,u2[r].Re,u2[r].Im);
getch();*/

		}

					/*      u1 = 1/((1/M) * u1)     */
		complex_equate(&u1[r],
		  complex_exponent(tmp+2,
		    complex_times(tmp+1,&u1[r],&one_over_M),&minus_one));

					/*	u2 = 1/((1/M) * u2)      */

		complex_equate(&u2[r],
		  complex_exponent(tmp+2,
		    complex_times(tmp+1,&u2[r],&one_over_M),&minus_one));

						     /* x = exp(abs(x**(a-1))) *
							    sgn(x)  */
printf("u1 = (%e,%e)  u2 = (%e,%e)\n",u1[r].Re,u1[r].Im,u2[r].Re,u2[r].Im);
getch();

		for (t = 0; t < M; t ++)
		{
			os = r*M + t;

printf("a-1 = %4f, x[r][t] = (%e,%e) os = %i ",a_minus_1.Re,(x+os)->Re,(x+os)->Im,os);
			x_tmp.Re = sign((x + os) -> Re), x_tmp.Im = 0.0;

/*printf("x-tmp = (%e,%e)\n",x_tmp.Re,x_tmp.Im);*/

			complex_equate(( x + os),
			  complex_times(tmp+2,&x_tmp,
			    complex_exponent(tmp+1,
			      complex_abs(tmp,(x + os)),&a_minus_1)));

printf("new estimate of x[r][t] = (%e,%e) (os = %i)\n",
	(x + os)->Re,(x +os)->Im,os);


/*	getch();*/

		}
/*printf("now have new estimation of x\n");
getch();*/

		fft(1,row_lngth,nrows,n2p,x,&X[0][0]);


	plot_complex_data(2,M/2,&X,
	"the FFt of the output matrix, X");



	plot_complex_data(1,M,x,
	"the new estimation of output matrix, x");



/* it works to here 12:00 7-mar-89 */


		for (w = 0; w < M; w ++)
		{
							/*		*
							   R = R + (u1*Y  Y)
							*/


printf("getting R: Y[%i][%i] = (%e,%e)\n",r,w,Y[r][w].Re,Y[r][w].Im);


		  complex_equate(&R[r][w],
		    complex_add(tmp+3,&R[r][w],
		      complex_times(tmp+2,&u1[r],
			complex_times(tmp+1,&Y[r][w],
			  complex_conjagate(tmp,&Y[r][w])))));

/*printf("u1 = (%e,%e)\n",u1[r].Re,u1[r].Im);*/
/*printf("R[%i] = (%e,%e)\n",w,R[r][w].Re,R[r][w].Im);*/

/*printf("getting G \nY[%i][%i] = (%e,%e) X = (%e,%e)\n",
	r,w,Y[r][w].Re,Y[r][w].Im,X[r][w].Re,X[r][w].Im);*/

		  complex_equate(&G[r][w],
		    complex_add(tmp+3,&G[r][w],
		      complex_times(tmp+2,&u2[r],
			complex_times(tmp+1,&X[r][w],
			  complex_conjagate(tmp,&Y[r][w])))));

/*printf("u2 = (%e,%e)\n",u2[r].Re,u2[r].Im);                */
/*printf("G[%i] = (%e,%e)\n",w,G[r][w].Re,G[r][w].Im);*/
/*getch();*/
/*if(X[r][w].Re != 0.0) getch();*/
		}


printf("got R & G");
getch();

G[r][0].Re = 0, G[r][0].Im = 0;


				/* process next trace */
	}

/* works to here 11:34 08-mar-89 */


				/* calculate new filter vector */

	for (r = 0; r < N; r ++)
	{
		sum_of_Rs.Re = 0, sum_of_Rs.Im = 0;


					/* get sum over all w of R[r] */

		for (w = 0; w < M; w ++)
		{
			complex_equate(&sum_of_Rs,
				complex_add(tmp,&sum_of_Rs,&R[r][w]));

/*printf("R[w] = (%e,%e), sum of R from 0 to %i = (%e,%e)\n",
	R[r][w].Re,R[r][w].Im,w,sum_of_Rs.Re,sum_of_Rs.Im);*/
/*getch();*/

		}
printf("equated sum of R: (%e,%e)\n",sum_of_Rs.Re,sum_of_Rs.Im);
getch();

					/* get unnormalized F[r][w] */

		for (w = 0; w < M; w ++)
		{
printf("R[%i] = (%e,%e)  ",w,R[r][w].Re,R[r][w].Im);
printf("G[%i] = (%e,%e)  \n",w,G[r][w].Re,G[r][w].Im);



			complex_equate(&F[r][w],
			  complex_divide(tmp+2,&G[r][w],
			    complex_add(tmp+1,&R[r][w],
			      complex_times(tmp,&k,&sum_of_Rs))));

/*printf("k = (%e,%e), F[%i] = (%e,%e)\n",k.Re,k.Im,w,F[r][w].Re,F[r][w].Im);
getch();*/
		}

					/* get the sum of F[r]
					   over all w */



		sum_of_Fs.Re = 0, sum_of_Fs.Im = 0;

		for (w = 0; w < M; w ++)
		{
			complex_equate(&sum_of_Fs,
			  complex_add(tmp+2,&sum_of_Fs,
			    complex_times(tmp+1,&F[r][w],
			      complex_conjagate(tmp,&F[r][w]))));

/*printf("F[%i] = (%e,%e), sum = (%e,%e)\n",
	w,F[r][w].Re,F[r][w].Im,sum_of_Fs.Re,sum_of_Fs.Im);
getch();*/
		}
					/* normalize F[r][w] */

printf("Sum( F(n) ) = (%e,%e)\n\n",sum_of_Fs.Re,sum_of_Fs.Im);

		for (w = 0; w < M; w ++)
		{

			complex_equate(&F[r][w],
				complex_divide(tmp,&F[r][w],&sum_of_Fs));

/*printf("F[%i] = (%e,%e)\n",w,F[r][w].Re,F[r][w].Im);
getch();
*/
		}
	}
getch();

	plot_complex_data(2,M/2,&X,
	"the FFt of the filter matrix, F");

	fft(-1,row_lngth,nrows,n2p,&F[0][0],&f);

/*for (w = 0; w < M; w ++)
	printf("x = (%e,%e), f = (%e,%e), y = (%e,%e)\n",
		(x + w) -> Re,(x + w) -> Im,
		f[0][w].Re,f[0][w].Im,
		(y + w) -> Re,(y + w) -> Im);
getch();*/

	plot_complex_data(1,M,&f,
	"the filter matrix, f");

				/* calcualte V norm and check for
				   convergence */
	V_norm = 1;

	for (r = 0; r < N; r ++)
	{                                             /* numberator is 1/u1
							 raised to the M/a */

		norm_numberator = pow(u1[r].Re, -(M_over_a));

printf("print |X|^a = %e, M/a = %e, norm top: %e\n",
	1/u1[r].Re,M_over_a,(exp(( -128.0/2.2 ) * log(u1[r].Re))),norm_numberator);
getch();

		norm_denominator = 0;

		for (t = 0; t < M; t ++)
		{
			complex_equate(&x_tmp,
			  complex_abs(tmp,&X[r][t]));

/*printf("X = (%e,%e), |X| = %e ",X[r][t].Re,X[r][t].Im,
			x_tmp.Re);                      */

			norm_denominator =
			  norm_denominator + pow(x_tmp.Re,2.0);

/*printf("sum of |X|^2 = %e\n",norm_denominator);
getch();*/

		}

		norm_denominator =
			pow((one_over_M.Re * norm_denominator),M/2.0);

printf("norm numberator = %e, norm denominator = %e\n",
	norm_numberator,norm_denominator);
getch();

		if (norm_denominator != 0)
		{
			V_norm = V_norm * (norm_numberator/norm_denominator);
		}
		else
		{
			printf("error in computing V-norm, ",
			"press any key to continue\n");
			getch();
			goto no_log;
		}

printf("e^vnorm[%i] = %e\n",r,V_norm);

	}

	V_norm = log(V_norm);


	printf("last V norm = %e, current v norm = %e\n",last_V_norm,V_norm);
getch();
	if(last_V_norm - V_norm <= convergence_criteria)
	{
	goto no_log;
		return;
	}

	last_V_norm = V_norm;
no_log:
	iteration ++;
printf("iteratation #%i\n",iteration);
getch();
	goto iterate;




/*	plot_all(M,&Y,"MEDD filter values, f");*/

}

init_pl()
{
	extern int plbuf[3072];

	int dmachn = 0;    	/* Use channel 0 for Register Mode transfers. */
	int port = 0x0318;	/* Port hexadecimal 0x0318. */




    plinit(dmachn, plbuf, sizeof(plbuf));
    plslib("\\c\\eellib\\pllib.15");
    plsprc(0, port, 0, 1024);



    plflsh(1);
    plflsh(-1);
}


/* This routine reads rows in, does real fft's or ifft's, and writes out. */

fft(int flag, int num, int num_rows, int n2p, double *in_data, double *out_data)
{
	unsigned int pl_array_a;
	unsigned int pl_array_b;

	int i;
/*print_complex(M,in_data);*/

	pl_array_a = 0x1000;
	pl_array_b = pl_array_a + 2*4*num;


				/* transfer the data array into
				   the PL processor.  Then convert
				   the IEEE double format to PL double
				   format. */

	plxfto(in_data, pl_array_b, num << 1);
	plwtxf();                            	/* wait on xfer */
	vfie64(pl_array_b, pl_array_a, num);
	plwtrn();				/* wait on run */
/*goto skip2;*/

					      /* Do the FFT or IFFT as
						 indicated by flag.
						 FFT on real 32 bit double
						 first 16 bits considered
						 Re part, last 16 bits are
						 the complex values. */

    if(flag > 0)
	rfftrw(pl_array_a, n2p, num_rows);            /* Do FFT for flag > 0. */
    else
	riftrw(pl_array_a, n2p, num_rows);            /* Do IFFT for flag <= 0. */

	plwtrn();                                  /* Wait on run */

						/* convert from PL doubles
						   to IEEE doubles and
						   xfer into the array
						   data
						*/
/*skip2:*/
	vtie64(pl_array_a, pl_array_b, num);
	plwtrn();
	plxffm(pl_array_b, out_data, num << 1);
	plwtxf();

/*print_complex(M,out_data);*/

/*	if (flag > 0)
	{
	for (i = num/2; i < num; i ++)
		*(out_data + i) = 0;
	}
	*/
    return;
}
/*print_complex(int lngth, double *data)
{
    int i;

    for (i = 0; i < lngth*2; i = i + 2)
	printf("data[%i] = ( %e , %e ) \n",i,*(data + i),*(data + i + 1));

	getch();


}
	*/
get_data(int lngth, complex *in_data)
{


	int i, j;
	double source_wavelet[M], spike_trace[M];
	complex b[M],sum_of_Bs,B[M],F[M],f[M];

	double real_i, real_j, sum;
	double ave_noise,ave_signal;

	double alpha,phi,time;
	double decay_rate_0,decay_rate_1;

	int reflect_1, reflect_2, reflect_3,reflect_4;

	complex noise[M];
	double sum_of_sums = 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;

	double noise_level = .001;

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

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

	for (j = 0; j < lngth; j ++)
	{
		real_j = j;
		time = real_j/lngth;
			noise[j].Re = (random(1000)/1000.0 - .5);
			noise[j].Im = (random(1000)/1000.0 - .5);

	}

/*	fft(-1,row_lngth,nrows,n2p,&noise[0],&noise[0]);*/

/*	alpha = -18.778;
	phi   = -.1960;
	decay_rate = 1.46;*/

	alpha = 240;
	phi   = 0;
	decay_rate_0 = 40;
	decay_rate_1 = -10.0;

	alpha = 175;
	phi   = 0;
	decay_rate_0 = 20;
	decay_rate_1 = -4.0;

	for (j = 0; j < lngth; j ++)
	{
		(in_data + j) -> Re = 0.0;
		(in_data + j) -> Im = 0.0;
	}

	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) * init_amp_0;

	if (source_wavelet[j] < 0)
	{
		source_wavelet[j] = source_wavelet[j] *
			exp(-decay_rate_0 * time);
	}
	else
	{
		source_wavelet[j] = source_wavelet[j] *
			((decay_rate_1 * time) + init_amp_0);
	}

		if (j > 23) source_wavelet[j] = 0.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 = 0; i <= j; i ++)
		{

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


			sum = sum + (sum * noise_level * .001) * noise[j].Re;
		}


		((in_data + j)->Re) = sum;
		((in_data + j)->Im) = 0.0;
/*		((in_data + j)->Re) = 1.0;
		((in_data + j)->Im) = 1.0;*/

		sum_of_sums = sum_of_sums + sum;

		b[j].Re = source_wavelet[j];
		b[j].Im = 0.0;

	}

	for (j = 0; j < lngth; j ++)
		((in_data + j)->Re) = ((in_data + j)->Re) / sum_of_sums;


	plot_all(lngth,&source_wavelet,"Source wavelet, w(i)");
	plot_all(lngth,&spike_trace,"Reflection signals, q(i)");
	plot_complex_data(1,lngth,in_data,"Input Trace");

	plot_complex_data(1,lngth,&b[0],"b");

	plot_complex_data(0,lngth,&noise[0],"noise");

	bp(1.0,1.0);

		fft(1,row_lngth,nrows,n2p,&b,&B);
	bp(1.0,2.0);

for (i = 0; i < lngth; i++)
	printf("(%e,%e) ",B[i].Re,B[i].Im);

	bp(1.0,3.0);

	plot_complex_data(3,lngth/2,&B,"B~B");

/*	for (j = 0; j < lngth; j ++)
	{
		complex_equate(&sum_of_Bs,
		  complex_add(tmp,
	*/
}
/* debug start
r = 0;
a = 4;
printf("check: norm numberator, denominator and norm itself (1st crack)\n");

		for (l = 0; l < M; l ++)
		{
			os = r*M + l;

			((x + os) -> Re) = l+2.0;
			((x + os) -> Im) = l+1.0;
			u1[r].Re = 3.0, u1[r].Im = 0.0;
			u2[r].Re = 2.0, u2[r].Im = -1.0;

			R[r][l].Re = l+1, R[r][l].Im = 0;
			G[r][l].Re = l+2, G[r][l].Im = l+1;

			X[r][l].Re = l+2, X[r][l].Im = l+1;
			Y[r][l].Re = l+2, Y[r][l].Im = l+1;
			F[r][l].Re = l+2, F[r][l].Im = l+1;
		}
u1[r].Re = 1/1135888.0, u1[r].Im = 0.0;
/* debug end            */
bp(float i,float j)
{
printf("%f %f\n ",i,j);
getch();
}