Remez exchange algorithm with flow chart

The Remez exchange algorithm (Parks–McClellan) designs an optimal equiripple linear-phase FIR filter. It iteratively adjusts a set of "extremal" frequencies until the weighted approximation error has equal-height ripples. It is the usual way to implement the Chebyshev (minimax) criterion.

Problem set-up

For a symmetric, odd-length filter (Type I), H(e^jω) = e^{−jωL} H_r(ω) with

H_r(ω) = Σ (k = 0 to L) a(k) cos kω, L = (M−1)/2

and the weighted error over the pass and stop bands (the set S):

E(ω) = W(ω) [H_d(ω) − H_r(ω)], ω ∈ S

The goal is to choose a(k) to minimise δ = max over S of |E(ω)|.

Derivation

By the alternation theorem, the best H_r(ω) is the one for which |E(ω)| reaches its maximum δ at at least L + 2 points ω₀ < ω₁ < … < ω_{L+1} in S, with alternating signs:

E(ω_i) = (−1)^i δ ⇒ H_r(ω_i) = H_d(ω_i) − (−1)^i δ / W(ω_i)

These are L + 2 linear equations in the L + 1 coefficients a(k) and δ. Using x_i = cos ω_i, they are solved without matrix inversion:

• δ = [Σ_{k=0}^{L+1} b_k H_d(ω_k)] / [Σ_{k=0}^{L+1} (−1)^k b_k / W(ω_k)], b_k = Π_{i≠k} 1/(x_k − x_i)
• For the first L + 1 points, C_k = H_d(ω_k) − (−1)^k δ/W(ω_k), the amplitude response at any x = cos ω is the Lagrange interpolation H_r(ω) = [Σ d_k C_k/(x − x_k)] / [Σ d_k/(x − x_k)], d_k = Π_{i≠k, i≤L} 1/(x_k − x_i)

Algorithm steps

1. Choose M and band edges; guess L + 2 extremal frequencies (equally spaced in the bands, including the band edges).
2. Compute δ from the formula above.
3. Compute H_r(ω) by interpolation and evaluate E(ω) on a dense grid (about 16 × M points).
4. Find the local maxima of |E(ω)| on the grid. If at some frequencies |E| > |δ|, take the L + 2 largest alternating extrema as the new set (this is the "exchange").
5. Repeat steps 2–4 until max |E| ≈ |δ| (the set no longer changes). Convergence is guaranteed.
6. Compute the impulse response by inverse DFT of the sampled H_r(ω): h(n) = (1/M)[H_r(0) + 2 Σ (k=1..L) H_r(2πk/M) cos(2πk(n−L)/M)].

Flow chart

        ┌───────────────────────────┐
        │ Inputs: M, edges, W(ω)    │
        └─────────────┬─────────────┘
                      ▼
        ┌───────────────────────────┐
        │ Initial guess of L+2      │
        │ extremal frequencies      │
        └─────────────┬─────────────┘
                      ▼
        ┌───────────────────────────┐ ◄─────────┐
        │ Compute δ                 │           │
        └─────────────┬─────────────┘           │
                      ▼                         │
        ┌───────────────────────────┐           │
        │ Interpolate H_r(ω); find  │           │
        │ error E(ω) on dense grid  │           │
        └─────────────┬─────────────┘           │
                      ▼                         │
        ┌───────────────────────────┐   no      │
        │ max|E(ω)| ≤ |δ| (+ tol)?  ├───────────┤
        └─────────────┬─────────────┘  new set  │
                      │ yes            of L+2   │
                      ▼                extrema ─┘
        ┌───────────────────────────┐
        │ Compute h(n) (IDFT); stop │
        └───────────────────────────┘

The result is an equiripple filter with the smallest possible maximum error for that M. For a fixed δ the required M can be estimated from M ≈ [−10 log₁₀(δ₁δ₂) − 13] / (14.6 Δf) + 1, where Δf = (ωs − ωp)/2π.