#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 1
	#define row_lngth M




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;
	int mode = 1;
	float real_row_lngth = row_lngth;
	float real_n2p;



	init_pl();

	getch();

	real_n2p = log(row_lngth)/log(2);
	n2p = real_n2p;

	printf("n2p = %f = %i\n",real_n2p,n2p);
	getch();


goto skip1;

	junk_data(M,&output[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]);

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

	junk_data(M,&data[0][0]);
	fft(-1,row_lngth,nrows,n2p,&output[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:
goto skip;
	get_data(M,0,&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,0,&data[0][0]);

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

for (i = 0; i < N; i ++)
	get_data(M,i,&data[i][0].Re);


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


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

    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[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 = {4.0,0.0};          /* Gray's alpha */
	complex a_minus_1;

	int r0 = 0;
	complex u1[N];
	complex u2[N];

	double   sum_of_xs;
	complex sum_of_Rs;
	complex sum_of_Fs;

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

	complex one_over_M = {1.0/M,0.0};

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

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

printf("1/m = (%f,%f)\n",one_over_M.Re,one_over_M.Im);

			/* Initialize filter estimate f^0 */


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

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

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


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

	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,n2p,(y + r*M),&Y[r][0].Re);

	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, 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].Re,
				complex_times(tmp,&F[w].Re, &Y[r][w].Re));
/*			X[r][w].Re = F[w].Re * Y[r][w].Re *
					(cos(F[w].Im - Y[r][w].Im));
			X[r][w].Im = 0;*/

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

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


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

	plot_complex_data(1,M,(x + r*M),
	"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].Re,
			   complex_add(tmp+2, &u1[r].Re,
				complex_exponent(tmp+1,
					complex_abs(tmp,(x + os)),&a.Re)));


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

		}

					/*      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,n2p,(x + r*M),&X[r][0].Re);


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



	plot_complex_data(1,M,(x + r*M),
	"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[w].Re,
		    complex_add(tmp+3,&R[w].Re,
		      complex_times(tmp+2,&u1[r].Re,
			complex_times(tmp+1,&Y[r][w].Re,
			  complex_conjagate(tmp,&Y[r][w].Re)))));

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

/*printf("u2 = (%e,%e)\n",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[N-1].Re = 0, G[N-1].Im = 0;*/
G[M/2].Re = 0, G[M/2].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[w] */

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

/*printf("R[w] = (%e,%e), sum of R from 0 to %i = (%e,%e)\n",
	R[w].Re,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[w] */

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



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

/*printf("k = (%e,%e), F[%i] = (%e,%e)\n",k.Re,k.Im,w,F[w].Re,F[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.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");

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

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

	plot_complex_data(1,M,&f[0].Re,
	"the 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.Re,
			  complex_abs(tmp,(x + os)));

/*printf("X = (%e,%e), |X| = %e ",(x + os)->Re,(x + os)->Im,
			x_tmp.Re);                      */

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

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

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

		}

		norm_numberator =
			pow((one_over_M.Re * norm_numberator),M/a.Re);

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




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

	num = num * 2;			/* there are num real and
					   num imaginary elements
					   makeing num * 2 elements */

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


					      /* 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)
	cfft(pl_array_a, n2p);            /* Do FFT for flag > 0. */
    else
	cifft(pl_array_a, n2p);            /* Do IFFT for flag <= 0. */

	plwtrn();                                  /* Wait on run */

						/* convert from PL doubles
						   to IEEE doubles and
						   xfer into the array
						   data
						*/

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

/*	if(flag > 0)
	{
	for (i = num/2; i < num; i ++)
	{
		*(out_data + i) = 0;
	}
	}*/

    return;
}

get_data(int lngth, int channel, 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];
	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[ 2] =  1.0;
		spike_trace[ 7] =  .5;
		spike_trace[12] =  -.3;
		spike_trace[19] =  .9;
		spike_trace[28] =  .5;

	}

	if(channel == 2)
	{
		spike_trace[10] = -.300;
		spike_trace[12] =  .582;
		spike_trace[30] =  -1.900;
	}

	if(channel == 3)
	{
		spike_trace[10] = .200;
		spike_trace[18] =  .582;
		spike_trace[22] =  .900;
		spike_trace[32] =  .440;
	}

	if(channel == 4)
	{
		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 */

	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;


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

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

	plot_complex_data(3,lngth,&B,"B~B");
	print_complex(11,&B);

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

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

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


	plot_complex_data(3,lngth,&F[0],"F~F");

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

	plot_complex_data(1,lngth,&f,"inverse wavelet, f");
		fft(1,row_lngth,nrows,n2p,in_data,&B[0]);


		for (j = 0; j < lngth; 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,n2p,&B[0],&f[0]);

	plot_complex_data(0,lngth,&f,"f = fft(FY)");

	print_complex(11,in_data);
	print_complex(30,&f);
}
/* 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, complex *f)
{
	int i;

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

}

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

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

}