remez.c

大名鼎鼎的CE下播放软件,TCPPMP的源代码!!!2410下可以流畅的解QVGA的H264,MPEG4等格式.

	for (i = 1; i < gridsize - 1; i++) {
		if (((E[i] >= E[i - 1]) && (E[i] > E[i + 1]) && (E[i] > 0.0)) ||
		   ((E[i] <= E[i - 1]) && (E[i] < E[i + 1]) && (E[i] < 0.0)))
			foundExt[k++] = i;
	}

	/* Check for extremum at 0.5 */
	j = gridsize - 1;
	if (((E[j] > 0.0) && (E[j] > E[j - 1])) ||
	   ((E[j] < 0.0) && (E[j] < E[j - 1])))
		foundExt[k++] = j;

	/* Remove extra extremals */
	extra = k - (r + 1);

	while (extra > 0) {
		if (E[foundExt[0]] > 0.0)
			up = 1;                /* first one is a maxima */
		else
			up = 0;                /* first one is a minima */

		l = 0;
		alt = 1;
		for (j = 1; j < k; j++) {
			if (fabs(E[foundExt[j]]) < fabs(E[foundExt[l]]))
				l = j;              /* new smallest error. */
			if ((up) && (E[foundExt[j]] < 0.0))
				up = 0;             /* switch to a minima */
			else if ((!up) && (E[foundExt[j]] > 0.0))
				up = 1;             /* switch to a maxima */
			else {
				alt = 0;
				break;              /* Ooops, found two non-alternating */
			}                       /* extrema.  Delete smallest of them */
		}  /* if the loop finishes, all extrema are alternating */

		/* If there's only one extremal and all are alternating,
		 * delete the smallest of the first/last extremals. */
		if ((alt) && (extra == 1)) {
			if (fabs(E[foundExt[k - 1]]) < fabs(E[foundExt[0]]))
				l = foundExt[k - 1];   /* Delete last extremal */
			else
				l = foundExt[0];       /* Delete first extremal */
		}

		/* Loop that does the deletion */
		for (j = l; j < k; j++) {
			foundExt[j] = foundExt[j+1];
		}
		k--;
		extra--;
	}

	/* Copy found extremals to Ext[] */
	for (i = 0; i <= r; i++) {
		Ext[i] = foundExt[i];
	}

	free(foundExt);
}


/*********************
 * FreqSample
 *============
 * Simple frequency sampling algorithm to determine the impulse
 * response h[] from A's found in ComputeA
 *
 *
 * INPUT:
 * ------
 * int      N        - Number of filter coefficients
 * double   A[]      - Sample points of desired response [N/2]
 * int      symm     - Symmetry of desired filter
 *
 * OUTPUT:
 * -------
 * double h[] - Impulse Response of final filter [N]
 *********************/
static void FreqSample(int N, const double A[], double h[], int symm)
{
	int n, k;
	double x, val, M;

	M = (N - 1.0) / 2.0;
	if (symm == POSITIVE) {
		if (N % 2) {
			for (n = 0; n < N; n++) {
				val = A[0];
				x = Pi2 * (n - M) / N;
				for (k = 1; k <= M; k++)
					val += 2.0 * A[k] * cos(x * k);
				h[n] = val / N;
			}
		}
		else {
			for (n = 0; n < N; n++) {
				val = A[0];
				x = Pi2 * (n - M)/N;
				for (k = 1; k <= (N / 2 - 1); k++)
					val += 2.0 * A[k] * cos(x * k);
				h[n] = val / N;
			}
		}
	}
	else {
		if (N % 2) {
			for (n = 0; n < N; n++) {
				val = 0;
				x = Pi2 * (n - M) / N;
				for (k = 1; k <= M; k++)
					val += 2.0 * A[k] * sin(x * k);
				h[n] = val / N;
			}
		}
		else {
			for (n = 0; n < N; n++) {
				val = A[N / 2] * sin(Pi * (n - M));
				x = Pi2 * (n - M) / N;
				for (k = 1; k <= (N / 2 - 1); k++)
					val += 2.0 * A[k] * sin(x * k);
				h[n] = val / N;
			}
		}
	}
}

/*******************
 * isDone
 *========
 * Checks to see if the error function is small enough to consider
 * the result to have converged.
 *
 * INPUT:
 * ------
 * int    r     - 1/2 the number of filter coeffiecients
 * int    Ext[] - Indexes to extremal frequencies [r+1]
 * double E[]   - Error function on the dense grid [gridsize]
 *
 * OUTPUT:
 * -------
 * Returns 1 if the result converged
 * Returns 0 if the result has not converged
 ********************/

static int isDone(int r, const int Ext[], const double E[])
{
	int i;
	double min, max, current;

	min = max = fabs(E[Ext[0]]);
	for (i = 1; i <= r; i++) {
		current = fabs(E[Ext[i]]);
		if (current < min)
			min = current;
		if (current > max)
			max = current;
	}
	if (((max - min) / max) < 0.0001)
		return 1;
	return 0;
}

/********************
 * remez
 *=======
 * Calculates the optimal (in the Chebyshev/minimax sense)
 * FIR filter impulse response given a set of band edges,
 * the desired reponse on those bands, and the weight given to
 * the error in those bands.
 *
 * INPUT:
 * ------
 * int     numtaps     - Number of filter coefficients
 * int     numband     - Number of bands in filter specification
 * double  bands[]     - User-specified band edges [2 * numband]
 * double  des[]       - User-specified band responses [numband]
 * double  weight[]    - User-specified error weights [numband]
 * int     type        - Type of filter
 *
 * OUTPUT:
 * -------
 * double h[]      - Impulse response of final filter [numtaps]
 ********************/

void remez(double h[], int numtaps, int numband, double bands[],
           const double des[], const double weight[], int type)
{
	double *Grid, *W, *D, *E;
	int    i, iter, gridsize, r, *Ext;
	double *taps, c;
	double *x, *y, *ad;
	int    symmetry;

	if (type == BANDPASS)
		symmetry = POSITIVE;
	else
		symmetry = NEGATIVE;

	r = numtaps / 2;                  /* number of extrema */
	if ((numtaps % 2) && (symmetry == POSITIVE))
		r++;

	/* Predict dense grid size in advance for memory allocation
	 *   .5 is so we round up, not truncate */
	gridsize = 0;
	for (i = 0; i < numband; i++) {
		gridsize += (int) (2 * r * GRIDDENSITY *
		                   (bands[2 * i + 1] - bands[2 * i]) + .5);
	}
	if (symmetry == NEGATIVE) {
		gridsize--;
	}

	/* Dynamically allocate memory for arrays with proper sizes */
	Grid = (double *) Util_malloc(gridsize * sizeof(double));
	D = (double *) Util_malloc(gridsize * sizeof(double));
	W = (double *) Util_malloc(gridsize * sizeof(double));
	E = (double *) Util_malloc(gridsize * sizeof(double));
	Ext = (int *) Util_malloc((r + 1) * sizeof(int));
	taps = (double *) Util_malloc((r + 1) * sizeof(double));
	x = (double *) Util_malloc((r + 1) * sizeof(double));
	y = (double *) Util_malloc((r + 1) * sizeof(double));
	ad = (double *) Util_malloc((r + 1) * sizeof(double));

	/* Create dense frequency grid */
	CreateDenseGrid(r, numtaps, numband, bands, des, weight,
	                &gridsize, Grid, D, W, symmetry);
	InitialGuess(r, Ext, gridsize);

	/* For Differentiator: (fix grid) */
	if (type == DIFFERENTIATOR) {
		for (i = 0; i < gridsize; i++) {
			/* D[i] = D[i] * Grid[i]; */
			if (D[i] > 0.0001)
				W[i] = W[i] / Grid[i];
		}
	}

	/* For odd or Negative symmetry filters, alter the
	 * D[] and W[] according to Parks McClellan */
	if (symmetry == POSITIVE) {
		if (numtaps % 2 == 0) {
			for (i = 0; i < gridsize; i++) {
				c = cos(Pi * Grid[i]);
				D[i] /= c;
				W[i] *= c;
			}
		}
	}
	else {
		if (numtaps % 2) {
			for (i = 0; i < gridsize; i++) {
				c = sin(Pi2 * Grid[i]);
				D[i] /= c;
				W[i] *= c;
			}
		}
		else {
			for (i = 0; i < gridsize; i++) {
				c = sin(Pi * Grid[i]);
				D[i] /= c;
				W[i] *= c;
			}
		}
	}

	/* Perform the Remez Exchange algorithm */
	for (iter = 0; iter < MAXITERATIONS; iter++) {
		CalcParms(r, Ext, Grid, D, W, ad, x, y);
		CalcError(r, ad, x, y, gridsize, Grid, D, W, E);
		Search(r, Ext, gridsize, E);
		if (isDone(r, Ext, E))
			break;
	}
#ifndef ASAP
	if (iter == MAXITERATIONS) {
		Aprint("remez(): reached maximum iteration count. Results may be bad.");
	}
#endif

	CalcParms(r, Ext, Grid, D, W, ad, x, y);

	/* Find the 'taps' of the filter for use with Frequency
	 * Sampling.  If odd or Negative symmetry, fix the taps
	 * according to Parks McClellan */
	for (i = 0; i <= numtaps / 2; i++) {
		if (symmetry == POSITIVE) {
			if (numtaps % 2)
				c = 1;
			else
				c = cos(Pi * (double) i / numtaps);
		}
		else {
			if (numtaps % 2)
				c = sin(Pi2 * (double) i / numtaps);
			else
				c = sin(Pi * (double) i / numtaps);
		}
		taps[i] = ComputeA((double) i / numtaps, r, ad, x, y) * c;
	}

	/* Frequency sampling design with calculated taps */
	FreqSample(numtaps, taps, h, symmetry);

	/* Delete allocated memory */
	free(Grid);
	free(W);
	free(D);
	free(E);
	free(Ext);
	free(taps);
	free(x);
	free(y);
	free(ad);
}