/*
 *	*************************************************************************
 *	*									*
 *	*	This confidential and proprietary software may be used only	*
 *	*  as authorized by a licensing agreement from the Alta Group of	*
 *	*  Cadence Design Systems, Inc.  In the event of publication, the	*
 *	*  following notice is applicable:					*
 *	*									*
 *	*    (c) COPYRIGHT 1995 ALTA GROUP OF CADENCE DESIGN SYSTEMS, INC.	*
 *	*			ALL RIGHTS RESERVED				*
 *	*									*
 *	*	The entire notice above must be reproduced on all authorized	*
 *	*  copies.								*
 *	*									*
 *	*************************************************************************
 *
 */
/*
 * FILE:  parts/spb/cgs/generic/ovector/dvgenw.c
 * DATE:  Wed Feb 19, 1992
 * RELATED FILES:  
 * AUTHOR:  John Lundell
 * DESCRIPTION:  'c' DvGenWindow function
 * NOTES/WARNINGS:  
 * REVISION HISTORY:
 *	Release		Who	Date	Comments
 */
#include "cgs.h"

/*
 *  Functions
 */
static double i0();

/*
 * FUNCTION:  DvGenWindow
 * DESCRIPTION:
 *	DvGenWindow(i_type, d_param, sp_dst)
 *
 *  Generate a window function of type 'i_type' using
 *  parameter 'd_parm', placing the result in the
 *  ovector 'sp_dst'.  Type should be one of the following:
 *	OvBARTLETT	OvBARTLETT_FFT
 *	OvBLACKMAN	OvBLACKMAN_FFT
 *	OvHAMMING 	OvHAMMING_FFT
 *	OvHANNING 	OvHANNING_FFT
 *	OvKAISER  	OvKAISER_FFT
 *	OvCHEBYSHEV	OvCHEBYSHEV_FFT
 *
 * RETURN VALUE: sp_dst
 * NOTES/WARNINGS: 
 * REVISION HISTORY:
 */
Ovector DvGenWindow(i_type, d_param, sp_dst)
int		i_type;
double		d_param;
Ovector		sp_dst;
{
	register int		i,j,st;
	register Data		frw, bac;
	double			s1, s2, s3, wv;
	int			length = OvGetLength(sp_dst);
	int			half;
/*
 *  Get pointer values
 */
	frw = OvGetVirtStart(sp_dst);
	st = sp_dst->stepsize;
	bac = frw + (length - 1) * st;
	if (i_type & OvFFT_WINDOW) length++;

	half = (length + 1) >> 1;
/*
 *  Loop over vector for the given table type.
 */
	switch(i_type) {
	      case OvBARTLETT:
		*(double *) bac = 0.0;
		bac -= st;
	      case OvBARTLETT_FFT:
		*(double *) frw = 0.0;
		frw += st;
		
		s1 = 2.0 / (length - 1.0);
       		for (i = 1; i < half; i++){
			*(double *) frw = *(double *) bac = s1 * i;
			frw += st;
			bac -= st;
		}
		break;

	      case OvBLACKMAN:
		*(double *) bac = 0.0;
		bac -= st;
	      case OvBLACKMAN_FFT:
		*(double *) frw = 0.0;
		frw += st;
		
		s1 = PI2 / (length - 1.0);
		s2 = 2.0 * PI2 / (length - 1.0);
       		for (i = 1; i < half; i++){
			*(double *) frw = *(double *) bac = 0.42 - 0.5*cos(s1 * i) + 0.08*cos(s2 * i);
			frw += st;
			bac -= st;
		}
		break;
		
	      case OvHAMMING:
		*(double *) bac = 0.08;
		bac -= st;
	      case OvHAMMING_FFT:
		*(double *) frw = 0.08;
		frw += st;
		
		s1 = PI2 / (length - 1.0);
		for (i = 1; i < half; i++) {
			*(double *) frw = *(double *) bac = 0.54 - 0.46*cos(s1 * i);
			frw += st;
			bac -= st;
		}
		break;
		
	      case OvHANNING:
		*(double *) bac = 0.0;
		bac -= st;
	      case OvHANNING_FFT:
		*(double *) frw = 0.0;
		frw += st;
		
		s1 = PI2 / (length - 1.0);
       		for(i=1; i < half; i++){
			*(double *) frw = *(double *) bac = 0.5 - 0.5*cos(s1 * i);
			frw += st;
			bac -= st;
	        }
		break;
		
	      case OvKAISER:
		s1 = 1.0 / i0(d_param);
		*(double *) bac = s1;
		bac -= st;
	      case OvKAISER_FFT:
		s1 = 1.0 / i0(d_param);
		*(double *) frw = s1;
		frw += st;

		s2 = 2.0 / (length - 1.0);
		for (i = 1; i < half; i++) {
			s3 = s2 * i - 1;
			*(double *) frw = *(double *) bac = s1 * i0(d_param * sqrt(1 - s3*s3));
			frw += st;
			bac -= st;
		}
		break;

	      case OvCHEBYSHEV:
		*(double *) bac = 1.0;
		bac -= st;
	      case OvCHEBYSHEV_FFT:
		*(double *) frw = 1.0;
		frw += st;
/*
 *  Calculate the factor beta
 */
		d_param = pow(10.0, -d_param / 20.0);
		d_param = log(d_param + sqrt(d_param * d_param - 1));
		d_param = exp(2.0 * d_param / (length - 1.0));
		d_param = (d_param - 1.0) * (d_param - 1.0) / ((d_param + 1.0) * (d_param + 1.0));
/*
 *  Factor which indicates when we are getting close enough
 */
		s3 = d_param * 1.0e-10;
/*
 *  Generate the window
 */
		for (i = 1; i < half; i++) {
			s1 = d_param;
			wv = 0.0;
			s2 = length - i - 1.0;
			for (j = 1; j <= i; j++) {
				wv += s1;
				s1 *= d_param * (i-j) * (s2 - j) / (j * (j + 1.0));
				if (s1 < s3) break; 
			}
			*(double *) frw  = *(double *) bac = (length - 1.0) * wv;
			frw += st;
			bac -= st;
		}
/*
 *  Scale so that the middle point is one
 */
		frw -= st;
		s1 = 1.0 / *(double *) frw;
		Dvscale(s1, sp_dst, sp_dst);
		break;
	      default:
		printf("Invalid table type to DvGenWindow.\n");

	}
	
	return sp_dst;
}


/*---------------------------------------------------------------
 * FUNCTION:  i0
 * DESCRIPTION:
 *	val = i0(d_x)
 *
 *  Simple approximation of a I0 bessel function. 
 * RETURN VALUE:  bessel function value
 * NOTES/WARNINGS: 
 * REVISION HISTORY:
 *	Release		Who	Date	Comments
 */
static double i0(d_x)
double	d_x;
{
	register double	sum = 0.0;
	register double val = 1.0;
	register double	stp;
	register int	i;

	d_x *= 0.5;
	stp = (1.0 + d_x) * 1.0e-10;
	
	for (i = 1; i < 30; i++) {
		sum += val * val;
		val *= d_x / ((double ) i);
		if (val < stp) break;
	}

	return sum;
}