/* filtered.c
```********************************************************************************
* filtered.c -- Routines to filter a statistical series objects
*
* User routines
*	Pro_FilteredStat()	creates a filtered statistical series object.
*
* A filtered statistical series object is a series where each sample is a
* weighted sum of the corresponding samples in another series.  The filter
* function is assumed to be constant, so the new series is, in effect, a
* convolution of the old series with a fixed filter.
*
* Mark Showalter, PDS Ring-Moon Systems Node, March 1998
*******************************************************************************/
#include <stdio.h>
#include "profile.h"
#include "fortran.h"

/********************
* Type definitions *
********************/

typedef struct ZPRO_FILTERED_STRUCT {
PRO_OBJECT	*stat;
RL_INT4	df1, df2, fsize;
RL_FLT8	*filter;
} ZPRO_FILTERED;

/********************************
* Internal function prototypes *
********************************/

static void    ZPro_FilteredKRange RL_PROTO((RL_INT4 knew,
RL_INT4 *k1, RL_INT4 *k2,
RL_VOID *params));
static RL_FLT8 ZPro_FilteredWeight RL_PROTO((RL_INT4 knew, RL_INT4 k,
RL_VOID *params));
static void    ZPro_PrintFiltered  RL_PROTO((RL_VOID *params));

/*
********************************************************************************
* EXPORTED USER ROUTINES
********************************************************************************
*\$ Component_name:
*	Pro_FilteredStat (filtered.c)
*\$ Abstract:
*	This routine creates and returns a filtered statistical series, where
*	each new sample is a weighted sum of samples in another series.
*\$ Keywords:
*	PROFILE, SERIES
*	C, PUBLIC, SUBROUTINE
*\$ Declarations:
*	PRO_OBJECT	*Pro_FilteredStat(stat, filter, df1, df2)
*	PRO_OBJECT	*stat;
*	RL_FLT8		*filter;
*	RL_INT4		df1, df2;
*\$ Inputs:
*	stat		pointer to the statistical series to resample.
*	filter[] 	array of filter coefficients.  The array should be of
*			size |df2-df1|+1.
*	df1, df2	index offset of each array endpoint from filter center.
*			Typically, df1<0 and df2>0.
*\$ Outputs:
*	none
*\$ Returns:
*	pointer to a new statistical series, or NULL on non-fatal error.
*\$ Detailed_description:
*	This routine creates and returns a filtered statistical series, where
*	each new sample is a weighted sum of samples in another series.
*
*	The new series is, in effect, a convolution of the old series with a
*	fixed filter.  The filter coefficients are always adjusted so that the
*	sum of the weights on each new sample is unity.
*\$ External_references:
*	Profile toolkit
*\$ Side_effects:
*	Memory is allocated.  The new object links to the original statistical
*	series.
*\$ Examples:
*	Suppose stat is a statistical series.  Then this code creates a new
*	statistical series in which each sample is an average of the three
*	nearest samples in the original series.
*
*	PRO_OBJECT	*stat3;
*	RL_FLT8		filter[3] = {1., 1., 1.};
*
*	stat3 = Pro_FilteredStat(stat, filter, -1, 1);
*\$ Error_handling:
*	Profile library error handling is in effect.
*
*	Conditions raised:
*	PRO_CLASS_ERROR		if stat is NULL or is not a stat series.
*	RL_MEMORY_ERROR		on memory allocation failure.
*\$ Limitations:
*	none
*\$ Author_and_institution:
*	Mark R. Showalter
*	NASA/Ames Research Center
*\$ Version_and_date:
*	1.0: March 1998
*\$ Change_history:
*	none
*******************************************************************************/

PRO_OBJECT	*Pro_FilteredStat(stat, filter, df1, df2)
PRO_OBJECT	*stat;
RL_FLT8		*filter;
RL_INT4		df1, df2;
{
ZPRO_FILTERED	*params;
PRO_OBJECT	*new;
RL_INT4		i, k1, k2;
RL_FLT8		x1, dx;

/* Verify object type */
if (XPro_StatPtr(stat) == NULL) return NULL;

/* Create new filtered statistical series structure */
params = (ZPRO_FILTERED *) XRL_Malloc(sizeof(ZPRO_FILTERED));
if (params == NULL) return NULL;

params->fsize = df2 - df1 + 1;
if (df2 < df1) params->fsize = df1 - df2 + 1;

params->filter = (RL_FLT8 *) XRL_Malloc(params->fsize * sizeof(RL_FLT8));
if (params->filter == NULL) {
XRL_Free((RL_VOID *) params);
return NULL;
}

/* Fill structure */
params->stat = stat;

params->df1 = df1;
params->df2 = df2;
for (i=0; i < params->fsize; i++)
params->filter[i] = filter[i];

/* Reverse coefficients if necessary */
if (df2 < df1) {
params->df1 = df2;
params->df2 = df1;
for (i=0; i < params->fsize; i++)
params->filter[i] = filter[params->fsize - i];
}

/* Create new weighted statistical series object */
Pro_SeriesSampling(stat, &x1, NULL, &dx);
Pro_SeriesIndices(stat, &k1, &k2);

new = Pro_WeightedStat(stat, k1, k2, x1, dx,
ZPro_FilteredKRange, ZPro_FilteredWeight,
params, sizeof(ZPRO_FILTERED), TRUE,
ZPro_PrintFiltered);
XRL_Free((RL_VOID *) params);

/* NOTE: It is not necessary to enslave old object because Pro_WeightedStat()
does that for us!
if (new != NULL && XPro_EnslaveObject(new,stat)) {
Pro_FreeObject(new);
return NULL;
}
*/

/* Transfer coordinate names to new object */
Pro_RenameObject(new, 1, Pro_ObjectName(stat,1));
Pro_RenameObject(new, 2, Pro_ObjectName(stat,2));

return new;
}

/*
********************************************************************************
* INTERNAL ROUTINES
********************************************************************************
* ZPro_FilteredKRange(knew, k1, k2, params)
*
* This internal routine returns the range of nonzero weights in the old series
* corresponding to a given index in the new series.
*
* Inputs:
*	knew		new index.
*	params		pointer to the ZPRO_FILTERED structure.
*
* Outputs:
*	k1, k2		range of indices in old series for which weight is
*			nonzero.
*******************************************************************************/

static void ZPro_FilteredKRange(knew, k1, k2, params)
RL_INT4		knew, *k1, *k2;
RL_VOID		*params;
{
ZPRO_FILTERED	*filtered;

filtered = (ZPRO_FILTERED *) params;

*k1 = knew + filtered->df1;
*k2 = knew + filtered->df2;
}

/*
********************************************************************************
* ZPro_FilteredWeight(knew, k, params)
*
* This internal routine returns the (un-normalized) weight of an old sample in
* the new series.
*
* Inputs:
*	knew		index in new series.
*	k		index in old series.
*	params		pointer to the ZPRO_FILTERED structure.
*
* Return:		weight value.
*******************************************************************************/

static RL_FLT8	ZPro_FilteredWeight(knew, k, params)
RL_INT4		knew, k;
RL_VOID		*params;
{
ZPRO_FILTERED	*filtered;
RL_INT4		index;

filtered = (ZPRO_FILTERED *) params;

index = k - knew - filtered->df1;
if (index < 0 || index >= filtered->fsize) return 0.;

return filtered->filter[index];
}

/*
********************************************************************************
* ZPro_PrintFiltered(params)
*
* This internal routine prints information about a filtered stat series.
*
* Inputs:
*	params		pointer to the ZPRO_FILTERED data structure.
*******************************************************************************/

static void	ZPro_PrintFiltered(params)
RL_VOID		*params;
{
ZPRO_FILTERED	*filtered;
RL_INT4		i;

filtered = (ZPRO_FILTERED *) params;

/* Print filtered statistical series info... */
printf("\nFiltered statistical series parameters...\n");

for (i = filtered->df1; i <= filtered->df2; i++) {
printf("filter[%d] = %#g\n", i, filtered->filter[i - filtered->df1]);
}
}

/*
********************************************************************************
* FORTRAN INTERFACE ROUTINES
********************************************************************************
*\$ Component_name:
*	FPro_FilteredStat (filtered.c)
*\$ Abstract:
*	This routine creates and returns a filtered statistical series, where
*	each new sample is a weighted sum of samples in another series.
*\$ Keywords:
*	PROFILE, SERIES
*	FORTRAN, PUBLIC, SUBROUTINE
*\$ Declarations:
*	integer*4 function FPro_FilteredStat(stat, filter, df1, df2)
*	integer*4	stat, df1, df2
*	real*8		filter(*)
*\$ Inputs:
*	stat		FORTRAN pointer to the statistical series to resample.
*	filter		array of filter coefficients.  The array should be of
*			size |df2-df1|+1.
*	df1, df2	index offset of each array endpoint from filter center.
*			Typically, df1<0 and df2>0.
*\$ Outputs:
*	none
*\$ Returns:
*	FORTRAN pointer to a new statistical series, or 0 on non-fatal error.
*\$ Detailed_description:
*	This routine creates and returns a filtered statistical series, where
*	each new sample is a weighted sum of samples in another series.
*
*	The new series is, in effect, a convolution of the old series with a
*	fixed filter.  The filter coefficients are always adjusted so that the
*	sum of the weights on each new sample is unity.
*\$ External_references:
*	Profile toolkit
*\$ Side_effects:
*	Memory is allocated.  The new object links to the original statistical
*	series.
*\$ Examples:
*	Suppose stat is a statistical series.  Then this code creates a new
*	statistical series in which each sample is an average of the three
*	nearest samples in the original series.
*
*	integer*4	stat3
*	real*8		filter(3)/ 1.d0, 1.d0, 1.d0 /
*
*	stat3 = FPro_FilteredStat(stat, filter, -1, 1)
*\$ Error_handling:
*	Profile library error handling is in effect.
*
*	Conditions raised:
*	PRO_CLASS_ERROR		if stat is NULL or is not a stat series.
*	RL_MEMORY_ERROR		on memory allocation failure.
*	FORTRAN_POINTER_ERROR	if stat is not a valid FORTRAN object pointer.
*\$ Limitations:
*	none
*\$ Author_and_institution:
*	Mark R. Showalter
*	NASA/Ames Research Center
*\$ Version_and_date:
*	1.0: March 1998
*\$ Change_history:
*	none
*******************************************************************************/

RL_INT4	FORTRAN_NAME(fpro_filteredstat) (stat, filter, df1, df2)
RL_INT4	*stat, *df1, *df2;
RL_FLT8	*filter;
{
RL_VOID *ptr;
RL_INT4 index;

/* Look up statistical series pointer */
ptr = FORT_GetPointer(*stat);
if (ptr == NULL) return 0;

/* Create filtered statistical series */
ptr = (RL_VOID *) Pro_FilteredStat((PRO_OBJECT *) ptr, filter, *df1, *df2);
if (ptr == NULL) return 0;

/* Return new pointer */