Frequency-sampling design

Choose a filter’s gain at a few frequencies, 1 in the pass band and 0 above it, and the inverse DFT gives the taps. Watch the curve between the dots: it passes through every one and swings in between.

Pick the samples, get the taps

In Window-method FIR design (19.1) I started in time: take the ideal impulse response, cut it to $N_h$ taps, and look at the gain that comes out. This page starts at the other end. I say what the gain should be at a few frequencies, and the inverse DFT finds the taps.

The idea comes from Sampling the spectrum (13.1). $N$ values of a spectrum at $\Omega_k=2\pi k/N$ describe the signal repeated every $N$ samples. When $N$ is at least the signal’s length, they describe the signal itself. So if I choose $N_h$ values of a gain and invert them, I get $N_h$ taps whose spectrum passes through those values.

Let’s take $N_h=33$, so the sample frequencies are $\Omega_k=2\pi k/33$ for $k=0$ to 32. I want a filter with linear phase, so I choose real numbers $H_\text{zp}(\Omega_k)$ for its zero-phase response from Linear-phase systems (17.3). That response is a sum of cosines, and a cosine has the same value at $-\Omega$ as at $\Omega$. By Frequency in discrete time (12.1), $\Omega_{33-k}=2\pi-\Omega_k$ is the same frequency as $-\Omega_k$, so the samples must be mirrored:

So only $k=0$ to 16 are free: 17 numbers on the range from 0 to $\pi$.

Next comes the delay. By 17.3, symmetric taps have $H(e^{j\Omega})=e^{-j\alpha\Omega}H_\text{zp}(\Omega)$ with $\alpha=(N_h-1)/2$, which is 16 here. So the DFT values I hand to the inverse DFT are

The inverse DFT of The DFT (13.2) turns them into taps:

Now pair each $k$ from 1 to 16 with its mirror $33-k$. The mirror has the same sample. Its arrow is $e^{j(2\pi-\Omega_k)(n-\alpha)}=e^{-j\Omega_k(n-\alpha)}$, because $n-\alpha$ is a whole number, and a whole number of turns changes nothing. Euler’s pair from Complex exponentials & phasors (3.4) adds the two arrows to $2\cos\big(\Omega_k(n-\alpha)\big)$, and the taps are

Every term is real and the same at $n=\alpha+m$ as at $n=\alpha-m$, so the taps are real and symmetric about $\alpha$. By 13.1, these 33 taps have exactly the 33 spectrum samples I chose. Designing a filter this way is called frequency-sampling design.

For a first low-pass I set the samples to 1 for $k=0$ to 4 and to 0 from $k=5$ on. Sample 4 sits at $\Omega_4=0.24\pi$ and sample 5 at $\Omega_5=0.30\pi$. One tap can be worked out by hand. At the centre, $n=\alpha$, every cosine is $\cos0=1$, so $h[16]=(1+2\cdot4)/33=9/33=0.273$.

Think of a rope tied to posts. At every post the rope is held at the height you chose. Between the posts it hangs however it likes. The picture at the top of the page shows the posts, the taps they give, and the rope.

Watch the curve between the dots. In the pass band it rises to 1.121 between two samples of exactly 1. Next to the jump from 1 to 0 the stop band only gets down to −15.9 dB, which is a gain of 0.16.

Why does the rope hang the way it does? Put the formula for $h[n]$ back into the zero-phase response of 17.3 and add up over $n$ first. Each sample then brings along one copy of a fixed curve:

The fraction is the Dirichlet kernel of The DTFT (12.2), the periodic sinc, moved to $\Omega_k$. Divided by $N_h$, it is 1 at its own sample and 0 at every other one, which is why the rope is held at each post. Between the posts the kernels’ ripples add. In a long run of equal samples they nearly cancel, but at a jump from 1 to 0 nothing on the far side cancels them. So the response rings, as 19.1’s plain cut did, and the cure is not to jump.

One sample in the gap sinks the ripple

Suppose the sample at the jump, $k=5$, is not 0 but a value between 0 and 1. The step from 1 to 0 then has a stair in the middle. I call such an in-between value a transition sample. Everything else stays as before, so the taps still come from the formula above.

The gentler step has a price. The stop band now starts one sample later, at $\Omega_6=0.36\pi$, so the transition is one sample wider. Think of a ramp beside a step: it is gentler, and it takes more room.

To compare fairly, I measure every version from $\Omega_6$. With the transition sample at 0 the filter is the one above, and from $\Omega_6$ on its highest point is −20.5 dB. The −15.9 dB peak sat at $0.33\pi$, before $\Omega_6$.

Below, the sample at $k=5$ moves from 0 up and back down. Watch the dotted level of the stop band’s highest point as it moves.

After the clip, drag the diamond up or down to set the transition sample yourself, or use the arrow keys. Try 0.25 and then 0.5. Both sit near −30 dB, and the best value lies between them, not at the halfway point. Then try values above 0.5: at 0.75 the stop band is back at −20.5 dB, as bad as the abrupt step, and at 1 it is −16.1 dB.

The best single value is 0.3934, which gives −42.5 dB. That is 22 dB lower than with no transition sample. The pass band pays a little too: at 0.39 it strays up to 0.051 from 1.

Two transition samples do better still. The best pair is 0.5979 at $k=5$ and 0.1106 at $k=6$, and the stop band, which now starts at $\Omega_7=0.42\pi$, peaks at −66.7 dB. The optimum is sharp: at 0.598 and 0.111 the peak is already −65.9 dB, and at 0.60 and 0.11 it is −61.7 dB.

How are such values found? Look again at the formula for the taps. Each tap is a sum of the samples times fixed numbers, so the gain at every $\Omega$ is too. Lowering the highest stop-band gain over a few free samples is then a problem called linear programming, and that is how Rabiner, Gold and McGonegal tabulated the best transition samples in 1970.

Worked example

Let’s meet the running spec of Filter specifications (18.1). At $f_s=8$ kHz, the gain must stay within $1\pm0.05$ from 0 to 1 kHz and at most 0.01, which is −40 dB, from 1.5 kHz to 4 kHz.

Take $N_h=31$. The samples are $8000/31=258.1$ Hz apart, so $k=0$ to 3 lie at 0, 258.1, 516.1 and 774.2 Hz. These are in the pass band, and each is set to 1. The next two, $k=4$ at 1032.3 Hz and $k=5$ at 1290.3 Hz, fall in the gap between 1 and 1.5 kHz, so they are two transition samples. From $k=6$, at 1548.4 Hz, every sample is 0.

Optimised by linear programming, the two transition samples come out as 0.909 and 0.296. The centre tap is then a hand calculation, as before. With $\alpha=15$ every cosine is 1 at $n=15$, so $h[15]=\big(1+2(3+0.909+0.296)\big)/31=0.304$.

On a fine grid, the pass band stays within 0.049 of 1, against the 0.05 allowed. The stop band peaks at −40.01 dB, a gain of 0.00999 against the 0.01 allowed. So 31 taps meet the spec, but only at the limit, with 0.01 dB to spare. Even letting the pass band use the full 0.05 only reaches −40.04 dB.

Two more taps give room. With $N_h=33$ the samples are 242.4 Hz apart, and the transition samples are 0.456 at 1212.1 Hz and 0.0226 at 1454.5 Hz. The pass band stays within 0.042 of 1, and the stop band peaks at −48.0 dB. Two fewer than 31, $N_h=29$, cannot meet the spec: the best stop band is −34.1 dB.

Compare that with Window-method FIR design (19.1), where Kaiser’s window needed 38 taps for the same spec. Frequency sampling gets there with fewer, but it depends on where the samples land against the band edges. With $N_h=35$ the stop band is back up to −40.3 dB, only just passing again, although 33 taps had 8 dB to spare.

Where you’ll meet this

Frequency sampling is the natural way to design a filter to a response that was measured or drawn by hand, such as the correction for a room or a loudspeaker. SciPy’s firwin2 takes a list of frequencies and gains and joins them by straight lines on a dense grid. It then inverts the DFT and multiplies the taps by a window, as in 19.1, to calm the ripple. Audio equalisers and biquads (20.5) returns to equalising a measured response.

The same samples also give a filter structure. A comb filter puts zeros at every $\Omega_k$, and one resonator per nonzero sample puts a pole back on top of its zero, as in Resonators, notches and combs (17.4). When only a few samples are nonzero, that is a cheap way to compute a narrow-band filter, and Filter banks (23.1) builds banks of them.

Here only one or two samples were free, and the rest were fixed at 1 or 0. Optimal FIR design (19.3) frees every tap at once, so that the largest error over both bands is as small as it can be, and compares every method’s shortest filter for the running spec.

Reference card