Optimal FIR design

What makes one filter the best? Below, two 27-tap designs face the same spec: one makes the average error smallest, the other the worst error. Watch which one stays out of the hatched zones.

Least squares or least worst

In Window-method FIR design (19.1), Kaiser’s window met the running spec with 38 taps. In Frequency-sampling design (19.2), two well-chosen transition samples did it with 31. Every tap costs one multiplication per output sample, so I want to know the shortest filter that can do it. To find it, I first have to say what “the best filter” means.

Here is the running spec of Filter specifications (18.1) once more. The sample rate is 8 kHz. From 0 to 1 kHz the gain must stay within $1\pm0.05$, and from 1.5 kHz to 4 kHz it must stay at or below 0.01, which is 40 dB down. In rad/sample the pass band ends at $\Omega_\text{pass}=0.25\pi$ and the stop band starts at $\Omega_\text{stop}=0.375\pi$.

A filter with an odd number of symmetric taps has the real zero-phase response of Linear-phase systems (17.3):

and its gain is $\lvert H_\text{zp}(\Omega)\rvert$. For 27 taps the centre is $\alpha=13$. Only $h[0]$ to $h[13]$ are free, because the other 13 taps are their mirror copies. So a 27-tap design is a choice of 14 numbers.

What do we want $H_\text{zp}$ to be? I call the target the desired response $H_\text{des}(\Omega)$: 1 in the pass band and 0 in the stop band. In the transition band I ask for nothing, and errors there do not count. The error of a design is then $H_\text{des}(\Omega)-H_\text{zp}(\Omega)$ in the two bands.

The spec is not equally strict in both bands. It allows 0.05 in the pass band but only 0.01 in the stop band, five times less. So I multiply the stop-band error by 5 before I judge it, and call this factor the weight $\gamma(\Omega)$: 1 in the pass band, 5 in the stop band. The weighted error is

With this weight, the whole spec becomes one condition: $\lvert\varepsilon(\Omega)\rvert\le0.05$ in both bands. In the pass band that says $\lvert1-H_\text{zp}\rvert\le0.05$. In the stop band it says $5\lvert H_\text{zp}\rvert\le0.05$, which is a gain of at most 0.01.

There are two natural ways to make $\varepsilon$ small. The first is to make its average size small. Least squares chooses the taps that make the mean square of $\varepsilon$ over the two bands as small as it can be:

This is the measure of Convergence and the Gibbs phenomenon (7.3), where a partial sum was the best approximation in mean square. The second way is to make the largest value small. Minimax chooses the taps that make the worst point as low as it can be:

You have already met one least-squares design. The plain cut of 19.1 kept 27 terms of the ideal $h_\text{des}$, and by 7.3 a partial sum is the best in mean square. So the plain cut is the least-squares answer with no transition band and $\gamma=1$. It was the best by that measure, and it still rang next to the edge.

Least squares is also easy to compute. 19.2 showed that the gain is linear in the design’s numbers, so $\varepsilon$ is linear in the 14 taps. Its mean square is then a quadratic in the 14 taps, and its lowest point comes from 14 linear equations, solved once. SciPy’s firls does this.

Which measure should we use? Think of two ways to get to work. One usually takes 20 minutes but sometimes 60, and the other always takes 30. If you must never be late, you take the second.

A spec is like that: it does not look at the average, only at the worst point. The picture at the top of the page draws the best 27-tap design by each measure against the running spec. Both use the weight 5.

After the clip, press either design’s button to draw it at full strength.

Compare the two “stop band, RMS” readouts first. Least squares wins on average, −50.4 dB against −44.9 dB, and it still fails: right at the stop edge its gain climbs to −31.7 dB. The equiripple design spreads its error evenly, so no ripple sticks out. For comparison, Kaiser’s window at 27 taps strays 0.069 in the pass band and reaches −22.4 dB.

Equal ripples

The equiripple design is not a lucky guess: it is what minimax gives. Why should the best worst point come with ripples that are all the same height? Suppose one ripple were lower than the others. Its spare room would be wasted, because the filter could give some of it up to push the highest ripple down.

A theorem makes this exact. Call the largest size of the weighted error $\delta$.

Take a type I filter (odd length, symmetric taps, 17.3) with $\alpha+1$ free numbers. The alternation theorem says it has the smallest possible $\delta$ if and only if its weighted error reaches $\pm\delta$ at least $\alpha+2$ times, alternating in sign. Only the best filter does this. A design whose ripples all have one height is called equiripple.

For 27 taps, $\alpha+2=15$, so the best design touches $+\delta,-\delta,+\delta,\dots$ at 15 frequencies. In terms of gain, the pass band strays at most $\delta_\text{pass}=\delta$ from 1, and the stop band reaches at most $\delta_\text{stop}=\delta/5$, because the weight there is 5. I will not prove the theorem here, but the next instrument lets you count the 15 touches.

Remez exchange: move the points to the peaks

How do we find the equiripple filter? The Parks–McClellan program does it with the Remez exchange, a search that repeats four steps:

1. Pick 15 trial frequencies $\Omega_0,\Omega_1,\dots,\Omega_{14}$ in the two bands.
2. Find the one filter whose weighted error is exactly $+\delta,-\delta,+\delta,\dots$ at them.
3. Look where its weighted error actually peaks.
4. Move the trial points to those peaks, and repeat until no peak is higher than $\delta$.

Step 2 is a set of linear equations. At each trial point, with $m$ counting the points,

$H_\text{zp}$ is linear in the taps, so these are 15 linear equations for 15 unknowns: the 14 numbers $h[0]$ to $h[13]$, and $\delta$.

Think of levelling a wobbly table by always fixing its highest corner. Each step fixes the worst place, and the next worst place becomes the target. Here, $\delta$ of a trial filter can never be larger than the best possible worst error, and the largest error of that filter can never be smaller. So the answer is caught between the two numbers, and the search stops when they meet.

The instrument below runs the search for the running spec at 27 taps, and shows the weighted error $\varepsilon(\Omega)$, where equiripple means literally equal ripples.

The two numbers at each iteration are these.

Watch the two readouts. The level $\delta$ grows at every step, and the largest error shrinks. In iteration 1 the dots are spread evenly, and between them the error is hundreds of times larger than $\delta$. Once the dots sit on the peaks, the two numbers close in fast, and in iteration 6 they agree to three decimals.

In the last frame, count the dots. Five sit in the pass band and ten in the stop band, and their signs alternate: 15 touches, as the theorem demands. The stop-band ripples look as tall as the pass-band ones only because of the weight. In gain they are five times smaller.

The weight moves the ripple between the bands

Why the weight 5? Call the weight’s value in the stop band $\gamma_\text{stop}$. The equiripple design makes the weighted error $\delta$ high in both bands, so the gain strays $\delta_\text{pass}=\delta$ in the pass band and $\delta_\text{stop}=\delta/\gamma_\text{stop}$ in the stop band. Their ratio is

So the weight decides how the error is shared between the bands, and the length decides how much error there is to share. It is like splitting a fixed budget between two teams: more for one means less for the other. The running spec allows $0.05/0.01=5$, so that is the weight to use.

The instrument below shows three 27-tap equiripple designs with different weights.

After the clip, set the weight yourself with the slider, from 1 to 50. Move the weight up and down and watch the two readouts. The pass-band ripple and the stop-band level always move in opposite directions.

Only weights 4 to 8 pass: the pass band strays 0.038 to 0.048, and the stop band sits at −40.5 to −44.5 dB. At weight 3 the stop band fails at −38.8 dB, and at weight 9 the pass band fails at 0.051. At every weight the two ripples keep the ratio $\gamma_\text{stop}$, within 1 %.

One spec, seven designs

Now we can answer the question from the start. For each method I increased the length until the design met the running spec. The delay of a linear-phase filter is $\alpha=(N_h-1)/2$ samples (17.3), and the last column gives it at 8 kHz.

The equiripple design is the shortest. SciPy’s remez also designs even lengths, which are linear phase too (17.3), and 26 taps are enough. Against Kaiser’s 38, that is a third fewer multiplications for every output sample.

Is there a formula for the length, as there was for the Kaiser window? Kaiser also gave an estimate for equiripple designs. With $\Delta f=f_\text{stop}-f_\text{pass}$ the transition width in hertz (18.1),

For the running spec it gives 22.9, against the true 26. Like Kaiser’s window formula in 19.1, it is an estimate to check.

Why not always use equiripple? It needs an iterative program, where a window design is one formula. And its stop band does not fall away.

The 27-tap design stays at −41.9 dB right up to 4 kHz. Kaiser’s 38 taps are at −41.7 dB next to the edge, but below −59.0 dB from 2500 Hz on and below −70.0 dB from 3500 Hz on. If most of the unwanted signal lies far from the edge, the window design removes more of it.

Worked example

Let’s design the running spec as an equiripple filter, step by step, and check the result.

The weight. The spec allows $\delta_\text{pass}=0.05$ and $\delta_\text{stop}=0.01$, so $\gamma_\text{stop}=0.05/0.01=5$.

A first length. In Kaiser’s estimate, $\sqrt{0.05\cdot0.01}=0.02236$, and $-20\log_{10}0.02236=33.01$ dB. The transition is $\Delta f/f_s=500/8000$, so the denominator is $14.6\cdot0.0625=0.9125$. Then

Check and adjust. Start at 23 taps and add taps until the design passes. At 23 taps the pass band strays 0.078 and the stop band reaches −36.1 dB, so both fail. At 25 taps it is 0.056 and −39.1 dB, still failing. At 27 taps it is 0.040 and −41.9 dB, and the design passes.

With even lengths allowed, 26 taps pass with 0.044 and −41.1 dB.

The result at 27 taps. The Remez exchange takes 7 iterations from the even start and ends at $\delta=0.04004$. So $\delta_\text{pass}=0.0400$ and $\delta_\text{stop}=0.0400/5=0.0080$, which is −41.9 dB. The weighted error touches $\pm\delta$ 15 times, alternating, so no 27-tap filter does better with this weight. The delay is $\alpha=13$ samples, which is 1.625 ms at 8 kHz.

The alternatives at 27 taps. Least squares with the same weight fails at −31.7 dB, and it needs 39 taps to pass (43 with weight 1). Equiripple with weight 1, the same ripple in both bands, needs 35 odd taps, or 34 with an even length. The weight alone is worth 8 taps here.

Where you’ll meet this

scipy.signal.remez and MATLAB’s firpm are the standard tools for filters in fixed-point processors and hardware, where every tap costs. At 8 kHz, the 26-tap design needs 208 000 multiplications a second, and Kaiser’s 38 taps need 304 000. On a chip, fewer taps also means less memory and less power.

Least squares has its place too. When the unwanted signal is broadband noise, the noise power that leaks through depends on the mean square of the stop-band gain. That is the quantity least squares keeps small, and SciPy’s firls designs it in one step.

The same equal-ripple idea returns for recursive filters: elliptic filters are the equiripple designs among IIR filters, in Analog prototype filters (20.1). Choosing FIR or IIR (20.6) compares them for one spec. Long equiripple low-pass filters are the usual choice in Downsampling and decimation (22.1), and Finite word-length effects (21.3) shows what happens to their ripples when the taps are rounded.

Reference card