This section describes the algorithm used by gmeteor to design a
filter. This algorithm was originally described in the paper
METEOR: a Constraint-based FIR Filter Design Program, by K.
Steiglitz, T. W. Parks, and J. F. Kaiser, published in IEEE
Trans. Signal Processing, vol. 40, no. 8, pp. 1901-1909, August
1992.
The inputs to the algorithm are the filter length l, the symmetry (either sine or cosine symmetry), the sampling frequency, the grid density, and a list of specifications. Each specification consists of a band B, a weight w(f), and a constraint that has one of the following six forms: H(f) <= u(f), H(f) >= u(f), H'(f) <= 0, H'(f) >= 0, H”(f) <= 0, or H”(f) >= 0. We think of each specification as a compact representation for an infinite number of constraints, one for each frequency f in the band B.
The algorithm begins by rewriting the specifications as follows. Constraints of the form H(f) <= u(f) are rewritten as H(f) + w(f) y <= u(f). Constraints of the form H(f) >= u(f) are rewritten as H(f) - w(f) y >= u(f). In the other four cases (specifications on the derivatives of H(f)) the constraint is left unchanged and the weight is ignored.
Next, from the filter length, symmetry, and sampling frequency,
gmeteor derives an analytic expression for H(f) in terms of
the filter coefficients h[i]. This analytic
expression is the discrete-time Fourier transform of the
h[i]'s, but without diving into too many details, we can
simply state that H(f) has the general form
where t is some function. The important property is that, for a
fixed f, H(f) is a linear combination of the filter
coefficients.
At this point, gmeteor solves the following optimization problem:
find the filter coefficients h[i]'s so as to maximize
y, subject to the constraints rewritten as discussed above. This
optimization problem has an infinite number of constraints. To solve
it, gmeteor samples the constraints on a finite grid of frequencies.
The grid is determined by the grid-density parameter, which
specifies the number of grid points in the range [0..F/2]. After
this sampling, the filter-design problem has been reduced a
linear-programming problem, which gmeteor solves by means of the dual
simplex algorithm. (I believe it is possible to solve the unsampled
continuous problem using symbolic algebra and a modification of the
revised dual simplex algorithm, but the finite algorithm works well
already.)