equiripple_mband.cpp

// ntaps: the number of filter coefficients
// M: Every Mth coefficient will be zero, this is the "M"
//   in M-band
// fp: Filter's passband cutoff.  Can be zero.
LPFIRLITE	equiripple_halfband_fl(int ntaps, const int M,
		const double fp) {
	// Assume that our filters will always have odd length
	ntaps = (ntaps & 1) ? ntaps : ntaps-1;
	int		L;
	{
		int	k;
		L = 0;
		for(k = 1; k<= ntaps/2; k++) {
			if ((k%M)==0)
				continue;
			L++;
		}
	}
	int		nz = L+1; // number of zeros in frequency
	int		ngrid = nextlg(ntaps*16)+1; // Grid points to examine
	printf("%d TAPS, M = %d, fp = %f, L = %d, nz = %d, ngrid=%d\n", ntaps, M, fp, L, nz, ngrid);
	double		dgrid;
	RLMATRIX	mx(nz,nz);
	RLVECTOR	x(nz), b(nz);
	int		*nodei;
	double		*fq, *Hf, *nodef, stopband, delta, old_delta;

	// fp ranges from 0 to 1/2M.
	stopband = 1./M -fp;  // Stopband is determined by the passband
	assert(fp < 1./M/2.); // Passband cutoff must be < 1/(2M)
	assert(1./M/2.< stopband);
	fq = new double[ngrid];
	Hf = new double[ngrid];
	nodei= new int[L+1];
	nodef= new double[L+1];

	// Create a frequecy grid
	dgrid = (0.5-stopband) / (ngrid-1);
	fq[0] = stopband;
	for(int k=1; k<ngrid; k++)
		fq[k] = fq[k-1] + dgrid;
	fq[ngrid-1] = 0.5;

	//
	// Select points for our linear equations across that grid
	nodef[0] = fq[0];
	nodei[0] = 0;
	for(int k=1; k<nz-1; k++) {
		nodei[k] = k*ngrid/nz;
		nodef[k] = fq[nodei[k]];
	}
	nodei[nz-1] = ngrid-1;
	nodef[nz-1] = fq[nodei[nz-1]];
	assert(nodei[nz-1] > nodei[nz-2]);

	// Iterate until "good enough"
	do {
		// Build our matrix
		for(int i=0; i<nz; i++) {
			// H(f0) = delta
			int k = 1;	// Coefficient location
			for(int j=0; j<L; j++) {
				mx[i][j] = 2.0*cos(nodef[i]*2.0*M_PI*k);
				k++;
				if ((k % M)==0)
					k++;
			} mx[i][L] = (i==0) ? 1. : -mx[i-1][L];
			b[i] = -1./M;	// Known -h[0] value
		}

		// Solve for our coefficients
		x = mx.solve(b);
		
		// Check if we are within our limits
		delta = fabs(x[nz-1]);
		printf("delta = %f -> %6.2f dB\n", delta,
			20*log(delta)/log(10.));
		if ((fabs(delta-old_delta)/delta)<1e-5) {
			break;
		}

		old_delta = delta;

		// Calculate the frequency response of the filter
		for(int i=0; i<ngrid; i++) {
			double	acc = 1./M;
			int k = 1;	// Coefficient location
			for(int j=0; j<L; j++) {
				acc += 2.0 * cos(fq[i]*2.0*M_PI*k)*x[j];
				k++;
				if ((k % M)==0)
					k++;
			}
			Hf[i] = acc;
		}

		// Calculate a new set of frequencies via Remez-Exchange

		// Treat the first node special: it *must* be at the limit
		nodei[0] = 0;
		nodef[0] = fq[nodei[0]];
		int	lastnode = 0;

		// Walk through the rest of the nodes
		//	First node, k=0 is fixed at the edge of the stopband
		//	Last node should be fixed at the other edge ... might
		//	  not be though.
		for(int k=1; k<nz; k++) {
			// Walk forward until the sign changes
			int	sgn;
			int	ni = lastnode + 1;

			if (lastnode == 0)
				sgn = 1;
			else
				sgn = (Hf[lastnode] > 0.0) ? 1 : -1;

			// Skip forward until the sign changes
			while((ni < ngrid)&&(Hf[ni]*sgn >= 0))
				ni++;

			// Now move forward until we find the maximum, before
			// the sign changes again
			sgn = -sgn;
			while((ni < ngrid)&&(Hf[ni]*sgn > Hf[ni-1]*sgn))
				ni++;

			nodei[k] = ni-1;
			assert(nodei[k] > nodei[k-1]);
			nodef[k] = fq[nodei[k]];
			lastnode = nodei[k];
		}
	} while(1);

	LPFIRLITE	result;
	{
		result = new FIRLITE(ntaps);

		int	eqn = 0, hloc=ntaps/2+1, lloc = hloc-2;

		for(int k=0; k<ntaps; k++)
			result[k] = 0.0;
		result[ntaps/2] = 1./M;
		for(int k=1; k<=ntaps/2; k++) {
			double	hk;
			if ((k%M) == 0) 
				hk = 0.0;
			else {
				printf("R[%4d] = R[%4d] = %11.5f\n",
					hloc, lloc, x[eqn]);
				hk = x[eqn++];
			}
			result[hloc] = hk;
			result[lloc] = hk;
			hloc++;
			lloc--;
		}
	}

	delete[]	fq;
	delete[]	Hf;
	delete[]	nodef;

	return result;
}