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

						/* parameter definitions */
	#define n 1
	#define N n
	#define m 128
	#define M m/2 + 1
	#define nrows 1
	#define row_lngth m

	#define view_plots 0			 /* view  mode on if != 0 */
	#define debug 0				 /* debug mode on if != 0 */
	#define pl 1


typedef struct
{
	double Re;
	double Im;

} complex;

						/* data input array and
						   output arrays
						   they are actually
						   real arrays but will
						   be treated as complex
						   by the MED algorithm
						   FFT implmentation used
						   by Gray */
	complex data[n][m];
	complex output[n][m];
					       /* variables used by
						  the pl array processor
						  initialized in the pl_init
						  function */
	unsigned int pl_array_a;
	unsigned int pl_array_b;
	unsigned int pl_array_c;
	unsigned int pl_array_d;

	int plbuf[3072];

	int m2p;				/* exponent of 2 of the
						   number (m) of samples per
						   channel, initialized in
						   pl_init */

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;
	int mode = 1;

	complex imag[n][m],true_filter[m];


	init_pl();

	getch();

	printf("m2p = %i\n",m2p);
	getch();


goto skip1;

	junk_data(m,&data[0][0]);
/*	more_data(m,&data[0][0]);*/

	plot_complex_data(1,m,&data[0][0],"test data: x");
	print_complex(8,&data[0][0]);
	print_complex(m,&data[0][0]);

	fft(1,row_lngth,nrows,m2p,&data[0][0],&imag[0][0]);
	plot_complex_data(3,M,&imag[0][0],"FFT(x) from pl processor");
	print_complex(8,&imag[0][0]);
	print_complex(M,&imag[0][0]);

	fft(-1,row_lngth,nrows,m2p,&imag[0][0],&data[0][0]);
	plot_complex_data(1,m,&data[0][0],"IFFT(x) from pl processor");
	print_complex(8,&data[0][0]);
	print_complex(m,&data[0][0]);

		    goto end;
skip1:
skip:

for (i = 0; i < n; i ++)
	get_data(m,i,&data[i][0],&true_filter);


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


	plot_complex_data(1,m,&output,
	"MEDD Output, Y(i) = f * x(i)");
end:

    return;
}

deconvolve(complex *y, complex *x, complex *true_filter)
{

/*	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 real (double) matrix (address given by *y)

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

*/
	#define convergence_criteria .001


	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[m];
	complex F[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[m];		/* frequency domain of autocorrelations */
	complex G[m];		/* "                 " crosscorrelations cubed */


	complex a = {6.0,0.0};         /* Gray's alpha */
	complex a_minus_1;

	int r0 = 0;

	complex u1[n];

	complex u2[n];

	complex sum_of_xs;
	complex sum_of_Rs;
	complex sum_of_Fs;

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

	double one_over_m = 1.0/m;
	double m_over_a   = m/a.Re;
	double m_over_2   = m/2.0;

	complex one_over_M = {1.0/m, 0.0};

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

	complex x_tmp;

	double norm_numberator;
	double norm_denominator;

	double V_norm;
	double last_V_norm = 0.0;

	int iteration = 0;

	complex Two = {2.0, 0.0};
	complex minus_One = {-1.0,0.0};

	complex cx;


	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.Re - 1.0, a_minus_1.Im  = 0.0;

printf("a-1 = (%f) and 1/m = \n",a_minus_1.Re);
print_complex(1,&one_over_M);
getch();

			/* Initialize filter estimate f^0 */


	for (t = 0; t < m; t ++)
	{
		f[t].Re = 0.0, f[t].Im = 0;
/*		f[t].Re = (true_filter + t)->Re;
		f[t].Im = (true_filter + t)->Im;*/
	}


	f[t0].Re = 1.0;


	if (view_plots != 0 ) plot_complex_data(1,m,&f[0],
	"first approx. to filter f");


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


	if (view_plots != 0) plot_complex_data(2,m,&F[0].Re,
	"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[w].Re = 0, R[w].Im = 0;
			G[w].Re = 0, G[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 + r*m),
	"the input matrix, y");

	fft(1,row_lngth,nrows,m2p,(y + r*m),&Y[r][0]);

	plot_complex_data(2,m,&Y[r][0].Re,
	"the FFT of the input matrix, Y");

				/* provide zero padding on Y */

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

			y[r][m + l].Re = 0.0, Y[r][m + l].Im = 0.0;
		} wrong!! */

				/* calculate the output matrix X
				   using the approximation F */

if (pl != 1)
{
		for (w = 0; w < m; w ++)
		{
			complex_equate(&X[r][w],
				complex_times(tmp,&F[w],&Y[r][w]));
		}
}
else
{
		pc_to_pl(m, pl_array_a, &F[0]);
		pc_to_pl(m, pl_array_b, &Y[r][0]);

		cvmul(pl_array_a, pl_array_b, pl_array_c, m);
		pl_to_pc(m, pl_array_c, &X[r][0]);

}
/*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,&X[r][0].Re,
	"the FFT of the output matrix, X");

	printf("X = FY:\n");
	print_complex(10,&X[r][0].Re);


		fft(-1,row_lngth,nrows,m2p,&X[r][0].Re,(x + r*m));

	plot_complex_data(1,m,(x + r*m),
	"the output matrix, x");

	print_complex(10,(x + r*m));


		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;

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

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

/*		     complex_equate(&cx,
				complex_exponent(tmp+1,
					complex_abs(tmp,(x + os)),&a.Re));
print_complex(1,&cx);
		     complex_equate(&cx,
					complex_abs(tmp,(x + os)));
print_complex(1,&cx);
	*/
		     complex_equate(&u2[r].Re,
			   complex_add(tmp+1,&u2[r].Re,
				complex_exponent(tmp,(x + os),&Two.Re)));

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

		}
print_complex(1,&u1[r]);
print_complex(1,&minus_One);


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

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

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

						     /* 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("x = (%e,%e) os = %i "(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.Re,
			    complex_exponent(tmp+1,
			      complex_abs(tmp,(x + os)),&a_minus_1.Re)));

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

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

		fft(1,row_lngth,nrows,m2p,(x + r*m),&X[r][0].Re);



	plot_complex_data(1,m,(x + r*m),
	"the new estimation of output matrix, x");
	print_complex(10,(x + r*m));

	plot_complex_data(2,m,&X[r][0].Re,
	"the FFt of the output matrix, X");
	print_complex(10,&X);



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


		for (w = 0; w < m; w ++)
		{
							/*		*
							   R = R + (u1*Y  Y)
							*/
			os = r*m + w;

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

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

/*printf("u1[%i] = (%e,%e)\n",r,u1[r].Re,u1[r].Im);*/
/*printf("R[%i] = (%e,%e)\n",w,R[w].Re,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[w],
		    complex_add(tmp+3,&G[w],
		      complex_times(tmp+2,&u2[r],
			complex_times(tmp+1,(x + os),
			  complex_conjagate(tmp,(y+os))))));

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


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

/*G[0].Re = 0, G[0].Im = 0;*/
/*G[M-1].Re = 0, G[M-1].Im = 0;*/
/*G[M].Re = 0, G[M].Im = 0;*/


				/* process next trace */
	}

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

					/* get unnormalized F[w] */

			pc_to_pl(m, pl_array_a, &G[0]);
			pc_to_pl(m, pl_array_b, &R[0]);

			cvdivr(pl_array_a, pl_array_b, pl_array_c, m);

			pl_to_pc(m, pl_array_c, &f[0]);



	plot_complex_data(1,m,&f[0],
	"the unnormalized filter matrix, f");

printf("unnormalized F =\n");

print_complex(10,&f);

					/* 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.Re,
			  complex_add(tmp+2,&sum_of_Fs.Re,
			    complex_times(tmp+1,&f[w].Re,
			      complex_conjagate(tmp,&f[w].Re))));

/*printf("F[%i] = (%e,%e), sum = (%e,%e)\n",
	w,f[w].Re,f[w].Im,sum_of_Fs.Re,sum_of_Fs.Im);
getch();*/
		}
					/* normalize f[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[w].Re,
				complex_divide(tmp,&f[w].Re,&sum_of_Fs.Re));

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

	getch();

/*	plot_complex_data(2,m,&F[0].Re,
	"the FFt of the filter matrix, F");
	print_complex(10,&F);

	fft(-1,row_lngth,nrows,m2p,&F[0].Re,&f[0]);

	printf("F = \n");
	print_complex(m,&F[0].Re);*/

	plot_complex_data(1,m,&f[0],
	"the filter matrix, f");

	fft(1,row_lngth,nrows,m2p,&f[0].Re,&F[0]);


	plot_complex_data(1,m,&F[0],
	"the fft of filter matrix, F");



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

	for (r = 0; r < n; r ++)
	{
		norm_numberator  = 0;
		norm_denominator = 0;

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

			complex_equate(&x_tmp,
			  complex_abs(tmp,(x + os)));

/*printf("norm rop: %e\n",norm_numberator);
printf("X = (%e,%e), |X| = (%e,%e) |x|^a = %e\n",(x + os)->Re,(x + os)->Im,
			x_tmp.Re,x_tmp.Im,exponent(x_tmp.Re,a.Re));
			getch();*/

			norm_numberator =
			 norm_numberator + exponent(x_tmp.Re,a.Re);

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


/*printf("os = %i, sum of |X|^a = %e\nsum of |X|^2 = %e\n",
	os, norm_numberator,norm_denominator);
if (norm_numberator < 0) getch();*/

		}
printf("total: sum of |X|^a = %e\nsum of |X|^2 = %e\n",
	norm_numberator,norm_denominator);
getch();

		norm_numberator =
			log (one_over_m * norm_numberator) * (m_over_a);

		norm_denominator =
			log(one_over_m * norm_denominator) * (m_over_2);

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

		if (norm_denominator != 0)
		{
			V_norm = V_norm * exp(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)
	{
	printf("would have bailed out!\n");

	goto no_log;
		return;
	}

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





}

init_pl()
{
	extern int plbuf[3072];

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

	float real_m2p;

				/* initialize the exponent of two of m
				   samples per channel */

	real_m2p = log(row_lngth)/log(2);
	m2p = real_m2p;

				/* initialize the pl processor array
				   addresses for 64 bit (double) complex
				   values of length m */

	pl_array_a = 0x1000;
	pl_array_b = pl_array_a + 2*4*m;
	pl_array_c = pl_array_b + 2*4*m;
	pl_array_d = pl_array_c + 2*4*m;



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

pc_to_pl(int num, unsigned int pl_array, double *in_data)
{

				/*   Transfer the array in_data into
				   the PL processor.  Then convert
				   the IEEE-64 double format to PL double
				   format. There are m * 2 32 bit words
				   since the array is double (64 bits/word)
				   and since there are two 64 bit words
				   (one comprising the real part, et al)
				   then there are m*2*2 32 bit words to xfer
				   (ie num << 2).
				     Once in the pl_processor there are 2*m
				   64 bit words to convert to pl native 32
				   bit words (one word for real part, et al)
				     The same logic is used to return the
				   processed array (in-data) back to the
				   host array (out_data).
				   */


	plxfto(in_data, pl_array_d, num << 2);
	plwtxf();                            	/* wait on xfer */
	vfie64(pl_array_d, pl_array, num << 1);
	plwtrn();				/* wait on run */

    return;
}

pl_to_pc(int num, unsigned int pl_array, double *out_data)
{

				/* convert from PL native
				   to IEEE-64 doubles and
				   xfer into the array
				   out_data
				*/

	vtie64(pl_array, pl_array_d, num << 1);
	plwtrn();
	plxffm(pl_array_d, out_data, num << 2);
	plwtxf();


    return;
}

fft(int flag, int num, int num_rows, int m2p, double *in_data, double *out_data)
{

	int i;
	complex scalar;


				/*   Transfer the array in_data into
				   the PL processor.  Then convert
				   the IEEE-64 double format to PL double
				   format. There are m * 2 32 bit words
				   since the array is double (64 bits/word)
				   and since there are two 64 bit words
				   (one comprising the real part, et al)
				   then there are m*2*2 32 bit words to xfer
				   (ie num << 2).
				     Once in the pl_processor there are 2*m
				   64 bit words to convert to pl native 32
				   bit words (one word for real part, et al)
				     The same logic is used to return the
				   processed array (in-data) back to the
				   host array (out_data).
				   */


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


				      /* Do the FFT or IFFT as
					 indicated by flag.
					 FFT on complex 32 bit pl
					 native floats, the first
					 32 bits are considered the
					 complex part and the
					 last 32 bits are the
					 imaginary part.
					 */

    if(flag > 0)
    {
	cfft(pl_array_a, m2p);  	        /* Do FFT for flag > 0. */
	cvconj(pl_array_a, pl_array_a, num);    /* take complex conjagate */
    }
    else
    {
	cvconj(pl_array_a, pl_array_a, num);	/* take complex conjagate */
	cifft(pl_array_a, m2p);            	/* Do IFFT for flag <= 0. */
/*	cereal(pl_array_a, pl_array_b, num);
	cvreal(pl_array_b, pl_array_a, num);*/
    }

	plwtrn();                               /* Wait on run */


				/* scale the data
				   by root(n)
				*/


	scalar.Re = sqrt(num);                 	/* if ifft the scale factor
	scalar.Im = 0;				   is root(m) */

	if (flag > 0)
	{
		scalar.Re = 1/sqrt(num);        /* otherwise, the scale
		scalar.Im = 0;                     factor is 1/root(m) */
	}

				/* convert the scaling factor in scalar
				   from IEEE-64 doubles to the pl_processor
				   native mode for scaling.
				*/

	plxfto(&scalar, pl_array_b, 4);         /* xfer to pl processor */
	plwtxf();                            	/* wait on xfer */
	vfie64(pl_array_b, pl_array_b, 2);
	plwtrn();				/* wait on run */

				/* perform the multiplication
				   of scalar on complex
				   array in a, place in a,
				*/

	cvcsml(pl_array_b, pl_array_a, pl_array_a, num);

				/* convert from PL native
				   to IEEE-64 doubles and
				   xfer into the array
				   out_data
				*/

	vtie64(pl_array_a, pl_array_b, num << 1);
	plwtrn();
	plxffm(pl_array_b, out_data, num << 2);
	plwtxf();


    return;
}

get_data(int lngth, int channel, complex *in_data, complex *true_filter)
{


	int i, j;
	double source_wavelet[m], spike_trace[m];
	complex sum_of_Bs,B[m],F[m];
	complex b[m],f[m];
	complex Tmp[5], *tmp = &Tmp[0];

	complex B_tmp;


	double real_i, real_j, sum;
	double ave_noise,ave_signal;
	double max_sig = 0;

	double alpha = 127.0;
	double phi = 0;
	double time;
	double decay_rate = -20.0;

	complex noise[m];
	double sum_of_sums = 0;

	double init_amp_0 =  1.0;

	double noise_level = .001;

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

	}

	for (j = 0; j < lngth; j ++)
	{
		(in_data + j) -> Re = 0.0;
		(in_data + j) -> Im = 0.0;
		source_wavelet[j] = 0;
		spike_trace[j] = 0;
	}
	/* create source wavelet (Cabrelli ex. 1) */

		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;

/*	for (j = 0; j < lngth; j ++)
	{
		real_j = j;
		time = real_j/lngth;
		source_wavelet[j] =
		sin(alpha * time + phi) * exp(decay_rate * time) * init_amp_0;
	}
*/

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

	if(channel == 0)
	{
/*		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;
*/
		spike_trace[ 2] = 1.00;
		spike_trace[ 7] =  .5;
		spike_trace[12] =  -.3;
		spike_trace[19] =   .9;
		spike_trace[28] =  .5;
	}

	if(channel == 1)
	{
/*		spike_trace[ 5] = .6;
		spike_trace[10] =  .5;
		spike_trace[15] =  -.900;
		spike_trace[20] =  1.0;
		spike_trace[40] =  -.420;
		spike_trace[45] = -.10;
*/
		spike_trace[ 3] =  .6;
		spike_trace[ 9] =  .3;
		spike_trace[11] =  .3;
		spike_trace[14] =  .3;
		spike_trace[28] =  .5;

	}

	if(channel == 2)
	{
		spike_trace[ 6] =  1.3;
		spike_trace[17] =  -.1;
	}

	if(channel == 3)
	{
		spike_trace[ 2] =  .1;
		spike_trace[ 5] =  -.2;
		spike_trace[10] =  .3;
		spike_trace[13] =  -.4;
		spike_trace[22] =  .5;
	}

	if(channel == 4)
	{
		spike_trace[ 4] =  2.3;
		spike_trace[ 7] =  .1;
		spike_trace[17] =  -.3;
		spike_trace[21] =  .3;
		spike_trace[28] =  -.5;
	}
	if(channel == 5)
	{
		spike_trace[ 2] =  2.3;
		spike_trace[ 8] =  .5;
		spike_trace[16] =  -.7;
		spike_trace[19] =  .3;
		spike_trace[24] =  -1.5;
	}


	/* create convolution of source and spike */

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

/*printf("source wavelet[%i] = %f\n",j,source_wavelet[j]);
if (source_wavelet[j] != 0.0) getch();*/

		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)->Re) = sum;
		((in_data + j)->Im) = 0.0;


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

		max_sig = max(max_sig,sum);

	}

	sum_of_sums = 0;

/*	for (j = 0; j < lngth; j ++)
	{
		ave_signal = ((in_data + j)->Re);
		ave_noise = (max_sig * noise_level) * noise[j].Re;
		((in_data + j)->Re) = ave_signal + ave_noise;

/*	printf("s = %f, n = %f, max(z) = %f, s/n = %f\n",
	ave_signal,ave_noise,max_sig,ave_signal/ave_noise);
	getch();*/
/*
		sum_of_sums = sum_of_sums + ((in_data + j)->Re);
	}

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

	if (view_plots != 0)
	{

	plot_real_data(0,lngth,&source_wavelet,"Source wavelet, w(i)");
	plot_real_data(0,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");*/

	}

		fft(1,row_lngth,nrows,m2p,&b[0],&B[0]);

	if (view_plots != 0)
		plot_complex_data(3,lngth/2,&B,"B~B");

	sum_of_Bs.Re = 0, sum_of_Bs.Im = 0;

	for (j = 0; j < m; j ++)
	{
		complex_equate(&sum_of_Bs,
		  complex_add(tmp+2,&sum_of_Bs,
		    complex_times(tmp+1,&B[j],
		      complex_conjagate(tmp,&B[j]))));
	}

	if (view_plots != 0)
		printf("sum(B) = %f\n",sum_of_Bs.Re);

	sum_of_Bs.Re = sum_of_Bs.Re * 1.0e-08;
	sum_of_Bs.Im = sum_of_Bs.Im * 1.0e-08;

	for (j = 0; j < m; j ++)
	{
		complex_equate(&B_tmp,
		  complex_conjagate(tmp,&B[j]));

		complex_equate(&F[j],
		  complex_divide(tmp+2,&B_tmp,
		    complex_add(tmp+1,&sum_of_Bs,
		      complex_times(tmp,&B[j],&B_tmp))));
	}


	if (view_plots != 0)
		plot_complex_data(3,m,&F[0],"F~F");

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

	if (view_plots != 0)
		plot_complex_data(1,lngth,&f,"inverse wavelet, f");

		fft(1,row_lngth,nrows,m2p,in_data,&B[0]);

		for (j = 0; j < m; j ++)
		{
			complex_equate((true_filter + j),&f[j]);
		}

		for (j = 0; j < m; j ++)
		{
			complex_equate(&B[j],
				complex_times(tmp,&F[j], &B[j]));

/*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();*/
		}

		fft(-1,row_lngth,nrows,m2p,&B[0],&f[0]);

	if (view_plots != 0)
		plot_complex_data(1,lngth,&f,"f = fft(FY)");





}
/* 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[l].Re = l+1, R[l].Im = 0;
			G[l].Re = l+2, G[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            */


data_fill(complex *f)
{
	complex X[m+1];
	int i;

/*	 This is the test data record used in the spectrum examples.
	from Marple "Digital Spectral Analysis with Applications",
	Appendix II
*/

	  X[ 1].Re =   1.349839091, X[ 1].Im =   2.011167288;
	  X[ 2].Re =  -2.117270231, X[ 2].Im =   0.817693591;
	  X[ 3].Re =  -1.786421657, X[ 3].Im =  -1.291698933;
	  X[ 4].Re =   1.162236333, X[ 4].Im =  -1.482598066;
	  X[ 5].Re =   1.641072035, X[ 5].Im =   0.372950256;
	  X[ 6].Re =   0.072213709, X[ 6].Im =   1.828492761;
	  X[ 7].Re =  -1.564284801, X[ 7].Im =   0.824533045;
	  X[ 8].Re =  -1.080565453, X[ 8].Im =  -1.869776845;
	  X[ 9].Re =   0.927129090, X[ 9].Im =  -1.743406534;
	  X[10].Re =   1.891979456, X[10].Im =   0.972347319;
	  X[11].Re =  -0.105391249, X[11].Im =   1.602209687;
	  X[12].Re =  -1.618367076, X[12].Im =   0.637513280;
	  X[13].Re =  -0.945704579, X[13].Im =  -1.079569221;
	  X[14].Re =   1.135566235, X[14].Im =  -1.692269921;
	  X[15].Re =   1.855816245, X[15].Im =   0.986030221;
	  X[16].Re =  -1.032083511, X[16].Im =   1.414613724;
	  X[17].Re =  -1.571600199, X[17].Im =   0.089229003;
	  X[18].Re =  -0.243143231, X[18].Im =  -1.444692016;
	  X[19].Re =   0.838980973, X[19].Im =  -0.985756695;
	  X[20].Re =   1.516003132, X[20].Im =   0.928058863;
	  X[21].Re =   0.257979959, X[21].Im =   1.170676708;
	  X[22].Re =  -2.057927608, X[22].Im =   0.343388647;
	  X[23].Re =  -0.578682184, X[23].Im =  -1.441192508;
	  X[24].Re =   1.584011555, X[24].Im =  -1.011150956;
	  X[25].Re =   0.614114344, X[25].Im =   1.508176208;
	  X[26].Re =  -0.710567117, X[26].Im =   1.130144477;
	  X[27].Re =  -1.100205779, X[27].Im =  -0.584209621;
	  X[28].Re =   0.150702029, X[28].Im =  -1.217450142;
	  X[29].Re =   0.748856127, X[29].Im =  -0.804411888;
	  X[30].Re =   0.795235813, X[30].Im =   1.114466429;
	  X[31].Re =  -0.071512341, X[31].Im =   1.017092347;
	  X[32].Re =  -1.732939839, X[32].Im =  -0.283070654;
	  X[33].Re =   0.404945314, X[33].Im =  -0.781708360;
	  X[34].Re =   1.293794155, X[34].Im =  -0.352723092;
	  X[35].Re =  -0.119905084, X[35].Im =   0.905150294;
	  X[36].Re =  -0.522588372, X[36].Im =   0.437393665;
	  X[37].Re =  -0.974838495, X[37].Im =  -0.670074046;
	  X[38].Re =   0.275279552, X[38].Im =  -0.509659231;
	  X[39].Re =   0.854210198, X[39].Im =  -0.008278057;
	  X[40].Re =   0.289598197, X[40].Im =   0.506233990;
	  X[41].Re =  -0.283553183, X[41].Im =   0.250371397;
	  X[42].Re =  -0.359602571, X[42].Im =  -0.135261074;
	  X[43].Re =   0.102775671, X[43].Im =  -0.466086507;
	  X[44].Re =  -0.009722650, X[44].Im =   0.030377999;
	  X[45].Re =   0.185930878, X[45].Im =   0.808869600;
	  X[46].Re =  -0.243692726, X[46].Im =  -0.200126961;
	  X[47].Re =  -0.270986766, X[47].Im =  -0.460243553;
	  X[48].Re =   0.399368525, X[48].Im =   0.249096692;
	  X[49].Re =  -0.250714004, X[49].Im =  -0.362990230;
	  X[50].Re =   0.419116348, X[50].Im =  -0.389185309;
	  X[51].Re =  -0.050458215, X[51].Im =   0.702862442;
	  X[52].Re =  -0.395043731, X[52].Im =   0.140808776;
	  X[53].Re =   0.746575892, X[53].Im =  -0.126762003;
	  X[54].Re =  -0.559076190, X[54].Im =   0.523169816;
	  X[55].Re =  -0.344389260, X[55].Im =  -0.913451135;
	  X[56].Re =   0.733228028, X[56].Im =  -0.006237417;
	  X[57].Re =  -0.480273813, X[57].Im =   0.509469569;
	  X[58].Re =   0.033316225, X[58].Im =   0.087501869;
	  X[59].Re =  -0.321229130, X[59].Im =  -0.254548967;
	  X[60].Re =  -0.063007891, X[60].Im =  -0.499800682;
	  X[61].Re =   1.239739418, X[61].Im =  -0.013479125;
	  X[62].Re =   0.083303742, X[62].Im =   0.673984587;
	  X[63].Re =  -0.762731433, X[63].Im =   0.408971250;
	  X[64].Re =  -0.895898521, X[64].Im =  -0.364855707;

	  for (i = 0; i < 64; i ++)
	  {
		(f + i) -> Re = X[i+1].Re;
		(f + i) -> Im = X[i+1].Im;
	  }

}
more_data(int lngth, double *f)
{
	int i;

	  for (i = 0; i < lngth; i ++)
	  {
		*(f + i) = 1;
	  }

}

junk_data(int lngth, complex *f)
{
	int i;

	  for (i = 0; i < lngth; i ++)
	  {
		(f + i) -> Re = 0;
		(f + i) -> Im = 0;
	  }
(f + 0)->Re = 1;
}