bessel.c

ofdm的完整系统模型,包含信道参数,多径模型,doppler频移等都可以自由修改!是您做仿真的有力帮助.c语言运行速度快!

/*
 *	*************************************************************************
 *	*									*
 *	*	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/cgsmisc/bessel.c
 * DATE:  Fri Feb 12, 1993
 * RELATED FILES:  
 * AUTHOR:  John Lundell
 * DESCRIPTION:
 *	Zeroth order Bessel and modified Bessl functions
 *
 * NOTES/WARNINGS:  
 * REVISION HISTORY:
 *	Release		Who	Date	Comments
 */
#include "cgs.h"
/*
 *  Globals & Defines
 */

double PIO4   =  7.85398163397448309616E-1;    /* pi/4 */
double SQ2OPI =  7.9788456080286535587989E-1;  /* sqrt( 2/pi ) */

#define GAM25	3.625609908221908
#define PI	3.141592653589793
#define EPS	1.0e-12

static double PP[7] = {
	7.96936729297347051624E-4,
	8.28352392107440799803E-2,
	1.23953371646414299388E0,
	5.44725003058768775090E0,
	8.74716500199817011941E0,
	5.30324038235394892183E0,
	9.99999999999999997821E-1,
};
static double PQ[7] = {
	9.24408810558863637013E-4,
	8.56288474354474431428E-2,
	1.25352743901058953537E0,
	5.47097740330417105182E0,
	8.76190883237069594232E0,
	5.30605288235394617618E0,
	1.00000000000000000218E0,
};
static double QP[8] = {
	-1.13663838898469149931E-2,
	-1.28252718670509318512E0,
	-1.95539544257735972385E1,
	-9.32060152123768231369E1,
	-1.77681167980488050595E2,
	-1.47077505154951170175E2,
	-5.14105326766599330220E1,
	-6.05014350600728481186E0,
};
static double QQ[7] = {
	6.43178256118178023184E1,
	8.56430025976980587198E2,
	3.88240183605401609683E3,
	7.24046774195652478189E3,
	5.93072701187316984827E3,
	2.06209331660327847417E3,
	2.42005740240291393179E2,
};
static double YP[8] = {
	1.55924367855235737965E4,
	-1.46639295903971606143E7,
	5.43526477051876500413E9,
	-9.82136065717911466409E11,
	8.75906394395366999549E13,
	-3.46628303384729719441E15,
	4.42733268572569800351E16,
	-1.84950800436986690637E16,
};
static double YQ[7] = {
	1.04128353664259848412E3,
	6.26107330137134956842E5,
	2.68919633393814121987E8,
	8.64002487103935000337E10,
	2.02979612750105546709E13,
	3.17157752842975028269E15,
	2.50596256172653059228E17,
};
static double DR1 = 5.78318596294678452118E0;
static double DR2 = 3.04712623436620863991E1;

static double RP[4] = {
	-4.79443220978201773821E9,
	1.95617491946556577543E12,
	-2.49248344360967716204E14,
	9.70862251047306323952E15,
};
static double RQ[8] = {
	4.99563147152651017219E2,
	1.73785401676374683123E5,
	4.84409658339962045305E7,
	1.11855537045356834862E10,
	2.11277520115489217587E12,
	3.10518229857422583814E14,
	3.18121955943204943306E16,
	1.71086294081043136091E18,
};

static double A[] =
{
	-4.41534164647933937950E-18,
	3.33079451882223809783E-17,
	-2.43127984654795469359E-16,
	1.71539128555513303061E-15,
	-1.16853328779934516808E-14,
	7.67618549860493561688E-14,
	-4.85644678311192946090E-13,
	2.95505266312963983461E-12,
	-1.72682629144155570723E-11,
	9.67580903537323691224E-11,
	-5.18979560163526290666E-10,
	2.65982372468238665035E-9,
	-1.30002500998624804212E-8,
	6.04699502254191894932E-8,
	-2.67079385394061173391E-7,
	1.11738753912010371815E-6,
	-4.41673835845875056359E-6,
	1.64484480707288970893E-5,
	-5.75419501008210370398E-5,
	1.88502885095841655729E-4,
	-5.76375574538582365885E-4,
	1.63947561694133579842E-3,
	-4.32430999505057594430E-3,
	1.05464603945949983183E-2,
	-2.37374148058994688156E-2,
	4.93052842396707084878E-2,
	-9.49010970480476444210E-2,
	1.71620901522208775349E-1,
	-3.04682672343198398683E-1,
	6.76795274409476084995E-1,
};
static double B[] =
{
	-7.23318048787475395456E-18,
	-4.83050448594418207126E-18,
	4.46562142029675999901E-17,
	3.46122286769746109310E-17,
	-2.82762398051658348494E-16,
	-3.42548561967721913462E-16,
	1.77256013305652638360E-15,
	3.81168066935262242075E-15,
	-9.55484669882830764870E-15,
	-4.15056934728722208663E-14,
	1.54008621752140982691E-14,
	3.85277838274214270114E-13,
	7.18012445138366623367E-13,
	-1.79417853150680611778E-12,
	-1.32158118404477131188E-11,
	-3.14991652796324136454E-11,
	1.18891471078464383424E-11,
	4.94060238822496958910E-10,
	3.39623202570838634515E-9,
	2.26666899049817806459E-8,
	2.04891858946906374183E-7,
	2.89137052083475648297E-6,
	6.88975834691682398426E-5,
	3.36911647825569408990E-3,
	8.04490411014108831608E-1,
};


/*
 *  Functions
 */
double cgsPolyEval(), cgsPolyEval1();
double cgsChebyEval();
static double	ascend(), hankel();

/*---------------------------------------------------------------
 * FUNCTION:  J0
 * DESCRIPTION:
 *	double J0(d_x)
 *
 *  Zero order Bessel function.
 *
 * RETURN VALUE:  function value
 * NOTES/WARNINGS: 
 * REVISION HISTORY:
 *	Release		Who	Date	Comments
 */
double J0(d_x)
double	d_x;
{
	double	p, q, w, z, xn;

	if (d_x < 0.0) d_x = -d_x;
	
	if (d_x <= 5.0) {
		z = d_x * d_x;
		if (d_x < 1.0e-5) return (1.0 - z/4.0);
		  
		p = (z - DR1) * (z - DR2);
		p = p * cgsPolyEval(z, RP, 3) / cgsPolyEval1(z, RQ, 8);
		return p;
	}
	
	w = 5.0 / d_x;

	q = 25.0 / (d_x * d_x);
	p = cgsPolyEval(q, PP, 6) / cgsPolyEval(q, PQ, 6);
	q = cgsPolyEval(q, QP, 7) / cgsPolyEval1(q, QQ, 7);
	xn = d_x - PIO4;
	p = p * cos(xn) - w * q * sin(xn);

	return p * SQ2OPI / sqrt(d_x);
}

/*---------------------------------------------------------------
 * FUNCTION:  I0
 * DESCRIPTION:  
 *	double IO(d_x)
 *
 *  Zero order modified Bessel function.
 *
 * RETURN VALUE:  function value
 * NOTES/WARNINGS: 
 * REVISION HISTORY:
 *	Release		Who	Date	Comments
 */
double I0(d_x)
double	d_x;
{
	double	y;
	
	if (d_x < 0.0) d_x = -d_x;

	if (d_x <= 8.0) {
		y = 0.5 * d_x - 2.0;
		return (exp(d_x) * cgsChebyEval(y, A, 30));
	}
	
	return (exp(d_x) * cgsChebyEval(32.0/d_x - 2.0, B, 25) / sqrt(d_x));
	
}

/*---------------------------------------------------------------
 * FUNCTION:  cgsPolyEval
 * DESCRIPTION:
 *	double cgsPolyEval(d_x, dp_coefs, i_order)
 *
 *  Evaluates the polynomial of order 'i_order' and
 *  coefficients 'dp_coefs' at the point 'd_x'.
 *
 * RETURN VALUE:  polynomial value
 * NOTES/WARNINGS: 
 * REVISION HISTORY:
 *	Release		Who	Date	Comments
 */
double cgsPolyEval(d_x, dp_coef, i_order)
register double	d_x;
register double	*dp_coef;
register int	i_order;
{
	register double	ans;
	
	ans = *dp_coef++;
	
	do {
		ans = ans * d_x  +  *dp_coef++;
	  } while(i_order-- > 1);
	
	return ans;
}

/*---------------------------------------------------------------
 * FUNCTION:  cgsPolyEval1
 * DESCRIPTION:
 *	double cgsPolyEval1(d_x, dp_coef, i_order)
 *
 *  Evaluates the polynomial of order 'i_order' and
 *  coefficients 'dp_coefs' at the point 'd_x'.  The
 *  highest order coefficient is assumed to be one.
 *
 * RETURN VALUE:  polynomial value
 * NOTES/WARNINGS: 
 * REVISION HISTORY:
 *	Release		Who	Date	Comments
 */
double cgsPolyEval1(d_x, dp_coef, i_order)
register double d_x;
register double *dp_coef;
register int	i_order;
{
	register double ans;
	
	ans = d_x + *dp_coef++;
	
	do {
		ans = ans * d_x  + *dp_coef++;
	} while(i_order-- > 2);
	
	return ans;
}

/*---------------------------------------------------------------
 * FUNCTION:  cgsChebyEval
 * DESCRIPTION:
 *	double cgsChebyEval(d_x, dp_coefs, i_order)
 *
 *  Evaluates the Chebyshev series 'dp_coefs' at the
 *  point 'd_x'.
 *
 * RETURN VALUE:  Series value
 * NOTES/WARNINGS: 
 * REVISION HISTORY:
 *	Release		Who	Date	Comments
 */
double cgsChebyEval(d_x, dp_coef, i_order)
double	d_x;
double	*dp_coef;
int	i_order;
{
	double	b2;
	double	b1 = 0.0;
	double	b0 = *dp_coef++;
	
	do {
		b2 = b1;
		b1 = b0;
		b0 = d_x * b1  -  b2  + *dp_coef++;
	} while (i_order-- > 2);
	
	return  0.5 * (b0-b2);
}

double J25(x)
double	 x;
{
	if (x < 18.0) {
		return ascend(x);
	}
	else {
		return hankel(x);
	}
}

static double ascend(x)
double	x;
{
	double	k = 1.0;
	double	z = -x*x/4.0;
	double	den = 1.25;
	double	v = 1.0;
	double	u = 1.0;
	
	do {
		u *= z / (k * den);
		k += 1.0;
		den += 1.0;
		v += u;
	} while (fabs(u/v) > EPS);

	v *= pow(0.5*x, 0.25) / (0.25*GAM25);
	
	return v;
}

static double hankel(x)
double	x;
{
	double	p = 1.0;
	double	q = 0.0;
	double	k = 1.0;
	double	sgn = 1.0;
	double	v = 1.0;
	double	z = 8*x;
	
	do {
		v *= (0.25 - k*k) / (z * 0.5 * (k+1.0));
		q += sgn * v;
		sgn *= -1.0;
		k += 2.0;
		v *= (0.25 - k*k) / (z * 0.5 * (k+1.0));
		p += sgn * v;
		k += 2.0;
	} while(fabs(v/p) > EPS || fabs(v/q) > EPS);

	z = x - 3.0*PI/8.0;
	
	v = sqrt(2.0/(PI*x)) * (p*cos(z) - q*sin(z));

	return v;
}