smooth.c

su 的源代码库

since it tends to honor them. If the peaks are too narrow with respect to the
filter width, however, some attenuation of the peaks occurs. Higher order 
(higher m) filters tend to perform better with narrower peaks at the expense of
the degree of smoothing of broader features. In general, the best results are
obtained when the full width (nl+nr+1) of the 4-degree filter is between 1 and
2 times the FWHM (full width at half of maximum) of the features that we
desire to smooth in the data.

Wider filters produce greater smoothing effect than narrower ones for filters
of a given order.

To within the roundoff error, the filter coefficients should add to one so that
the area under the curve is preserved after the smoothing. The sign of the 
samples is not preserved in general.

In order to compute filtered numerical derivatives (and not just to smooth 
the data), the parameter ld should be given the value of the desired
derivative (everything else is the same). 
	
A typical call to the convolution subroutine to apply the filter will be:
conv (nl+nr+1,-nl,filter,nd,0,data,nd,0,result);  where data[nd] is the input
data to smooth and result[nd] is the smoothed data

References:
Press, Teukolsky, Vettering, Flannery, Numerical Recipes in C:
		the art of scientific programming. Cambridge University
		Press. Second Edition (1992).
Ziegler, Applied Spectroscopy, Vol. 35, pp. 88-92.
Bromba and Ziegler, Analytical Chemistry, Vol. 53, pp 1582-1586.
***************************************************************************
Credits:
**************************************************************************
Credits
slightly modified from Numerical Recipes in C by Gabriel Alvarez (1995)
*************************************************************************/

{
	int imj,ipj,j,k,kk,mm;
	int *idx;			/* vector of indices */
	float d,fac,sum;		/* auxiliary variables */
	float **a;			/* matrix of .. */
	float *b;			/* vector of ... */

	/* defensive programming, check input */
	if (np < nl+nr+1 || nl < 0 || nr < 0 || ld > m || nl+nr < m)
		err(" error in smoothing function, check input parameters");

	/* allocate working space */
	idx = alloc1int(m+1);
	a = alloc2float(m+1,m+1);
	b = alloc1float(m+1);

	/* set up normal equations of the desired least-squares fit */
	for (ipj=0; ipj<=(m << 1); ipj++) {
		sum = (ipj ? 0.0 : 1.0);
		for (k=1; k<=nr; k++)
			sum +=pow((double)k,(double)ipj);
		for (k=1; k<=nl; k++)
			sum +=pow((double)-k,(double)ipj);
		mm = MIN(ipj,2*m-ipj);
		for (imj= -mm; imj<=mm; imj+=2)
			a[(ipj+imj)/2][(ipj-imj)/2]=sum;
	}

	/* solve normal equation via LU decomposition */
	LU_decomposition(m+1, a, idx, &d);

	/* initialize vector */
	for (j=0; j<=m; j++)
		b[j]=0.0;
	b[ld]=1.0;		/* right-hand-side vector is unit depending
					on what derivative we want */

	/* get one row of the inverse matrix */
	backward_substitution(m+1, a, idx, b);

	/* initialize output vector */
	for (kk=0; kk<np; kk++)
		filter[kk]=0.0;

	/* compute Savitzky-Golay coefficients */
	for (k= -nl; k<=nr; k++) {
		sum = b[0];
		fac=1.0;

		/* each coefficient is the dot product of powers of an 
		integer	with the inverse matrix row */
		for (mm=0; mm<m; mm++)
			sum +=b[mm+1]*(fac *=k);	
		/*kk=((np-k) % np);  store coefficients in wrap around order*/
		kk=k+nl;	/*store coefficients in normal order */
		filter[kk]=sum;	
	}

	/* free working space */
	free2float(a);
	free1float(b);
	free1int(idx);
}
/************************************************************************
	Subroutine to compute one-dimensional smoothing filter
		via running window average. 

************************************************************************/
void rwa_smoothing_filter (int flag, int nl, int nr, float *filter)
/************************************************************************
Input:
flag		=1 for rectangular window. =2 for triangular (weighted) window
nl		number of left (past) data points used
nr		number of right (future) data points used

Output:
filter		vector[nl+nr+1] of filter points to be convolved with the data

Notes:
The rectangular window should only be used when the data to be smoothed is
fairly smooth already; larger windows are not recommended unless extreme
smoothing is desired.

The triangular window will give more weight to the points that are closer
to the one to smooth. Although this is normally desirable, the degree of
smoothing for a given filter length is much less than with the rectangular
windom.

Both of these windows preserve the are under the curve and the sign of the 
samples to be smoothed.

A typical call to the convolution subroutine to apply the filter will be:
conv (nl+nr+1,-nl,filter,nd,0,data,nd,0,result);  where data[nd] is the input
data to smooth and result[nd] is the smoothed data
************************************************************************
Credits:
Gabriel Alvarez (1995).
************************************************************************/
{
	int i;				/* loop counter */
	int np=nl+nr+1;			/* number of filter points */

	if (flag==1) {
		float scale=1.0/np;	/* scale for rectangular window */

		for (i=0; i<np; i++) filter[i]=scale;

	} else if (flag==2) {
		float scale=0.0;	/* scale for triangular window */
		
		for (i= -nl; i<0; i++) filter[i+nl]=1+(float)i/nl;
		for (i=1; i<nr; i++) filter[i+nl]=1-(float)i/nr;
		filter[nl]=1.0;
	
		/* normalize */
		for (i=0;i<np;i++) scale+=filter[i];
		for (i=0;i<np;i++) filter[i] /=scale;
	} else err("error in rwa_smoothing filter, flag should 1 or 2\n");
}

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

		Subroutine to apply 2-D gaussian smoothing
	
******************************************************************************/	
void gaussian2d_smoothing (int nx, int nt, int nsx, int nst, float **data)
/******************************************************************************
Input:
nx		number of horizontal samples (traces)
nt		number of vertical (time) samples 
nsx		number of samples for horizontal smoothing (0 for no smoothing)
nst		number of samples for vertical smoothing (0 for no smoothing)
data		2-D array[nx][nt] of data to be smoothed

Output:
data		2-D array[nx][nt] of smoothed data
******************************************************************************/
{
	int ix,it;			/* loop counters */	
	float **temp;			/* temporary (scratch) array */
	
	/* allocate temporary space */
	temp=alloc2float(nx,nt);

	/* smooth over time */
	if (nst>0 && nt>1) {
		
		/* apply 1-D gaussian smoothing over time */
		for (ix=0; ix<nx; ix++)
			 gaussian1d_smoothing (nt, nst, data[ix]);
	}

	/*out1=fopen ("data/Stoz49","w");
	fwrite (data[49],sizeof(float),nt,out1);
	fclose(out1);*/

	/* if parameters are not reasonable, return */
	if (nsx<=0 || nx<=1) return;
	
	/* smooth over traces */
	for (it=0; it<nt; it++) {

		/* copy time smoothed data to temporary array */
		for (ix=0; ix<nx; ix++) temp[it][ix]=data[ix][it];

		/*if (it==300) {
			out1=fopen ("data/xoz300","w");
			fwrite(temp[300],sizeof(float),nx,out1);
			fclose(out1);
		}*/
	
		/* apply 1-D gaussian smoothing over traces */
		gaussian1d_smoothing (nx, nsx, temp[it]);

		/*if (it==300) {
			out1=fopen ("data/Sxoz300","w");
			fwrite(temp[300],sizeof(float),nx,out1);
			fclose(out1);
		}*/

		/* copy smoothed data back to output array */
		for (ix=0; ix<nx; ix++) data[ix][it]=temp[it][ix];
	}
	
	/* free allocated space */
	free2float(temp);
}

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

	Subroutine to apply a one-dimensional gaussian smoothing 

******************************************************************************/
void gaussian1d_smoothing (int ns, int nsr, float *data)
/******************************************************************************
Input:
ns		number of samples in the input data 
nsr		width (in samples) of the gaussian for which 
		amplitude > 0.5*max amplitude
data		1-D array[ns] of data to smooth

Output:
data		1-D array[ns] of smoothed data
******************************************************************************/
{
	int is;				/* loop counter */
	float sum=0.0;
	float fcut;
	float r;
	float fcutr=1.0/nsr;
	static int n;
	static int mean;
	static float fcutl;
	static float s[401];		/* smoothing filter array */
	float *temp;			/* temporary array */

	/* allocate space */
	temp=alloc1float(ns);

	/* save input fcut */
	fcut=fcutr;

	/* don't smooth if nsr equal to zero */
	if (nsr==0 || ns<=1) return;

	/* if halfwidth more than 100 samples, truncate */
	if (nsr>100) fcut=1.0/100;

	/* initialize smoothing function if not the same as the last one used */
	if (fcut !=fcutl) {
		fcutl=fcut;

		/* set span of 3, at width of 1.5*exp(-PI*1.5**2)=1/1174 */
		n=3.0/fcut+0.5;
		n=2*n/2+1;		/* make it odd for symmetry */

		/* mean is the index of the zero in the smoothing wavelet */
		mean=n/2;

		/* s(n) is the smoothing gaussian */
		for (is=1; is<=n; is++) {
			r=is-mean-1;
			r= -r*r*fcut*fcut*PI;
			s[is-1]=exp(r);
		}
			
		/* normalize to unit area, will preserve DC frequency at full 
		amplitude. Frequency at fcut will be half amplitude */
		for (is=0; is<n; is++) sum +=s[is];
		for (is=0; is<n; is++) s[is] /=sum;
	}

	/* convolve by gaussian into buffer */
	if (1.01/fcutr>(float)ns) {

		/* replace drastic smoothing by averaging */
		sum=0.0;
		for (is=0; is<ns; is++) sum +=data[is];
		sum /=ns;
		for (is=0; is<ns; is++) data[is]=sum;	

	} else {

		/* convolve with gaussian */
		conv (n, -mean, s, ns, -mean, data, ns, -mean, temp);

		/* copy filtered data back to output array */
		for (is=0; is<ns; is++) data[is]=temp[is];
	}

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

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

	Subroutine to apply gaussian smoothing to a histogram

******************************************************************************/
void smooth_histogram (int nintlh, float *pdf)
/******************************************************************************
Input:
nintlh		number of histogram bins
pdf		1-D array[nintlh] of input histogram to be smoothed

Output:
pdf		1-D array[nintlh] of smoothed histogram
******************************************************************************/
{
	int i;			/* loop counter */
	int ng=11;		/* number of samples in smoothing filter */
	int mg=6;		/* index of sample with zero lag */
	float rv=1.0;		/* ?? */
	float r;		/* auxiliary variable */
	float sum=0.0;		/* auxiliary variable */
	float *filter;		/* gaussian filter */
	float *spdf;		/* smoothed histogram */

	/* allocate working space */
	filter=alloc1float(ng);
	spdf=alloc1float(nintlh);

	/* compute gaussian filter */	
	for (i=1; i<=ng; i++) {
	
		r=i-mg;
		r= -r*r/rv/rv*PI;
		filter[i-1]=exp(r);	
	}

	/* apply filter */
	conv (ng, -mg+1, filter, nintlh, 0, pdf, nintlh, 0, spdf);

	/* normalize output histogram */
	for (i=0; i<nintlh; i++) sum +=spdf[i]; 
	for (i=0; i<nintlh; i++) pdf[i]=spdf[i]/sum; 

	/* free allocated space */
	free1float(filter);
	free1float(spdf);
}