suharlan.c

su 的源代码库

mean index	Index of the mean value sample in data histogram 
niter		Number of iterations for golden search (usually 10 will suffice)
pd		Array [nsamples] of data probability density function
pn		Array [nsamples] of noise probability density function

Output:
ps		Array [nsamples] of signal probability density function
*******************************************************************************
Note:
The deconvolved signal probability density funcion is computed by an
optimization algorithm in which the "best" ps(x) is computed to maximize the
probability (the fit) of the data histogram via minization of a suitable
cross-entropy indicator. The minization is performed for ps(x) with 
constraints of positivity adn unit area.
*******************************************************************************
Credits:
Translated to C by Gabriel Alvarez (1995) from a FORTRAN IV version
	written by Dr. Bill Harlan (1982)   
******************************************************************************/
{
	int j;			/* loop counter */
	int ndiv=20;		/* number of iteration for golden search */
	int iter;		/* index of iterations */
	int idiv;		/* index of golden search */
	float r1,r2;		/* proportions for ps updating */
	float rmax=1000.;	/* clip value for division of histograms */
	float alpha=0.5;	/* optimum golden search value */
	float xvalue=FLT_MIN;	/* optimum golden search x-value */
	float fvalue;		/* value of cross-entropy estimator */
	float scale=1.0;	/* some sort of scaling factor */
	float xgold;		/* x-value obtained from the golden search */
	float agold;		/* alpha obtained from the golden search */
	float sum;
	float *psp;		/* array of trial perturbed ps */
	float *grad;		/* array of computed gradient */

	/* allocate working space */
	psp = alloc1float(nsamples);
	grad = alloc1float(nsamples);

	/* initialize ps array to a spike at the mean amplitude of histograms */
	for (j=0; j<nsamples; j++) ps[j]=0.0;
	ps[mean_index]=1.0;

	/* main loop */
	for (iter=0; iter<niter; iter++) {
		
		/* scale down alphas if last search was too large 
		xvalue should have the smallest non-zero value previously used */
		scale = 2*alpha;
		if (alpha==0.0) scale = 2*xvalue;
		if (scale>1.0) scale = 1.0;

		/* compute normalized gradient of the cross-entropy estimator */
		gradient (nsamples, mean_index, rmax, pd, pn, ps, grad);

		/* do golden search to come up with optimum alpha value */
		idiv=1;
		while (idiv > 0 /*TRUE*/) {
			golden_search (fvalue, &idiv, &xgold, &agold);

			/* update xvalue and alpha */
			xvalue = xgold*scale;
			alpha = agold*scale;

			/* test to check loop exiting condition */
			if (idiv>=ndiv) break;

			/* compute trial ratios */
			r1 = 1.0 - xvalue;
			r2 = xvalue;

			/* update trial ps using trial ratios r1 and r2 */
			sum=0.0;
			for (j=0; j<nsamples; j++) {
				psp[j]=r1*ps[j]+r2*grad[j];
				sum +=psp[j];
			}

			/* normalize updated ps */
			for (j=0; j<nsamples; j++) psp[j] /=ABS(sum);

			/* compute cross-entropy estimator for updated ps */
			cross_entropy (nsamples, mean_index, rmax, pd, pn, 
				psp, &fvalue);
		}

		/* compute optimum ratios from optimum alpha */
		r1 = 1.0 - alpha;
		r2 = alpha;

		/* update ps using optimum ratios r1 and r2 */
		sum=0.0;
		for (j=0; j<nsamples; j++) {
			ps[j]=r1*ps[j]+r2*grad[j];
			sum +=ps[j];
		}

		/* normalize updated ps */
		for (j=0; j<nsamples; j++) ps[j] /=ABS(sum);
	}

	/* free allocated space */
	free1float(grad);
	free1float(psp);
}

/******************************************************************************

	Subroutine to compute the gradient of the cross-entropy estimator 
	  for constrained deconvolution of probability density
				functions

******************************************************************************/
void gradient (int ns, int si, float rmax, float *pd, float *pn,
	float *ps, float *grad)
/******************************************************************************
Input parameters:
ns		number of samples (intervals) in histograms
si		sample index of mean value in histograms
rmax		clip value for division of histograms
pd		data probability density function
pn		noise probability density function
ps		estimate of signal probability density function

Output value
grad		array of gradient of cross-entropy estimator
*******************************************************************************
Credits:
Translated to C by Gabriel Alvarez (1995) from a Fortran IV  version 
	written by Dr. Bill Harlan (1983)
******************************************************************************/
{
	int j;			/* loop variable */
	float sum;		/* variable to hold sum */
	float *pspn;		/* array for convolution of ps and pn */
	
	/* allocate space */
	pspn = alloc1float(ns);

	/* do convolution of pn with estimate of ps */
	conv (ns, -si, pn, ns, -si, ps, ns, -si, pspn); 

	/* divide pd and pspn term by term */
	for (j=0; j<ns; j++)
		pspn[j] = divide (rmax, pd[j], pspn[j]);

	/* cross-correlate pn with  pd/(ps*pn) */
	xcor (ns, -si, pn, ns, -si, pspn, ns, -si, grad); 

	/* normalize gradient */
	sum=0.0;
	for (j=0; j<ns; j++) sum +=grad[j];
	for (j=0; j<ns; j++) grad[j] /=ABS(sum);

	/* free allocated space */
	free1float(pspn);
}

/******************************************************************************

	Subroutine to compute cross-entropy estimator 

******************************************************************************/
void cross_entropy (int ns, int si, float rmax, float *pd, float *pn,
	float *ps, float *fvalue)
/******************************************************************************
Input parameters:
ns		number of samples (intervals) in histograms
si		sample index of mean value in histograms
rmax		clip value for division of histograms
pd		data probability density function
pn		noise probability densityu function
ps		estimate of signal probability function

Output value
fvalue		pointer to computed cross-entropy value
*******************************************************************************
Credits:
Translated to C by Gabriel Alvarez (1995) from a Fortran IV  version 
	written by Dr. Bill Harlan (1983)
******************************************************************************/
{
	int j;			/* loop counter */
	float sum=0.0;		/* to store the sum */
	float *pspn;		/* variable for convolution of ps and pn */

	/* allocate working space */
	pspn = alloc1float(ns);

	/* convolve ps and pn */
	conv (ns, -si, ps, ns, -si, pn, ns, -si, pspn);

	/* compute (ps*pn)/pd */
	for (j=0; j<ns; j++)
		pspn[j]=divide (rmax, pspn[j], pd[j]);
	
	/* compute the sum */
	for (j=0; j<ns; j++) {
		
		/* compute pd/(ps*pn) as the inverse of (ps*pn)/pd */
		pspn[j]=divide (rmax, 1.0, pspn[j]);

		/* multiply pd by log(pd/(ps*pn)) */
		sum +=pd[j]*log(pspn[j]);
	}

	/* output the result */
	*fvalue=sum;
} 

/******************************************************************************

	Subroutine to compute a=b/c  with a given clip value

******************************************************************************/
float divide(float rmax, float a, float b) 
/******************************************************************************
Input:
rmax		clip value
a		numerator
b		denominator

Output:	
		quotient=a/b
*******************************************************************************
Credits:
Translated to C by Gabriel Alvarez (1995) from a Fortran IV  version 
	written by Dr. Bill Harlan (1983)
******************************************************************************/
#define BIG 1000000000.0
{
	float quotient;			/* variable to hold the quotient a/b */

	/* check to see if both a and b are larger than BIG */ 	
	if (ABS(a)>BIG && ABS(b)>BIG) {
		a /=BIG;
		b /=BIG;
	} 
	
	/* check again */
	if (ABS(a)>BIG && ABS(b)>BIG) {
		a /=BIG;
		b /=BIG;
	} 

	rmax = ABS(rmax);		
	if (ABS(a)<(rmax*ABS(b))) {
		quotient  = a/b;
	
		/* this test will be true only if a and b are very small */
		if (ABS(quotient)>rmax) quotient=0.0;

	} else {
		quotient = rmax;

		/* check to see if quotient should be negative */
		if ((a>0.0 && b<0.0)||(a<0.0 && b>0.0)) quotient *= -1.0;
	}

	/* if numerator is zero, so is the quotient */	
	if (a==0.0) quotient=0.0;

	/* output quotient */
	return(quotient);
}

/******************************************************************************

	Subroutine to perform a one dimensional golden search
		for a functional value in the interval [0,1] 
	
******************************************************************************/
void golden_search (float fvalue, int *iter, float *xvalue, float *alpha)
/******************************************************************************
Input:
iter		index of iteration
fvalue		value of the function to find the minimum

Output:
xvalue		value of x that minimizes de function
alpha		optimum value
*******************************************************************************
Note:
The golden search algorithm used here is not optimum in the sense that it
does not make use of the ability of C to pass pointers to functions in a
subroutine. It was devised this way to cope with Fortran's lack of this
feature. A parabolic interpolation could perhaps be a better option.
*******************************************************************************
Credits:
Translated to C by Gabriel Alvarez (1995) from a Fortran IV  version 
	written by Dr. Bill Harlan (1983)
******************************************************************************/
{
	int j;				/* loop counter */
	static int ifill=0;		/* auxiliary function index */
	int im;				/* auxiliary index */
	float fm;			/* auxiliary variable */
	static float rmag=1.0;		/* golden ratio */
	static float x[4];		/* array to store searching points */
	static float f[4];		/* array to store function values */

	/* test condition for iteration index */
	if (*iter<1) *iter=1;			/* this should never happen */
	if (*iter==1) {
	
		/* define golden ratio */	
		*alpha = 0.0;
		rmag = 0.5*(sqrt(5.0)-1.0);

		/* initialize points for the search */
		x[0] = 0.0;
		x[1] = 1.0-rmag;
		x[2] = rmag;
		x[3] = 1.0;
		f[ifill] = fvalue;
		ifill = *iter-1;

	} else if (*iter<=4) {
	
		/* initialize function values */
		f[ifill] = fvalue;
		ifill = *iter-1;
	} else { 
	
		/* start the search */
		f[ifill] = fvalue;
		im = 0;
		fm = f[0];
		for (j=1; j<4; j++) {
			if (f[j]>=fm) continue;
			im = j;
			fm = f[j];
		}

		/* update bracketing interval */
		*alpha = x[im];
		if ((im==0)||(im==1)) {
			x[3] = x[2];
			f[3] = f[2];
			x[2] = x[1];
			f[2] = f[1];
			x[1] = (x[2] - x[0])*rmag +x[0];
			ifill = 1;
		} else {
			x[0] = x[1];
			f[0] = f[1];
			x[1] = x[2];
			f[1] = f[2];
			x[2] = x[3] - (x[3] - x[1])*rmag;
			ifill = 2;
		}
	}

	/* output values */
	*xvalue = x[ifill];
	*iter=*iter+1;
}

/* for graceful interrupt termination */
static void closefiles(void)
{
	efclose(headerfp);
	efclose(tracefp);
	eremove(headerfile);
	eremove(tracefile);
	exit(EXIT_FAILURE);
}