Make your own free website on Tripod.com
/*******************************/
/*			       */
/*			       */
/*			       */
/*	Listing 13.1	       */
/*			       */
/*	gridFreq()	       */
/*			       */
/*			       */
/*******************************/
#include "globDefs.h"

real gridFreq( real gridParam[], int gI)
{
real work;
static real incP, incS, freqP, freqS;
static int r, gridDensity, mP, mS, gP;

if (gridParam[0] == 1.0)
{
	gridParam[0] = 0.0;

	freqP = gridParam[1];
	freqS = gridParam[2];
	r = gridParam[3];
	gridDensity = gridParam[4];
	work = (0.5 + freqP - freqS)/r;
	mP = floor(0.5 + freqP/work);
	gridParam[5] = mP;
	gP = mP * gridDensity;
	gridParam[7] = gP;
	mS = r + 1 - mP;
	gridParam[6] = mS;
	incP = freqP / gP;
	incS = (0.5 - freqS) / ((mS - 1) * gridDensity);
}
else
{
	work = (gI <= gP) ? (gI*incP) : (freqS + (gI - (gP + 1)) * incS);
}
return(work);
}



/*************************/
/*	Listing 13.2	 */
/*			 */
/*			 */
/*			 */
/* 	deslpfResp()     */
/*			 */
/*			 */
/*			 */
/*************************/




real desLpfResp(real freqP, real freq)
{
real result;
result = 0.0;
if(freq <= freqP)
	result = 1.0;
return (result);
}


/****************/
/*		*/
/*		*/
/* Listing 13.3	*/
/* weightLp()	*/
/*		*/
/*		*/
/*		*/
/****************/

real weightLp(real kk, real freqP, real freq)
{
real result;
result=1.0;
if(freq<=freqP)
	result = 1.0/kk;
return (result);
}



/***********************/
/*		       */
/* Listing 13.5	       */
/*		       */
/*		       */
/* ComputeRemezA()     */
/*		       */
/*		       */
/*		       */
/***********************/


real computeRemezA(real gridParam[], int gridMax, int r, real kk, real freqP, int iFF[], int initFlag, real contFreq)
{
static int i, j, k, sign;
static real freq, denom, numer, alpha, delta;
static real absDelta, xCont, term;
static real x[50], beta[50], gamma[50];
real aa;

if (initFlag)
{
	for(j=0;j<=r;j++)
	{
		freq=gridFreq(gridParam,iFF[j]);
		x[j] = cos(TWO_PI * freq);
	}

/* compute delta */

denom=0.0;
numer=0.0;
sign=-1;
for(k=0;k<=r;k++)
{
	sign=-sign;
	alpha=1.0;
	for(i=0;i<=(r-1);i++)
	{
		if(i==k)
			continue;
		alpha = alpha / (x[k] - x[i]);
	}

beta[k] = alpha;
if(k != r)
	alpha = alpha/(x[k]-x[r]);

freq = gridFreq(gridParam,iFF[k]);
numer = numer + alpha * desLpfResp(freqP,freq);
denom = denom + sign*(alpha/weightLp(kk,freqP,freq));
}
delta = numer/denom;
absDelta = fabs(delta);

sign = -1;
for (k=0;k<=r-1; k++)
	{
	sign = -sign;
	freq = gridFreq(gridParam, iFF[k]);
	gamma[k] = desLpfResp(freqP, freq) - sign * delta / weightLp(kk, freqP, freq);
	}
}
else
{
	xCont = cos(TWO_PI * contFreq);
	numer = 0.0;
	denom = 0.0;
	for(k=0;k<r;k++)
	{
		term = xCont - x[k];
		if(fabs(term)<1.0e-7)
		{
			aa = gamma[k];
			goto done;
		}
		else
		{
			term = beta[k]/(xCont - x[k]);
			denom += term;
			numer += gamma[k]*term;
		}
	}
	aa = numer/denom;
}
done:
return(aa);
}


/************************/
/*			*/
/*	Listing 13.4	*/
/*			*/
/*			*/
/*	remezError()	*/
/*			*/
/*			*/
/*			*/
/************************/


void remezError(real gridParam[], int gridMax, int r, real kk, real freqP, int iFF[], real ee[])
{
	int j;
	real freq,aa;
	aa = computeRemezA( gridParam, gridMax, r, kk, freqP, iFF, 1, 0.0);
	for (j=0; j<=gridMax; j++)
	{
		freq = gridFreq(gridParam,j);
		aa = computeRemezA( gridParam, gridMax, r, kk, freqP, iFF, 0, freq);
		ee[j] = weightLp(kk, freqP, freq) * (desLpfResp(freqP, freq) - aa);
	}
	return;
}

/***********************************/
/*				   */
/*				   */
/*				   */
/*				   */
/*     listing 13.6		   */
/*				   */
/*     remezsearch()		   */
/*				   */
/*				   */
/*				   */
/*				   */
/*				   */
/***********************************/


void remezSearch(real ee[] ,real absDelta, int gP, int iFF[], int gridMax, int r, real gridParam[])
{
	int i,j,k,extras, indexOfSmallest;
	real smallestVal;

	k=0;

/* test for extremum at f=0 */

	if (((ee[0]>0.0) && (ee[0]>ee[1]) && (fabs(ee[0])>=absDelta)) || ((ee[0]<0.0) && (ee[0]=absDelta)))
	{
		iFF[k]=0;
		k++;
	}

/* search for extrema in passband */

	for (j=1; j=ee[j-1]) && (ee[j]>ee[j+1]) && (ee[j]>0.0)) || ((ee[j]<=ee[j-1]) && (ee[j]=ee[j-1]) && (ee[j]>ee[j+1]) && (ee[j]>0.0)) || ((ee[j]<=ee[j-1]) && (ee[j]0.0) && (ee[j]>ee[j-1]) && (fabs(ee[j])>=absDelta)) || ((ee[j]<0.0) && (ee[j]=absDelta)))
	{
		iFF[k]=gridMax;
		k++;
	}

/*------------------------------------------------------------*/

/* find and remove superfluous extermal frequencies  */

	if (k>r+1)
	{
		extras = k - (r+1);
		for (i=1; i<=extras; i++)
		{
			smallestVal = fabs(ee[iFF[0]]);
			indexOfSmallest = 0;
			for (j=1; j<k; j++)
			{
				if(fabs(ee[iFF[j]]) >= smallestVal) 
					continue;
				smallestVal = fabs(ee[iFF[j]]);
				indexOfSmallest = j;
			}
			k--;
			for(j=indexOfSmallest; j<k; j++)
				iFF[j]=iFF[j+1];
		}
	}
	return;
}

/********************************/
/*				*/	
/*				*/	
/*				*/	
/* Listing 13.7			*/	
/*				*/	
/* remezstop()			*/
/*				*/	
/********************************/


int remezStop( int iFF[], int r)
{
	static int oldiFF[50];
	int j,result;

	result = 1;
	for (j=0; j<=r; j++)
	{
		if(iFF[j] != oldiFF[j]) 
			result = 0;
		oldiFF[j] = iFF[j];
	}
	return(result);
}

/********************************/
/*				*/
/*				*/
/*				*/
/*	Listing 13.8		*/
/*				*/
/*				*/
/*      remezStop2()		*/
/*				*/
/*				*/
/*				*/
/********************************/
	
int remezStop2(real ee[], int iFF[], int r)
{
	real biggestVal, smallestVal, qq;
	int j,result;
	result=0;
	biggestVal = fabs(ee[iFF[0]]);
	smallestVal = fabs(ee[iFF[0]]);
	for (j=1; j<=r; j++)
	{
		if (fabs(ee[iFF[j]]) < smallestVal)
			smallestVal = fabs(ee[iFF[j]]);
		if (fabs(ee[iFF[j]]) > biggestVal)
			biggestVal = fabs(ee[iFF[j]]);
	}
	qq = (biggestVal - smallestVal)/biggestVal;
	if(qq<0.01)
		result=1;
	return(result);
}

/****************************************/
/*					*/
/*					*/
/*					*/
/*					*/
/*	Listing 13.9			*/
/*					*/
/*					*/
/*	remezFinish()			*/
/*					*/
/*					*/
/****************************************/

void remezFinish(real extFreq[], int nn, int r, real freqP, real kk, real aa[], real h[])
{
	int k,n,gridMax,iFF[1];
	real freq,sum;
	static real gridParam[1];

	for(k=0; k<r; k++)
	{
		freq = (real) k / (real) nn;
		aa[k] = computeRemezA( gridParam, gridMax, r, kk, freqP, iFF, 0, freq);
	}
	fsDesign(nn, 1, aa, h);
	return;
}

/****************************************/ 
/*					*/
/*					*/
/*		Listing 13.10		*/
/*					*/
/*		remez()			*/
/*					*/
/*					*/
/****************************************/ 


void remez( int nn, int r, int gridDensity, real kk, real freqP, real freqS, real extFreq[], real h[])
{
	int m,gridMax,j,mP,gP,mS;
	real absDelta,freq;
	static real gridParam[10];
	static int iFF[50];
	static real ee[1024];

	/*-------------------------------*/

	/* set up frequency grid         */

	gridParam[0] = 1.0;
	gridParam[1] = freqP;
	gridParam[2] = freqS;
	gridParam[3] = r;
	gridParam[4] = gridDensity;
	freq = gridFreq(gridParam,0);
	mP = gridParam[5];
	mS = gridParam[6];
	gP = gridParam[7];
	freqP = freqP + (freqP/(2.0*gP));
	gridMax = 1 + gridDensity*(mP+mS-1);
	/*--------------------------------*/
	
	/*  make initial guess of extremal frequencies */

	for (j=0; j<mP; j++)
		iFF[j] = (j+1)*gridDensity;
	for (j=0; j<mS; j++)
		iFF[j+mP] = gP + 1 + j * gridDensity;

	/*---------------------------------*/

	/*  find optimal locations for extremal frequencies */

	for (m=1; m<=20; m++)
	{
		remezError( gridParam, gridMax, r, kk, freqP, iFF, ee);
		remezSearch(ee, absDelta, gP, iFF, gridMax, r, gridParam);
		remezStop2(ee,iFF,r);
		if (remezStop(iFF,r))
			break;
	}
	
	
	for (j=0; j<=r; j++)
	{
		extFreq[j] = gridFreq(gridParam, iFF[j]);
	}
	remezFinish( extFreq, nn, r, freqP, kk, ee, h);

	return;
}

/********************************/
/*				*/		
/*				*/		
/*				*/		
/*	Listing 12.1		*/
/*				*/
/*	fsdesign()		*/
/*				*/
/*				*/
/*				*/
/********************************/


int fsDesign( int N, int firType, real A[], real h[])
{
	int n,k,status;
	real x, M;
	M=(N-1.0)/2.0;
	status=0;
	switch(firType)
	{
		case 1:
			if(N%2)
			{
				for(n=0;n<N;n++)
				{
					h[n]=A[0];
					x=TWO_PI*(n-M)/N;
					for(k=1;k<=M;k++)
					{
						h[n]=h[n]+2.0*A[k]*cos(x*k);
					}
					h[n]=h[n]/N;
				}
			}
			else
				{status = 1;}
			break;
/**************************************************************************/

		case 2:
			if(N%2)
				{status=2;}
			else {
				for(n=0;n<N;n++){
					h[n]=A[0];
					x=TWO_PI*(n-M)/N;
				for(k=1;k<=(N/2-1);k++){
					h[n]=h[n]+2.0*A[k]*cos(x*k);
					}
				h[n]=h[n]/N;
				}
				}
				break;
/**************************************************************************/

	case 3:
		if(N%2) {
			for (n=0;n<N;n++) {
				h[n]=0;
				x=TWO_PI*(M-n)/N;
				for (k=1;k<=M;k++) {
					h[n]=h[n]+2.0*A[k]*sin(x*k);
					}
				h[n]=h[n]/N;
				}
			}
		else 
			{status = 3;}
		break;
/*************************************************************************/
	case 4:
		if(N%2)
			{status = 4;}
		else {
			for(n=0;n<N;n++) {
				h[n]=A[N/2]*sin(PI*(M-n));
				x=TWO_PI*(n-M)/N;
				for(k=1;k<=(N/2-1);k++){
					h[n]=h[n]+2.0*A[k]*sin(x*k);
					}
				h[n]=h[n]/N;
				}
			}
			break;
	}
	return(status);
}
/*
main()
{
int i,n,r,l;
real k,fp,fs;
real e_freq[100],h_samp[100];
printf("Enter the number of coefficients (n)\n");
scanf("%d",&n);
printf("Enter the Number of approximating functions (r)\n");
scanf("%d",&r);
printf("Enter the grid density (L)\n");
scanf("%d",&l);
printf("Enter the ripple ratio \"delta1/delta2\" (k)\n");
scanf("%lf",&k);
printf("Enter the Passband Frequency\n");
scanf("%lf",&fp);
printf("Enter the Stopband Frequency\n");
scanf("%lf",&fs);
remez(n,r,l,k,fp,fs,e_freq,h_samp);
printf("COEFFICIENTS\n");
printf("------------\n\n");
for (i=0;i<n;i++)
	printf("h[%d] = %lf\n",i,h_samp[i]);
printf("\n\n");
printf("EXTREMAL FREQUENCIES\n");
printf("------------\n\n");
for (i=0;i<=r;i++)
	printf("f[%d] = %lf\n",i,e_freq[i]);
}
*/