holbergd1.c

来自「该程序是用vc开发的对动态数组进行管理的DLL」· C语言 代码 · 共 134 行

/* Copyright (c) Colorado School of Mines, 2003.*/
/* All rights reserved.                       */

/*********************** self documentation **********************/
/*****************************************************************************
HOLBERG1D - Compute coefficients of Holberg's 1st derivative filter

holberg1d	comput the coefficients of Holberg's 1st derivative filter

******************************************************************************
Function Prototype:
void holbergd1 (float e, int n, float d[]);

******************************************************************************
Input:
e		maximum relative error in group velocity
n		number of coefficients in filter (must be 2, 4, 6, or 8)

Output:
d		array[n] of coefficients

******************************************************************************
Notes:
Coefficients are output in a form suitable for convolution.  The
derivative is centered halfway between coefficients d[n/2-1] and d[n/2].

Coefficients are computed via the power series method of Kindelan et al.,
1990, On the construction and efficiency of staggered numerical
differentiators for the wave equation:  Geophysics 55, 107-110.
See also, Holberg, 1987, Computational aspects of the choice of
operator and sampling interval for numerical differentiation in
large-scale simulation of wave phenomena:  Geophys. Prosp., 35, 629-655

******************************************************************************
Reference:
Kindelan et al., 1990, 
On the construction and efficiency of staggered numerical
differentiators for the wave equation:  Geophysics 55, 107-110.

See also, Holberg, 1987, Computational aspects of the choice of
operator and sampling interval for numerical differentiation in
large-scale simulation of wave phenomena:  Geophys. Prosp., 35, 629-655

******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 06/06/91
*****************************************************************************/
/**************** end self doc ********************************/

#include "cwp.h"

void holbergd1 (float e, int n, float d[])
/*****************************************************************************
Compute coefficients of Holberg's 1st derivative filter
******************************************************************************
Input:
e		maximum relative error in group velocity
n		number of coefficients in filter (must be 2, 4, 6, or 8)

Output:
d		array[n] of coefficients
******************************************************************************
Notes:
Coefficients are output in a form suitable for convolution.  The
derivative is centered halfway between coefficients d[n/2-1] and d[n/2].

Coefficients are computed via the power series method of Kindelan et al.,
1990, On the construction and efficiency of staggered numerical
differentiators for the wave equation:  Geophysics 55, 107-110.
See also, Holberg, 1987, Computational aspects of the choice of
operator and sampling interval for numerical differentiation in
large-scale simulation of wave phenomena:  Geophys. Prosp., 35, 629-655
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 06/06/91
*****************************************************************************/
{
	static float b[4][4][6] = {
	{
		{1.0f,1.0f,0.0f,0.0f,0.0f,0.0f},
		{0.0f,0.0f,0.0f,0.0f,0.0f,0.0f},
		{0.0f,0.0f,0.0f,0.0f,0.0f,0.0f},
		{0.0f,0.0f,0.0f,0.0f,0.0f,0.0f}
	},
	{
		{1.125f,0.4330f,-0.4583f,0.2566f,0.0f,0.0},
		{-0.04167f,-0.1443f,-0.1806f,-0.0855f,0.0f,0.0},
		{0.0f,0.0f,0.0f,0.0f,0.0f,0.0f},
		{0.0f,0.0f,0.0f,0.0f,0.0f,0.0f}
	},
	{
		{1.172f,0.2742f,-0.3006f,0.2637f,-0.2391f,0.0f},
		{-0.06510f,-0.1371f,-0.0100f,0.1077f,-0.0086f,0.0f},
		{0.004688f,0.0274f,0.0661f,0.0826f,0.0530f,0.0f},
		{0.0f,0.0f,0.0f,0.0f,0.0f,0.0f}
	},
	{
		{1.1963f,0.1987f,-0.2195f,0.2096f,-0.2038f,0.1779f},
		{-0.07975f,-0.1192f,0.03929f,0.04049f,-0.05841f,0.04389f},
		{0.009570f,0.03975f,0.04853f,-0.00967f,-0.05507f,-0.01117f},
		{-0.0006975f,-0.00568f,-0.02014f,-0.04039f,-0.04937f,-0.03625f}
	} };
	int i,j;
	
	/* compute coefficients via power series */
	for (i=0; i<n; ++i)
		for (j=0,d[i]=0.0; j<6; ++j)
			d[i] += (float)(b[n/2-1][i][j]*pow(e,2.0*j/n));
			
	/* arrange coefficients for convolutional derivative filter */
	for (i=n/2; i<n; ++i)
		d[i] = -d[i-n/2];
	for (i=0; i<n/2; ++i)
		d[i] = -d[n-i-1];	
}

#ifdef TEST
void main()
{
	int i,n;
	char s[100];
	float e,d[100],dd[100];
	
	while(1) {
		printf("Enter e, n:  ");
		gets(s);  sscanf(s,"%f %d",&e,&n);
		holberg1(e,n,d);
		for (i=0; i<n; ++i)
			printf("d[%d]=%g\n",i,d[i]);
		conv(n,0,d,n,0,d,n*2-1,0,dd);
		for (i=0; i<n*2-1; ++i)
			printf("dd[%d]=%g\n",i,dd[i]);
	}
}
#endif /* TEST */