Filter Structures

Realizing a transfer function: direct, transposed, cascade, parallel, lattice

Chapter 6 ended with a set of coefficients: a transfer function $H(z)$ that meets a specification. That is the what. This chapter is the how: given $H(z)$, how do you actually wire up the delays, multipliers, and adders that compute it sample by sample?

The surprise is that there is no single answer. The same $H(z)$ can be realized by many different block diagrams, all mathematically identical in exact arithmetic, yet very different once you build them. They differ in how much memory they need, how many operations per sample, and (most importantly) how gracefully they degrade when the coefficients and signals are stored in finite-precision numbers. Choosing a structure is the last design decision before the filter becomes code or silicon. It is the real start of the “math to metal” pipeline.

We work with the general difference equation from Chapter 2,

\[y[n] = \sum_{k=0}^{M} b_k\, x[n-k] \;-\; \sum_{k=1}^{N} a_k\, y[n-k],\]

whose transfer function is

\[H(z) = \frac{\sum_{k=0}^{M} b_k\, z^{-k}}{1 + \sum_{k=1}^{N} a_k\, z^{-k}} = \frac{B(z)}{A(z)}.\]

Every structure in this chapter realizes this same $H(z)$. They are different circuits for one function.

Why structure matters

In exact arithmetic, rearranging a signal-flow graph that computes $H(z)$ cannot change its output. Real hardware is not exact. Coefficients are rounded to a finite word length, and intermediate sums are quantized. Two effects make the choice of structure consequential:

Coefficient sensitivity. The poles of a high-order $A(z)$ can be extremely sensitive to small changes in the polynomial coefficients. Rounding a single $a_k$ by one least-significant bit can move a pole enough to ripple the passband or even push it outside the unit circle, making a “stable” design unstable. Breaking the filter into small sections makes each section’s poles far less sensitive.
Dynamic range and limit cycles. Different structures produce different intermediate signal magnitudes for the same input. A structure whose internal nodes swing widely is prone to overflow in fixed point; feedback combined with quantization can also trap an IIR filter in a small self-sustaining oscillation called a limit cycle.

So the questions a structure answers are practical: How many delay registers? How many multiply-accumulates per sample? How robust is it to coefficient rounding? Does it overflow? The rest of the chapter introduces the structures in roughly increasing order of numerical robustness.

Direct forms

The most literal realization writes the difference equation out term by term, one multiplier per coefficient and one delay per $z^{-1}$. This is the direct form. Chapter 2 already drew the generalized FIR transversal structure and the IIR Direct Form I (DF-I), which keeps separate delay lines for the input history and the output history.

DF-I uses $M + N$ delay elements. You can do better: because the all-zero section $B(z)$ and the all-pole section $1/A(z)$ are in series, and series LTI systems commute, you can swap their order and then share one delay line between them. That is Direct Form II (DF-II), the canonical form, which uses only $\max(M, N)$ delays, the theoretical minimum. The biquad page works through DF-I, DF-II, and the transposed form for the second-order case in detail, with diagrams. Here we focus on what is general about structures, starting with the transformation that turns one structure into another.

The transposition theorem

Every structure has a twin that realizes the same transfer function. You obtain it by transposing the signal-flow graph, following three steps (the transposition theorem):

Reverse the direction of every branch.
Replace every branch (pickoff) node with an adder, and every adder with a branch node.
Interchange the input and the output.

Reversing a flow graph leaves its transfer function unchanged, so the transposed structure computes exactly the same $H(z)$ with the same coefficients, but it distributes the arithmetic differently and therefore has different numerical behaviour.

Applied to DF-II, this gives the transposed Direct Form II (TDF-II). It is the workhorse of practical IIR filtering: its state update accumulates smaller intermediate values than DF-II, it pipelines naturally, and it is the form that SciPy’s sosfilt and ARM’s CMSIS-DSP arm_biquad_cascade_df2T use internally. The biquad page shows its diagram and state equations.

The point of transposition

Transposition is free: it costs nothing to apply and never changes $H(z)$. It is purely a tool for trading numerical properties. When someone says “use the transposed form,” they mean “same filter, friendlier arithmetic.”

Cascade form

A single high-order $(b, a)$ pair is the most fragile way to implement a filter, precisely because of the coefficient sensitivity described above. The standard fix is to factor the transfer function into a product of low-order sections:

\[H(z) = g \prod_{k=1}^{L} H_k(z), \qquad H_k(z) = \frac{b_{0k} + b_{1k}\,z^{-1} + b_{2k}\,z^{-2}}{1 + a_{1k}\,z^{-1} + a_{2k}\,z^{-2}}.\]

Each $H_k(z)$ is a second-order section (a biquad), and the sections run in series, the output of one feeding the next.

Figure 1: Cascade form. The transfer function is factored into a series of low-order sections; in practice each H_k(z) is a biquad.

An order-$N$ filter needs $\lceil N/2 \rceil$ sections. The win is robustness: a biquad’s two poles are easy to place accurately, and rounding its handful of coefficients barely moves them, whereas rounding the coefficients of one big polynomial can move its poles dramatically. Cascading also makes stability trivial to check section by section.

In Python, scipy.signal.tf2sos performs the factorization and pole-zero pairing for you, and sosfilt runs the cascade (each section in transposed Direct Form II):

Show the code

# A 4th-order lowpass, expressed two ways
b, a = butter(4, 0.25)                 # single (b, a) pair
sos = tf2sos(b, a)                     # cascade of second-order sections
print("sections:", sos.shape[0])
print(np.round(sos, 4))

# The two realizations agree to machine precision (at this modest order)
rng = np.random.default_rng(0)
x = rng.standard_normal(2000)
err = np.max(np.abs(sosfilt(sos, x) - lfilter(b, a, x)))
print(f"max |sosfilt - lfilter| = {err:.2e}")

sections: 2
[[ 0.0102  0.0204  0.0102  1.     -0.8554  0.2097]
 [ 1.      2.      1.      1.     -1.113   0.5741]]
max |sosfilt - lfilter| = 5.77e-15

Each row of the sos matrix is one section, laid out as $[\,b_0, b_1, b_2, 1, a_1, a_2\,]$. At fourth order both realizations agree, but as the order grows the single-polynomial path loses accuracy fast while the cascade stays clean. The filter design chapter and the biquad page show that divergence directly. This is why every serious IIR implementation ships as second-order sections.

Section ordering matters

The factorization is not unique: you can pair any pole with any zero and order the sections however you like. These choices do not change $H(z)$ but they do change the internal signal levels, which matters for overflow in fixed point. tf2sos uses a sensible default (pair each pole with its nearest zero, order to keep intermediate gains moderate). For most work the default is fine.

Parallel form

Instead of factoring $H(z)$ into a product, you can expand it into a sum. A partial-fraction expansion writes the transfer function as

\[H(z) = \sum_{k} H_k(z) \;+\; (\text{optional direct term}),\]

where each $H_k(z)$ is again a low-order (typically second-order) section, but now the sections run in parallel: the input feeds all of them at once and their outputs are summed.

Figure 2: Parallel form. The input drives every section simultaneously and the section outputs are summed.

scipy.signal.residuez computes the expansion (residues $r_k$, poles $p_k$, and any direct terms). Grouping each complex-conjugate pole pair into one real second-order section gives a realizable parallel structure:

Show the code

def parallel_sections(b, a):
    """Partial-fraction expansion grouped into real second-order sections."""
    r, p, k = residuez(b, a)
    sections, used = [], np.zeros(len(p), bool)
    for i in range(len(p)):
        if used[i]:
            continue
        # find the conjugate partner of a complex pole
        j = next((jj for jj in range(len(p)) if jj != i and not used[jj]
                  and np.isclose(p[jj], np.conj(p[i])) and abs(p[i].imag) > 1e-9), -1)
        if j >= 0:                                   # conjugate pair -> real biquad
            used[i] = used[j] = True
            bsec = np.real([r[i] + r[j], -(r[i]*p[j] + r[j]*p[i]), 0.0])
            asec = np.real([1.0, -(p[i] + p[j]), p[i]*p[j]])
        else:                                        # lone real pole -> first order
            used[i] = True
            bsec, asec = np.real([r[i], 0.0, 0.0]), np.real([1.0, -p[i], 0.0])
        sections.append((bsec, asec))
    return sections, np.real(k)

b, a = butter(4, 0.25)
sections, direct = parallel_sections(b, a)
y_par = sum(lfilter(bs, as_, x) for bs, as_ in sections)
if len(direct):
    y_par = y_par + lfilter(direct, [1.0], x)
print(f"parallel sections: {len(sections)}")
print(f"max |parallel - lfilter| = {np.max(np.abs(y_par - lfilter(b, a, x))):.2e}")

parallel sections: 2
max |parallel - lfilter| = 4.44e-15

The cascade and parallel forms realize the same filter and share the same robustness advantage over a single high-order section. They differ in detail: the cascade preserves the exact zeros of $H(z)$ in each section (good when the zeros matter, as in notch filters), while the parallel form keeps the sections fully independent (good for parallel hardware, and the natural form when you think of a filter as a sum of resonances).

Lattice structures

The lattice is a different shape entirely. Instead of a delay line tapped by coefficients, it is a ladder of identical two-input, two-output stages, each parametrized by a single reflection coefficient $k_m$. Lattices are prized in speech coding (linear predictive coding), adaptive filtering, and anywhere coefficient robustness and a cheap stability test matter.

Each stage carries a forward signal $f_m[n]$ and a backward signal $g_m[n]$. For an all-zero (FIR) filter the stages run from $m = 1$ to $M$:

\[f_0[n] = g_0[n] = x[n],\] \[f_m[n] = f_{m-1}[n] + k_m\, g_{m-1}[n-1],\] \[g_m[n] = k_m\, f_{m-1}[n] + g_{m-1}[n-1],\]

and the output is $y[n] = f_M[n]$. A single stage is the cross-coupled “butterfly” below: the forward signal passes straight through and also feeds the backward adder through $k_m$, while the delayed backward signal passes straight through and feeds the forward adder through $k_m$.

Figure 3: One stage of an FIR (all-zero) lattice. The forward paths $f_m$ and backward paths $g_m$ are cross-coupled through the reflection coefficient $k_m$; the delay sits on the backward path.

The reflection coefficients are not the filter coefficients $b_k$; you obtain them from the polynomial by a step-down (Levinson) recursion. Starting from the monic polynomial $A(z) = 1 + a_1 z^{-1} + \cdots + a_M z^{-M}$, peel off one stage at a time:

\[k_m = a_m^{(m)}, \qquad a_i^{(m-1)} = \frac{a_i^{(m)} - k_m\, a_{m-i}^{(m)}}{1 - k_m^2}.\]

The same lattice stages, rewired so the forward recursion runs backward, realize the all-pole filter $1/A(z)$:

\[f_M[n] = x[n], \qquad f_{m-1}[n] = f_m[n] - k_m\, g_{m-1}[n-1], \qquad g_m[n] = k_m\, f_{m-1}[n] + g_{m-1}[n-1],\]

with output $y[n] = f_0[n]$. The code below builds both lattices from their reflection coefficients and confirms they reproduce the ordinary direct-form output to machine precision (the near-zero errors are the verification; this chapter introduces lattices, which scipy.signal does not provide ready-made):

Show the code

def poly_to_reflection(a):
    """Step-down (Levinson) recursion: monic polynomial -> reflection coefficients."""
    a = np.asarray(a, float)
    a = a / a[0]
    M = len(a) - 1
    k = np.zeros(M)
    for m in range(M, 0, -1):
        k[m - 1] = a[m]
        a = (a[:m + 1] - a[m] * a[:m + 1][::-1])[:m] / (1 - a[m] ** 2)
    return k

def fir_lattice(k, x):
    """All-zero FIR lattice; returns f_M[n]."""
    M = len(k)
    g_delayed = np.zeros(M)                       # g_{m-1}[n-1]
    y = np.zeros_like(x, float)
    for n in range(len(x)):
        f, g = x[n], np.empty(M + 1)
        g[0] = x[n]
        for m in range(1, M + 1):
            f_next = f + k[m - 1] * g_delayed[m - 1]
            g[m] = k[m - 1] * f + g_delayed[m - 1]
            f = f_next
        y[n] = f
        g_delayed = g[:M].copy()
    return y

def allpole_lattice(k, x):
    """All-pole IIR lattice realizing 1/A(z); returns f_0[n]."""
    M = len(k)
    g_delayed = np.zeros(M)
    y = np.zeros_like(x, float)
    for n in range(len(x)):
        f, g = x[n], np.empty(M + 1)
        for m in range(M, 0, -1):
            f = f - k[m - 1] * g_delayed[m - 1]
            g[m] = k[m - 1] * f + g_delayed[m - 1]
        g[0] = f
        y[n] = f
        g_delayed = g[:M].copy()
    return y

# FIR lattice vs direct-form FIR
b_fir = np.array([1.0, -0.6, 0.3, 0.2])
k_fir = poly_to_reflection(b_fir)
err_fir = np.max(np.abs(b_fir[0] * fir_lattice(k_fir, x) - lfilter(b_fir, [1.0], x)))
print(f"FIR  reflection coeffs = {np.round(k_fir, 4)}")
print(f"FIR  max |lattice - lfilter| = {err_fir:.2e}")

# All-pole lattice vs direct-form IIR
a_iir = np.array([1.0, -0.5, 0.25])
k_iir = poly_to_reflection(a_iir)
err_iir = np.max(np.abs(allpole_lattice(k_iir, x) - lfilter([1.0], a_iir, x)))
print(f"IIR  reflection coeffs = {np.round(k_iir, 4)}")
print(f"IIR  max |lattice - lfilter| = {err_iir:.2e}")

FIR  reflection coeffs = [-0.4783  0.4375  0.2   ]
FIR  max |lattice - lfilter| = 8.88e-16
IIR  reflection coeffs = [-0.4   0.25]
IIR  max |lattice - lfilter| = 4.44e-16

Two properties make the lattice attractive on hardware. First, stability is trivial to test: an all-pole lattice is stable if and only if every reflection coefficient satisfies $|k_m| < 1$, a single comparison per stage rather than a root-finding problem. Second, the lattice is modular and well-conditioned: increasing the order just appends a stage without recomputing the earlier ones, and the structure tolerates coefficient rounding gracefully, which is why it dominates speech LPC.

From math to metal

ARM’s CMSIS-DSP library ships lattice filters ready to drop onto a Cortex-M: arm_fir_lattice_f32 and arm_iir_lattice_f32 (with q15/q31 fixed-point variants). The FIR function takes exactly the reflection coefficients computed above; the IIR function realizes a general pole-zero response, so it needs the reflection coefficients plus a set of ladder (tap) coefficients for the numerator, with the all-pole case here as its special instance. The transposed Direct Form II biquad cascade from earlier is arm_biquad_cascade_df2T_f32. Picking a structure here is picking the library function you will call on the metal.

Transposition on an FPGA: where the pipeline registers go

The transposition theorem said the two forms compute the same $H(z)$ with the same coefficients, and that the choice is about arithmetic rather than mathematics. On an FPGA that choice has a name: the critical path, the longest run of combinational logic between two registers, which is what sets the maximum clock frequency.

Take a direct-form FIR. Every tap’s product has to be summed before the output register, so the adders form a chain (or a tree) between the multipliers and the output: the critical path grows with the number of taps, and a 64-tap filter clocks slower than an 8-tap one. Transpose it and the delays land between the adders instead of ahead of them. Each stage becomes one multiplier and one adder feeding its own register, so the longest path is a single multiply plus a single add whatever the tap count is. A 64-tap transposed FIR clocks at the same rate as an 8-tap one; it just uses more of the chip.

VHDL: transposed direct-form FIR

library ieee;
use ieee.std_logic_1164.all;
use ieee.numeric_std.all;

-- Transposed direct-form FIR. One register after every adder, so the longest
-- combinational path is ONE multiply plus ONE add regardless of N.
entity fir_transposed is
  generic (
    N     : positive := 8;    -- number of taps
    B_IN  : positive := 16;   -- input word width
    B_C   : positive := 16;   -- coefficient width
    B_ACC : positive := 34    -- = B_IN + B_C - 1 + ceil(log2 N); never overflows
  );
  port (
    clk   : in  std_logic;
    rst   : in  std_logic;                       -- synchronous, active high
    x_in  : in  signed(B_IN-1 downto 0);         -- one sample per clk
    y_out : out signed(B_ACC-1 downto 0)         -- full-precision output
  );
end entity;

architecture rtl of fir_transposed is
  type coeff_array is array (0 to N-1) of signed(B_C-1 downto 0);
  -- Symmetric low-pass, for illustration. Any tap set works unchanged, but
  -- this aggregate has to be restated if you change N.
  constant COEFF : coeff_array := (
    to_signed(  -412, B_C), to_signed(  -906, B_C),
    to_signed(  1247, B_C), to_signed(  9058, B_C),
    to_signed(  9058, B_C), to_signed(  1247, B_C),
    to_signed(  -906, B_C), to_signed(  -412, B_C));

  type acc_array is array (0 to N-1) of signed(B_ACC-1 downto 0);
  signal acc : acc_array := (others => (others => '0'));
begin
  process(clk)
  begin
    if rising_edge(clk) then
      if rst = '1' then
        acc <= (others => (others => '0'));
      else
        -- Last stage has nothing upstream of it: it starts a fresh partial sum.
        acc(N-1) <= resize(x_in * COEFF(N-1), B_ACC);
        -- Every other stage adds its own product to the partial sum arriving
        -- from the stage to its right. acc(k+1) reads its PRE-EDGE value, which
        -- is exactly the pipeline register that makes the path short.
        for k in N-2 downto 0 loop
          acc(k) <= acc(k+1) + resize(x_in * COEFF(k), B_ACC);
        end loop;
      end if;
    end if;
  end process;

  y_out <= acc(0);
end architecture;

Note what is not here: no delay line. The direct form needs an $N$-deep shift register holding past inputs; the transposed form holds partial sums instead, and the same sample is broadcast to every multiplier at once. The storage is wider (B_ACC rather than B_IN per stage) but the structure is otherwise identical in cost, which is the hardware statement of “transposition is free”.

Size B_ACC as $B_\text{in} + B_c - 1 + \lceil \log_2 N \rceil$ and no stage can overflow: each product needs $B_\text{in} + B_c - 1$ bits and summing $N$ of them adds $\lceil \log_2 N \rceil$. The model below checks that bound and the transposition claim at the same time, bit for bit.

import numpy as np

N, B_IN, B_C = 8, 16, 16
B_ACC = B_IN + B_C - 1 + int(np.ceil(np.log2(N)))          # = 34 bits
COEFF = np.array([-412, -906, 1247, 9058, 9058, 1247, -906, -412])

def fir_transposed(x):
    """Bit-accurate model of the VHDL entity: N accumulator registers, one clock per sample."""
    acc = [0] * N
    out = []
    for xn in x:
        nxt = [0] * N
        nxt[N - 1] = int(xn) * int(COEFF[N - 1])           # fresh partial sum
        for k in range(N - 2, -1, -1):
            nxt[k] = acc[k + 1] + int(xn) * int(COEFF[k])  # pre-edge read
        acc = nxt
        out.append(acc[0])
    return np.array(out, dtype=object), acc

rng = np.random.default_rng(0)
x = rng.integers(-2**(B_IN-1), 2**(B_IN-1), size=500)

y_hw, _ = fir_transposed(x)
y_ref = np.convolve(x, COEFF)[:len(x)]                     # the direct form, exactly

print(f"B_ACC needed = {B_ACC} bits, largest |acc| seen = {int(max(abs(y_hw)))}")
print(f"transposed == direct : {np.array_equal(y_hw.astype(np.int64), y_ref)}")

# Transposition changes the arithmetic, not the answer: bit-exact, not merely close.
assert np.array_equal(y_hw.astype(np.int64), y_ref), "transposed form must equal the direct form"
# And the width bound must hold with room to spare for the worst case, not just this input.
worst = N * 2**(B_IN-1) * 2**(B_C-1)
assert worst <= 2**(B_ACC-1), "B_ACC must cover the worst-case sum"
assert int(max(abs(y_hw))) < 2**(B_ACC-1), "this run must stay inside the accumulator"

B_ACC needed = 34 bits, largest |acc| seen = 608218647
transposed == direct : True

Testbench: check it against the direct form

The Python above proves the equivalence outside the simulator. Here is the same check inside it, which is the version you can run. The reference model is a direct-form FIR written the obvious way – keep the last $N$ inputs, multiply and sum them all in one clock – so the testbench compares the two structures from the transposition theorem against each other, cycle by cycle, and fails on the first sample where they part. Stimulus is the Galois LFSR from the noise generation page, which keeps the whole thing self-contained: no vector files, no scripts, nothing outside these two entities.

library ieee;
use ieee.std_logic_1164.all;
use ieee.numeric_std.all;
use std.env.all;                                 -- VHDL-2008, for stop

entity tb_fir_transposed is
end entity;

architecture sim of tb_fir_transposed is
  constant N         : positive := 8;
  constant B_IN      : positive := 16;
  constant B_C       : positive := 16;
  constant B_ACC     : positive := 34;
  constant N_SAMPLES : positive := 2000;

  type coeff_array is array (0 to N-1) of signed(B_C-1 downto 0);
  constant COEFF : coeff_array := (
    to_signed(  -412, B_C), to_signed(  -906, B_C),
    to_signed(  1247, B_C), to_signed(  9058, B_C),
    to_signed(  9058, B_C), to_signed(  1247, B_C),
    to_signed(  -906, B_C), to_signed(  -412, B_C));

  signal clk   : std_logic := '0';
  signal rst   : std_logic := '1';
  signal x_in  : signed(B_IN-1 downto 0)  := (others => '0');
  signal y_dut : signed(B_ACC-1 downto 0);
  signal y_ref : signed(B_ACC-1 downto 0) := (others => '0');
begin

  clk <= not clk after 5 ns;                     -- 100 MHz; the rate is arbitrary here

  dut : entity work.fir_transposed
    generic map (N => N, B_IN => B_IN, B_C => B_C, B_ACC => B_ACC)
    port map (clk => clk, rst => rst, x_in => x_in, y_out => y_dut);

  -- Direct form: the long adder chain, all N products summed in one clock.
  -- Correct, and exactly what the transposed form must reproduce.
  reference : process(clk)
    type sample_array is array (0 to N-2) of signed(B_IN-1 downto 0);
    variable dl  : sample_array := (others => (others => '0'));
    variable acc : signed(B_ACC-1 downto 0);
  begin
    if rising_edge(clk) then
      if rst = '1' then
        dl    := (others => (others => '0'));
        y_ref <= (others => '0');
      else
        -- Sum from the CURRENT sample and the delay line as it stands before
        -- this edge's shift, so that y_ref carries the same one-clock latency
        -- as the DUT.
        acc := resize(x_in * COEFF(0), B_ACC);
        for k in 1 to N-1 loop
          acc := acc + resize(dl(k-1) * COEFF(k), B_ACC);
        end loop;
        y_ref <= acc;

        for k in N-2 downto 1 loop               -- then shift, newest at 0
          dl(k) := dl(k-1);
        end loop;
        dl(0) := x_in;
      end if;
    end if;
  end process;

  -- Both structures reset to all zeros, so they must agree from the very first
  -- clock after reset: there is no priming window to skip past.
  checker : process(clk)
    variable cycle : natural := 0;
  begin
    if rising_edge(clk) and rst = '0' then
      cycle := cycle + 1;
      assert y_dut = y_ref
        report "MISMATCH at cycle " & integer'image(cycle) &
               ": transposed /= direct"
        severity failure;
    end if;
  end process;

  stimulus : process
    variable lfsr : std_logic_vector(15 downto 0) := x"ACE1";
  begin
    wait until rising_edge(clk);
    wait until rising_edge(clk);
    rst <= '0';

    for i in 1 to N_SAMPLES loop
      wait until rising_edge(clk);
      if lfsr(0) = '1' then                      -- Galois step, x"B400"
        lfsr := ('0' & lfsr(15 downto 1)) xor x"B400";
      else
        lfsr := '0' & lfsr(15 downto 1);
      end if;
      x_in <= signed(lfsr);                      -- full-scale, exercises B_ACC
    end loop;

    wait until rising_edge(clk);
    report "PASS: " & integer'image(N_SAMPLES) &
           " samples, transposed = direct on every cycle" severity note;
    stop;
  end process;

end architecture;

Two lines to run it, assuming both files are in the working directory:

ghdl -a --std=08 fir_transposed.vhd tb_fir_transposed.vhd
ghdl -r --std=08 tb_fir_transposed

A pass prints the note; a failure names the cycle it first diverged on and stops there. The severity failure matters: without it a mismatch is a printed line the run scrolls past, and a testbench you have to read the output of is a testbench that will eventually be believed while broken.

The same harness works for the other three entities on this site – the CIC decimator, the moving average and the LFSR – with the reference process swapped for whatever independent computation of the same quantity you trust. For the LFSR that reference is not a second implementation at all but the period itself: assert that the state returns to its seed after exactly $2^{16}-1$ clocks and never before.

Not simulated

The entity and its testbench are faithful transcriptions of the structure checked in Python, but neither has been through a simulator here – including the testbench itself, which has not been run through ghdl – so both are verified by inspection and against the model above. Run the two commands above before trusting any of it on hardware.

Choosing a structure

All of these realize the same $H(z)$. The differences are practical:

Structure	Delays	Robustness to coefficient rounding	Typical use
Direct Form I	$M + N$	Poor at high order; good per low-order section	Fixed-point biquads (separate delay lines resist overflow)
Direct Form II	$\max(M,N)$	Poor at high order; wide internal range	Rarely used directly at high order
Transposed DF-II	$\max(M,N)$	Good per section; small internal values	Default for floating-point (SciPy `sosfilt`, CMSIS `df2T`)
Cascade (SOS)	$2\lceil N/2\rceil$	Very good (low-order sections, exact zeros)	The standard production IIR form
Parallel	$2\lceil N/2\rceil$	Very good (independent sections)	Parallel hardware; sum-of-resonances designs
Lattice	$M$ or $N$	Very good; trivial $\lvert k_m\rvert < 1$ stability test	Speech LPC, adaptive filters

The headline advice is simple. Never implement a high-order IIR filter as a single $(b, a)$ pair. Factor it into second-order sections and run them in transposed Direct Form II, which is exactly what scipy.signal.sosfilt and CMSIS-DSP give you by default. Reach for lattice structures when you need a cheap stability guarantee or you are doing linear prediction.

Exercises

Exercise 1: Read a structure

A first-order structure is described by the difference equation \[y[n] = b_0\, x[n] + b_1\, x[n-1] - a_1\, y[n-1].\] Write its transfer function $H(z) = Y(z)/X(z)$. Is this Direct Form I or Direct Form II if the input and output histories use separate delay registers?

Answer

Taking the $z$-transform of both sides and collecting terms, \[H(z) = \frac{b_0 + b_1 z^{-1}}{1 + a_1 z^{-1}}.\] Separate delay registers for the input ($x[n-1]$) and the output ($y[n-1]$) make this Direct Form I. Sharing a single register between the feedforward and feedback paths would make it Direct Form II.

Exercise 2: Biquad difference equation and resonance

A biquad section has \[H(z) = \frac{1 + 0.1\,z^{-1} - 0.7199\,z^{-2}}{1 - 1.1786\,z^{-1} + 0.7246\,z^{-2}}.\] (a) Write its Direct Form II difference equations. (b) With $f_s = 100$ kHz, find the resonance frequency from the pole angle.

Answer

(a) With the shared state $w[n]$ and the denominator coefficients $a_1 = -1.1786$, $a_2 = 0.7246$: \[w[n] = x[n] + 1.1786\,w[n-1] - 0.7246\,w[n-2],\] \[y[n] = w[n] + 0.1\,w[n-1] - 0.7199\,w[n-2].\]

(b) Write the denominator as $A(z) = 1 - 2R\cos\theta\,z^{-1} + R^2 z^{-2}$. Then $R = \sqrt{0.7246} = 0.851$ and $\cos\theta = 1.1786 / (2R) = 0.692$, so $\theta = 0.806$ rad. The resonance frequency is \[f_c = \frac{\theta}{2\pi}\, f_s = \frac{0.806}{2\pi}\,(100\ \text{kHz}) \approx 12.8\ \text{kHz}.\] The poles have magnitude $R = 0.851 < 1$, so the section is stable.

Exercise 3: Cascade decomposition in Python

Take the 6th-order lowpass b, a = butter(6, 0.2). Decompose it into second-order sections with tf2sos, and compare sosfilt against lfilter on a noise input. Which one do you trust at this order, and why?

Answer

Show the code

b6, a6 = butter(6, 0.2)
sos6 = tf2sos(b6, a6)
xt = rng.standard_normal(4000)
diff = np.max(np.abs(sosfilt(sos6, xt) - lfilter(b6, a6, xt)))
print(f"{sos6.shape[0]} sections, max |sosfilt - lfilter| = {diff:.2e}")

3 sections, max |sosfilt - lfilter| = 9.91e-14

The two agree here because the input is well-scaled and double precision is generous, but the single $(b, a)$ path is the fragile one: its 6th-order denominator has tightly clustered poles whose locations are hypersensitive to coefficient rounding. In fixed point, or at higher order, lfilter on the flat polynomial degrades while the three-section sosfilt cascade stays accurate. Trust the SOS form.

Exercise 4: Transposition preserves the transfer function

The transposed Direct Form II second-order section has the state equations \[y[n] = b_0\,x[n] + s_1[n-1], \quad s_1[n] = b_1\,x[n] - a_1\,y[n] + s_2[n-1], \quad s_2[n] = b_2\,x[n] - a_2\,y[n]\] (this is the structure you derive by hand in Exercise 5). Show by direct calculation, taking $z$-transforms and eliminating the internal states $S_1(z)$ and $S_2(z)$, that it has the same transfer function as the direct form it was transposed from.

Answer

Take the $z$-transform of each equation, where a one-sample delay becomes a factor $z^{-1}$: \[Y(z) = b_0 X(z) + z^{-1} S_1(z), \quad S_1(z) = b_1 X(z) - a_1 Y(z) + z^{-1} S_2(z), \quad S_2(z) = b_2 X(z) - a_2 Y(z).\] Substitute $S_2(z)$ into the expression for $S_1(z)$: \[S_1(z) = (b_1 + b_2 z^{-1})\,X(z) - (a_1 + a_2 z^{-1})\,Y(z),\] then substitute that into the equation for $Y(z)$: \[Y(z) = b_0 X(z) + z^{-1}\bigl[(b_1 + b_2 z^{-1})\,X(z) - (a_1 + a_2 z^{-1})\,Y(z)\bigr].\] Collecting the $Y$ terms on the left, \[Y(z)\,\bigl(1 + a_1 z^{-1} + a_2 z^{-2}\bigr) = \bigl(b_0 + b_1 z^{-1} + b_2 z^{-2}\bigr)\,X(z),\] so \[H(z) = \frac{Y(z)}{X(z)} = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}},\] which is exactly the transfer function of the Direct Form II section it was transposed from. The transpose rearranges where the additions happen, not what is computed: same $H(z)$, different numerical behaviour. That is why the transposition theorem is “free.”

Exercise 5: Transpose Direct Form II by hand

Starting from a second-order Direct Form II section, apply the three transposition rules (reverse every branch, swap branch nodes with adders, exchange input and output) and describe the resulting transposed Direct Form II. Where do the delays and the multiplies by $b_k$ and $-a_k$ end up?

Answer

Reversing the graph turns the input fan-out node (where $x[n]$ splits into the $b_k$ taps) into the output summing node, and turns the output adder into the input fan-out. The two delay elements stay, but now each one sits between an adder and the next, accumulating a partial sum rather than holding a raw sample. The result is the structure \[y[n] = b_0\,x[n] + s_1[n-1], \quad s_1[n] = b_1\,x[n] - a_1\,y[n] + s_2[n-1], \quad s_2[n] = b_2\,x[n] - a_2\,y[n].\] Each $b_k$ multiplies the input and each $-a_k$ multiplies the output, with the products injected into the delay chain. This is transposed Direct Form II, the form sosfilt and CMSIS df2T use, and it is preferred because the intermediate states $s_1, s_2$ stay smaller than the raw accumulator of plain DF-II.

Summary

The same $H(z)$ has many structures (realizations). In exact arithmetic they are identical; in finite precision they differ in memory, operation count, and robustness.
Direct forms implement the difference equation literally. DF-II uses the minimum number of delays; the transposed DF-II has the best numerical behaviour and is the floating-point default.
Transposition turns any structure into a twin with the same $H(z)$ but different arithmetic, for free.
Cascade (product of biquads) and parallel (sum of sections) forms break a fragile high-order polynomial into robust low-order sections. Cascade SOS is the production standard.
Lattice structures realize the filter as a ladder of reflection-coefficient stages, with a trivial $|k_m| < 1$ stability test, favoured in speech and adaptive filtering.
The practical rule: factor into second-order sections and run them in transposed Direct Form II.

You now have the full design toolchain: sampling and quantization, discrete-time systems, the z-domain and frequency domain, filter design, and the structures that turn a design into something a processor can run. The chapters that follow put it to work on real signals.

Going deeper

The next chapters apply these structures directly: biquad filters builds cascaded second-order sections for audio equalizers and shows them on real microcontrollers, and multirate systems uses polyphase and CIC structures for efficient sample-rate conversion. In the Topics section, adaptive filtering lets a filter’s coefficients adjust themselves in real time. And the capstone exercise ties the first six chapters together in one end-to-end pipeline; this chapter’s structures are how its final design would actually be run.

References

For the transposition theorem, the cascade and parallel forms, and the lattice (including the step-down recursion), see Oppenheim and Schafer (Oppenheim and Schafer 2010, ch. 6), Proakis and Manolakis (Proakis and Manolakis 2007), and Mitra (Mitra 2006). The cascade-of-biquads argument is developed further on the biquad page, and the numerical-stability motivation appears in the filter design chapter.

Mitra, Sanjit K. 2006. Digital Signal Processing: A Computer-Based Approach. 3rd ed. McGraw-Hill.

Oppenheim, Alan V., and Ronald W. Schafer. 2010. Discrete-Time Signal Processing. 3rd ed. Pearson.

Proakis, John G., and Dimitris G. Manolakis. 2007. Digital Signal Processing: Principles, Algorithms, and Applications. 4th ed. Pearson.

--- title: "Filter Structures" subtitle: "Realizing a transfer function: direct, transposed, cascade, parallel, lattice" description-meta: "Realizing a transfer function in practice: direct, transposed, cascade, parallel and lattice structures, and what finite-precision arithmetic does to each." --- ```{python} #| echo: false import numpy as np from scipy.signal import lfilter, tf2sos, sosfilt, residuez, butter ``` [Chapter 6](06-filter-design.qmd) ended with a set of coefficients: a transfer function $H(z)$ that meets a specification. That is the *what*. This chapter is the *how*: given $H(z)$, how do you actually wire up the delays, multipliers, and adders that compute it sample by sample? The surprise is that there is no single answer. The same $H(z)$ can be realized by many different block diagrams, all mathematically identical in exact arithmetic, yet very different once you build them. They differ in how much memory they need, how many operations per sample, and (most importantly) how gracefully they degrade when the coefficients and signals are stored in finite-precision numbers. Choosing a **structure** is the last design decision before the filter becomes code or silicon. It is the real start of the "math to metal" pipeline. We work with the general difference equation from [Chapter 2](02-discrete-time.qmd), $$y[n] = \sum_{k=0}^{M} b_k\, x[n-k] \;-\; \sum_{k=1}^{N} a_k\, y[n-k],$$ whose transfer function is $$H(z) = \frac{\sum_{k=0}^{M} b_k\, z^{-k}}{1 + \sum_{k=1}^{N} a_k\, z^{-k}} = \frac{B(z)}{A(z)}.$$ Every structure in this chapter realizes this same $H(z)$. They are different *circuits* for one *function*. <hr> ## Why structure matters In exact arithmetic, rearranging a signal-flow graph that computes $H(z)$ cannot change its output. Real hardware is not exact. Coefficients are rounded to a finite word length, and intermediate sums are quantized. Two effects make the choice of structure consequential: - **Coefficient sensitivity.** The poles of a high-order $A(z)$ can be extremely sensitive to small changes in the polynomial coefficients. Rounding a single $a_k$ by one least-significant bit can move a pole enough to ripple the passband or even push it outside the unit circle, making a "stable" design unstable. Breaking the filter into small sections makes each section's poles far less sensitive. - **Dynamic range and limit cycles.** Different structures produce different intermediate signal magnitudes for the same input. A structure whose internal nodes swing widely is prone to overflow in fixed point; feedback combined with quantization can also trap an IIR filter in a small self-sustaining oscillation called a limit cycle. So the questions a structure answers are practical: How many delay registers? How many multiply-accumulates per sample? How robust is it to coefficient rounding? Does it overflow? The rest of the chapter introduces the structures in roughly increasing order of numerical robustness. <hr> ## Direct forms The most literal realization writes the difference equation out term by term, one multiplier per coefficient and one delay per $z^{-1}$. This is the **direct form**. [Chapter 2](02-discrete-time.qmd#direct-form-structures) already drew the generalized FIR transversal structure and the IIR **Direct Form I** (DF-I), which keeps separate delay lines for the input history and the output history. DF-I uses $M + N$ delay elements. You can do better: because the all-zero section $B(z)$ and the all-pole section $1/A(z)$ are in series, and series LTI systems commute, you can swap their order and then *share* one delay line between them. That is **Direct Form II** (DF-II), the canonical form, which uses only $\max(M, N)$ delays, the theoretical minimum. The [biquad page](09-biquad/index.qmd#filter-structures) works through DF-I, DF-II, and the transposed form for the second-order case in detail, with diagrams. Here we focus on what is *general* about structures, starting with the transformation that turns one structure into another. <hr> ## The transposition theorem Every structure has a twin that realizes the **same** transfer function. You obtain it by **transposing** the signal-flow graph, following three steps (the transposition theorem): 1. **Reverse the direction of every branch.** 2. **Replace every branch (pickoff) node with an adder, and every adder with a branch node.** 3. **Interchange the input and the output.** Reversing a flow graph leaves its transfer function unchanged, so the transposed structure computes exactly the same $H(z)$ with the same coefficients, but it distributes the arithmetic differently and therefore has different numerical behaviour. Applied to DF-II, this gives the **transposed Direct Form II** (TDF-II). It is the workhorse of practical IIR filtering: its state update accumulates smaller intermediate values than DF-II, it pipelines naturally, and it is the form that SciPy's `sosfilt` and ARM's CMSIS-DSP `arm_biquad_cascade_df2T` use internally. The [biquad page shows its diagram and state equations](09-biquad/index.qmd#direct-form-ii-transposed). ::: {.callout-note title="The point of transposition"} Transposition is free: it costs nothing to apply and never changes $H(z)$. It is purely a tool for trading numerical properties. When someone says "use the transposed form," they mean "same filter, friendlier arithmetic." ::: <hr> ## Cascade form A single high-order $(b, a)$ pair is the most fragile way to implement a filter, precisely because of the coefficient sensitivity described above. The standard fix is to **factor** the transfer function into a product of low-order sections: $$H(z) = g \prod_{k=1}^{L} H_k(z), \qquad H_k(z) = \frac{b_{0k} + b_{1k}\,z^{-1} + b_{2k}\,z^{-2}}{1 + a_{1k}\,z^{-1} + a_{2k}\,z^{-2}}.$$ Each $H_k(z)$ is a second-order section (a **biquad**), and the sections run in series, the output of one feeding the next. ```{dot} //| echo: false //| label: fig-cascade //| fig-cap: "Cascade form. The transfer function is factored into a series of low-order sections; in practice each H_k(z) is a biquad. " graph { layout=neato node [fontname="sans-serif" fontsize=14] edge [arrowsize=0.7 dir=forward] x [label="x[n]" shape=plaintext fontname="serif" pos="0,0!"] h1 [label=<H1(z)> shape=box width=0.95 height=0.62 pos="1.4,0!"] h2 [label=<H2(z)> shape=box width=0.95 height=0.62 pos="2.9,0!"] dd [label="···" shape=plaintext pos="4.1,0!"] hL [label=<HL(z)> shape=box width=0.95 height=0.62 pos="5.3,0!"] y [label="y[n]" shape=plaintext fontname="serif" pos="6.7,0!"] x--h1 h1--h2 h2--dd dd--hL hL--y } ``` An order-$N$ filter needs $\lceil N/2 \rceil$ sections. The win is robustness: a biquad's two poles are easy to place accurately, and rounding its handful of coefficients barely moves them, whereas rounding the coefficients of one big polynomial can move its poles dramatically. Cascading also makes stability trivial to check section by section. In Python, `scipy.signal.tf2sos` performs the factorization and pole-zero pairing for you, and `sosfilt` runs the cascade (each section in transposed Direct Form II): ```{python} # A 4th-order lowpass, expressed two ways b, a = butter(4, 0.25) # single (b, a) pair sos = tf2sos(b, a) # cascade of second-order sections print("sections:", sos.shape[0]) print(np.round(sos, 4)) # The two realizations agree to machine precision (at this modest order) rng = np.random.default_rng(0) x = rng.standard_normal(2000) err = np.max(np.abs(sosfilt(sos, x) - lfilter(b, a, x))) print(f"max |sosfilt - lfilter| = {err:.2e}") ``` Each row of the `sos` matrix is one section, laid out as $[\,b_0, b_1, b_2, 1, a_1, a_2\,]$. At fourth order both realizations agree, but as the order grows the single-polynomial path loses accuracy fast while the cascade stays clean. The [filter design chapter](06-filter-design.qmd#second-order-sections) and the [biquad page](09-biquad/index.qmd#cascading-biquads) show that divergence directly. This is why every serious IIR implementation ships as second-order sections. ::: {.callout-tip title="Section ordering matters"} The factorization is not unique: you can pair any pole with any zero and order the sections however you like. These choices do not change $H(z)$ but they do change the internal signal levels, which matters for overflow in fixed point. `tf2sos` uses a sensible default (pair each pole with its nearest zero, order to keep intermediate gains moderate). For most work the default is fine. ::: <hr> ## Parallel form Instead of *factoring* $H(z)$ into a product, you can *expand* it into a sum. A **partial-fraction expansion** writes the transfer function as $$H(z) = \sum_{k} H_k(z) \;+\; (\text{optional direct term}),$$ where each $H_k(z)$ is again a low-order (typically second-order) section, but now the sections run in **parallel**: the input feeds all of them at once and their outputs are summed. ```{dot} //| echo: false //| label: fig-parallel //| fig-cap: "Parallel form. The input drives every section simultaneously and the section outputs are summed. " graph { layout=neato node [fontname="sans-serif" fontsize=14] edge [arrowsize=0.65 dir=forward] x [label="x[n]" shape=plaintext fontname="serif" pos="-0.7,1.5!"] bt [label="" shape=none width=0 height=0 margin=0 pos="0.4,3!"] bm [label="" shape=point width=0.05 pos="0.4,1.5!"] bd [label="⋮" shape=plaintext pos="0.4,0.45!"] bb [label="" shape=none width=0 height=0 margin=0 pos="0.4,-0.6!"] h0 [label=<H1(z)> shape=box width=1.0 height=0.55 pos="1.7,3!"] h1 [label=<H2(z)> shape=box width=1.0 height=0.55 pos="1.7,1.5!"] hd [label="⋮" shape=plaintext pos="1.7,0.45!"] hL [label=<HL(z)> shape=box width=1.0 height=0.55 pos="1.7,-0.6!"] c0 [label="+" shape=circle width=0.34 fixedsize=true pos="3.2,3!"] c1 [label="+" shape=circle width=0.34 fixedsize=true pos="3.2,1.5!"] cd [label="⋮" shape=plaintext pos="3.2,0.45!"] cL [label="" shape=none width=0 height=0 margin=0 pos="3.2,-0.6!"] y [label="y[n]" shape=plaintext fontname="serif" pos="4.3,3!"] x--bm [arrowhead=none] bt--bm [arrowhead=none] bm--bd [arrowhead=none] bd--bb [arrowhead=none] bt--h0 bm--h1 bb--hL h0--c0 h1--c1 hL--cL [arrowhead=none] cL--cd [arrowhead=none] cd--c1 c1--c0 c0--y } ``` `scipy.signal.residuez` computes the expansion (residues $r_k$, poles $p_k$, and any direct terms). Grouping each complex-conjugate pole pair into one real second-order section gives a realizable parallel structure: ```{python} def parallel_sections(b, a): """Partial-fraction expansion grouped into real second-order sections.""" r, p, k = residuez(b, a) sections, used = [], np.zeros(len(p), bool) for i in range(len(p)): if used[i]: continue # find the conjugate partner of a complex pole j = next((jj for jj in range(len(p)) if jj != i and not used[jj] and np.isclose(p[jj], np.conj(p[i])) and abs(p[i].imag) > 1e-9), -1) if j >= 0: # conjugate pair -> real biquad used[i] = used[j] = True bsec = np.real([r[i] + r[j], -(r[i]*p[j] + r[j]*p[i]), 0.0]) asec = np.real([1.0, -(p[i] + p[j]), p[i]*p[j]]) else: # lone real pole -> first order used[i] = True bsec, asec = np.real([r[i], 0.0, 0.0]), np.real([1.0, -p[i], 0.0]) sections.append((bsec, asec)) return sections, np.real(k) b, a = butter(4, 0.25) sections, direct = parallel_sections(b, a) y_par = sum(lfilter(bs, as_, x) for bs, as_ in sections) if len(direct): y_par = y_par + lfilter(direct, [1.0], x) print(f"parallel sections: {len(sections)}") print(f"max |parallel - lfilter| = {np.max(np.abs(y_par - lfilter(b, a, x))):.2e}") ``` The cascade and parallel forms realize the same filter and share the same robustness advantage over a single high-order section. They differ in detail: the cascade preserves the exact zeros of $H(z)$ in each section (good when the zeros matter, as in notch filters), while the parallel form keeps the sections fully independent (good for parallel hardware, and the natural form when you think of a filter as a sum of resonances). <hr> ## Lattice structures The **lattice** is a different shape entirely. Instead of a delay line tapped by coefficients, it is a ladder of identical two-input, two-output stages, each parametrized by a single **reflection coefficient** $k_m$. Lattices are prized in speech coding (linear predictive coding), adaptive filtering, and anywhere coefficient robustness and a cheap stability test matter. Each stage carries a **forward** signal $f_m[n]$ and a **backward** signal $g_m[n]$. For an all-zero (FIR) filter the stages run from $m = 1$ to $M$: $$f_0[n] = g_0[n] = x[n],$$ $$f_m[n] = f_{m-1}[n] + k_m\, g_{m-1}[n-1],$$ $$g_m[n] = k_m\, f_{m-1}[n] + g_{m-1}[n-1],$$ and the output is $y[n] = f_M[n]$. A single stage is the cross-coupled "butterfly" below: the forward signal passes straight through and also feeds the backward adder through $k_m$, while the delayed backward signal passes straight through and feeds the forward adder through $k_m$. ```{dot} //| echo: false //| label: fig-fir-lattice //| fig-cap: "One stage of an FIR (all-zero) lattice. The forward paths $f_m$ and backward paths $g_m$ are cross-coupled through the reflection coefficient $k_m$; the delay sits on the backward path." graph { layout=neato node [fontname="sans-serif" fontsize=13] edge [arrowsize=0.6 dir=forward] fin [label=<fm−1> shape=plaintext fontname="serif" pos="-0.5,2.0!"] fb [label="" shape=point width=0.05 pos="0.7,2.0!"] At [label="+" shape=circle width=0.34 fixedsize=true pos="3.3,2.0!"] fo [label=<fm> shape=plaintext fontname="serif" pos="4.3,2.0!"] gin [label=<gm−1> shape=plaintext fontname="serif" pos="-0.5,0!"] z [label="z⁻¹" shape=box width=0.5 height=0.38 pos="0.7,0!"] gb [label="" shape=point width=0.05 pos="1.6,0!"] Bb [label="+" shape=circle width=0.34 fixedsize=true pos="3.3,0!"] go [label=<gm> shape=plaintext fontname="serif" pos="4.3,0!"] kd [label=<km> shape=triangle orientation=232 width=0.56 height=0.5 fixedsize=true pos="1.55,1.35!"] ku [label=<km> shape=triangle orientation=320 width=0.56 height=0.5 fixedsize=true pos="1.95,0.41!"] fin--fb [arrowhead=none] fb--At At--fo gin--z z--gb [arrowhead=none] gb--Bb Bb--go fb--kd [arrowhead=none] kd--Bb gb--ku [arrowhead=none] ku--At } ``` The reflection coefficients are not the filter coefficients $b_k$; you obtain them from the polynomial by a **step-down (Levinson) recursion**. Starting from the monic polynomial $A(z) = 1 + a_1 z^{-1} + \cdots + a_M z^{-M}$, peel off one stage at a time: $$k_m = a_m^{(m)}, \qquad a_i^{(m-1)} = \frac{a_i^{(m)} - k_m\, a_{m-i}^{(m)}}{1 - k_m^2}.$$ The same lattice stages, rewired so the forward recursion runs backward, realize the **all-pole** filter $1/A(z)$: $$f_M[n] = x[n], \qquad f_{m-1}[n] = f_m[n] - k_m\, g_{m-1}[n-1], \qquad g_m[n] = k_m\, f_{m-1}[n] + g_{m-1}[n-1],$$ with output $y[n] = f_0[n]$. The code below builds both lattices from their reflection coefficients and confirms they reproduce the ordinary direct-form output to machine precision (the near-zero errors *are* the verification; this chapter introduces lattices, which `scipy.signal` does not provide ready-made): ```{python} def poly_to_reflection(a): """Step-down (Levinson) recursion: monic polynomial -> reflection coefficients.""" a = np.asarray(a, float) a = a / a[0] M = len(a) - 1 k = np.zeros(M) for m in range(M, 0, -1): k[m - 1] = a[m] a = (a[:m + 1] - a[m] * a[:m + 1][::-1])[:m] / (1 - a[m] ** 2) return k def fir_lattice(k, x): """All-zero FIR lattice; returns f_M[n].""" M = len(k) g_delayed = np.zeros(M) # g_{m-1}[n-1] y = np.zeros_like(x, float) for n in range(len(x)): f, g = x[n], np.empty(M + 1) g[0] = x[n] for m in range(1, M + 1): f_next = f + k[m - 1] * g_delayed[m - 1] g[m] = k[m - 1] * f + g_delayed[m - 1] f = f_next y[n] = f g_delayed = g[:M].copy() return y def allpole_lattice(k, x): """All-pole IIR lattice realizing 1/A(z); returns f_0[n].""" M = len(k) g_delayed = np.zeros(M) y = np.zeros_like(x, float) for n in range(len(x)): f, g = x[n], np.empty(M + 1) for m in range(M, 0, -1): f = f - k[m - 1] * g_delayed[m - 1] g[m] = k[m - 1] * f + g_delayed[m - 1] g[0] = f y[n] = f g_delayed = g[:M].copy() return y # FIR lattice vs direct-form FIR b_fir = np.array([1.0, -0.6, 0.3, 0.2]) k_fir = poly_to_reflection(b_fir) err_fir = np.max(np.abs(b_fir[0] * fir_lattice(k_fir, x) - lfilter(b_fir, [1.0], x))) print(f"FIR reflection coeffs = {np.round(k_fir, 4)}") print(f"FIR max |lattice - lfilter| = {err_fir:.2e}") # All-pole lattice vs direct-form IIR a_iir = np.array([1.0, -0.5, 0.25]) k_iir = poly_to_reflection(a_iir) err_iir = np.max(np.abs(allpole_lattice(k_iir, x) - lfilter([1.0], a_iir, x))) print(f"IIR reflection coeffs = {np.round(k_iir, 4)}") print(f"IIR max |lattice - lfilter| = {err_iir:.2e}") ``` Two properties make the lattice attractive on hardware. First, **stability is trivial to test**: an all-pole lattice is stable if and only if every reflection coefficient satisfies $|k_m| < 1$, a single comparison per stage rather than a root-finding problem. Second, the lattice is **modular and well-conditioned**: increasing the order just appends a stage without recomputing the earlier ones, and the structure tolerates coefficient rounding gracefully, which is why it dominates speech LPC. ::: {.callout-note title="From math to metal"} ARM's CMSIS-DSP library ships lattice filters ready to drop onto a Cortex-M: `arm_fir_lattice_f32` and `arm_iir_lattice_f32` (with `q15`/`q31` fixed-point variants). The FIR function takes exactly the reflection coefficients computed above; the IIR function realizes a general pole-zero response, so it needs the reflection coefficients *plus* a set of ladder (tap) coefficients for the numerator, with the all-pole case here as its special instance. The transposed Direct Form II biquad cascade from earlier is `arm_biquad_cascade_df2T_f32`. Picking a structure here *is* picking the library function you will call on the metal. ::: <hr> ## Transposition on an FPGA: where the pipeline registers go The transposition theorem said the two forms compute the same $H(z)$ with the same coefficients, and that the choice is about arithmetic rather than mathematics. On an FPGA that choice has a name: the **critical path**, the longest run of combinational logic between two registers, which is what sets the maximum clock frequency. Take a direct-form FIR. Every tap's product has to be summed *before* the output register, so the adders form a chain (or a tree) between the multipliers and the output: the critical path grows with the number of taps, and a 64-tap filter clocks slower than an 8-tap one. Transpose it and the delays land *between* the adders instead of ahead of them. Each stage becomes one multiplier and one adder feeding its own register, so the longest path is a single multiply plus a single add **whatever the tap count is**. A 64-tap transposed FIR clocks at the same rate as an 8-tap one; it just uses more of the chip. ### VHDL: transposed direct-form FIR ```vhdl library ieee; use ieee.std_logic_1164.all; use ieee.numeric_std.all; -- Transposed direct-form FIR. One register after every adder, so the longest -- combinational path is ONE multiply plus ONE add regardless of N. entity fir_transposed is generic ( N : positive := 8; -- number of taps B_IN : positive := 16; -- input word width B_C : positive := 16; -- coefficient width B_ACC : positive := 34 -- = B_IN + B_C - 1 + ceil(log2 N); never overflows ); port ( clk : in std_logic; rst : in std_logic; -- synchronous, active high x_in : in signed(B_IN-1 downto 0); -- one sample per clk y_out : out signed(B_ACC-1 downto 0) -- full-precision output ); end entity; architecture rtl of fir_transposed is type coeff_array is array (0 to N-1) of signed(B_C-1 downto 0); -- Symmetric low-pass, for illustration. Any tap set works unchanged, but -- this aggregate has to be restated if you change N. constant COEFF : coeff_array := ( to_signed( -412, B_C), to_signed( -906, B_C), to_signed( 1247, B_C), to_signed( 9058, B_C), to_signed( 9058, B_C), to_signed( 1247, B_C), to_signed( -906, B_C), to_signed( -412, B_C)); type acc_array is array (0 to N-1) of signed(B_ACC-1 downto 0); signal acc : acc_array := (others => (others => '0')); begin process(clk) begin if rising_edge(clk) then if rst = '1' then acc <= (others => (others => '0')); else -- Last stage has nothing upstream of it: it starts a fresh partial sum. acc(N-1) <= resize(x_in * COEFF(N-1), B_ACC); -- Every other stage adds its own product to the partial sum arriving -- from the stage to its right. acc(k+1) reads its PRE-EDGE value, which -- is exactly the pipeline register that makes the path short. for k in N-2 downto 0 loop acc(k) <= acc(k+1) + resize(x_in * COEFF(k), B_ACC); end loop; end if; end if; end process; y_out <= acc(0); end architecture; ``` Note what is *not* here: no delay line. The direct form needs an $N$-deep shift register holding past inputs; the transposed form holds partial sums instead, and the same sample is broadcast to every multiplier at once. The storage is wider (`B_ACC` rather than `B_IN` per stage) but the structure is otherwise identical in cost, which is the hardware statement of "transposition is free". Size `B_ACC` as $B_\text{in} + B_c - 1 + \lceil \log_2 N \rceil$ and no stage can overflow: each product needs $B_\text{in} + B_c - 1$ bits and summing $N$ of them adds $\lceil \log_2 N \rceil$. The model below checks that bound and the transposition claim at the same time, bit for bit. ```{python} #| label: fir-transposed-model #| code-fold: false import numpy as np N, B_IN, B_C = 8, 16, 16 B_ACC = B_IN + B_C - 1 + int(np.ceil(np.log2(N))) # = 34 bits COEFF = np.array([-412, -906, 1247, 9058, 9058, 1247, -906, -412]) def fir_transposed(x): """Bit-accurate model of the VHDL entity: N accumulator registers, one clock per sample.""" acc = [0] * N out = [] for xn in x: nxt = [0] * N nxt[N - 1] = int(xn) * int(COEFF[N - 1]) # fresh partial sum for k in range(N - 2, -1, -1): nxt[k] = acc[k + 1] + int(xn) * int(COEFF[k]) # pre-edge read acc = nxt out.append(acc[0]) return np.array(out, dtype=object), acc rng = np.random.default_rng(0) x = rng.integers(-2**(B_IN-1), 2**(B_IN-1), size=500) y_hw, _ = fir_transposed(x) y_ref = np.convolve(x, COEFF)[:len(x)] # the direct form, exactly print(f"B_ACC needed = {B_ACC} bits, largest |acc| seen = {int(max(abs(y_hw)))}") print(f"transposed == direct : {np.array_equal(y_hw.astype(np.int64), y_ref)}") # Transposition changes the arithmetic, not the answer: bit-exact, not merely close. assert np.array_equal(y_hw.astype(np.int64), y_ref), "transposed form must equal the direct form" # And the width bound must hold with room to spare for the worst case, not just this input. worst = N * 2**(B_IN-1) * 2**(B_C-1) assert worst <= 2**(B_ACC-1), "B_ACC must cover the worst-case sum" assert int(max(abs(y_hw))) < 2**(B_ACC-1), "this run must stay inside the accumulator" ``` ### Testbench: check it against the direct form The Python above proves the equivalence outside the simulator. Here is the same check *inside* it, which is the version you can run. The reference model is a direct-form FIR written the obvious way -- keep the last $N$ inputs, multiply and sum them all in one clock -- so the testbench compares the two structures from the transposition theorem against each other, cycle by cycle, and fails on the first sample where they part. Stimulus is the [Galois LFSR](../topics/noise-generation/embedded.qmd#the-lfsr-on-an-fpga-16-flip-flops-and-four-xor-gates) from the noise generation page, which keeps the whole thing self-contained: no vector files, no scripts, nothing outside these two entities. ```vhdl library ieee; use ieee.std_logic_1164.all; use ieee.numeric_std.all; use std.env.all; -- VHDL-2008, for stop entity tb_fir_transposed is end entity; architecture sim of tb_fir_transposed is constant N : positive := 8; constant B_IN : positive := 16; constant B_C : positive := 16; constant B_ACC : positive := 34; constant N_SAMPLES : positive := 2000; type coeff_array is array (0 to N-1) of signed(B_C-1 downto 0); constant COEFF : coeff_array := ( to_signed( -412, B_C), to_signed( -906, B_C), to_signed( 1247, B_C), to_signed( 9058, B_C), to_signed( 9058, B_C), to_signed( 1247, B_C), to_signed( -906, B_C), to_signed( -412, B_C)); signal clk : std_logic := '0'; signal rst : std_logic := '1'; signal x_in : signed(B_IN-1 downto 0) := (others => '0'); signal y_dut : signed(B_ACC-1 downto 0); signal y_ref : signed(B_ACC-1 downto 0) := (others => '0'); begin clk <= not clk after 5 ns; -- 100 MHz; the rate is arbitrary here dut : entity work.fir_transposed generic map (N => N, B_IN => B_IN, B_C => B_C, B_ACC => B_ACC) port map (clk => clk, rst => rst, x_in => x_in, y_out => y_dut); -- Direct form: the long adder chain, all N products summed in one clock. -- Correct, and exactly what the transposed form must reproduce. reference : process(clk) type sample_array is array (0 to N-2) of signed(B_IN-1 downto 0); variable dl : sample_array := (others => (others => '0')); variable acc : signed(B_ACC-1 downto 0); begin if rising_edge(clk) then if rst = '1' then dl := (others => (others => '0')); y_ref <= (others => '0'); else -- Sum from the CURRENT sample and the delay line as it stands before -- this edge's shift, so that y_ref carries the same one-clock latency -- as the DUT. acc := resize(x_in * COEFF(0), B_ACC); for k in 1 to N-1 loop acc := acc + resize(dl(k-1) * COEFF(k), B_ACC); end loop; y_ref <= acc; for k in N-2 downto 1 loop -- then shift, newest at 0 dl(k) := dl(k-1); end loop; dl(0) := x_in; end if; end if; end process; -- Both structures reset to all zeros, so they must agree from the very first -- clock after reset: there is no priming window to skip past. checker : process(clk) variable cycle : natural := 0; begin if rising_edge(clk) and rst = '0' then cycle := cycle + 1; assert y_dut = y_ref report "MISMATCH at cycle " & integer'image(cycle) & ": transposed /= direct" severity failure; end if; end process; stimulus : process variable lfsr : std_logic_vector(15 downto 0) := x"ACE1"; begin wait until rising_edge(clk); wait until rising_edge(clk); rst <= '0'; for i in 1 to N_SAMPLES loop wait until rising_edge(clk); if lfsr(0) = '1' then -- Galois step, x"B400" lfsr := ('0' & lfsr(15 downto 1)) xor x"B400"; else lfsr := '0' & lfsr(15 downto 1); end if; x_in <= signed(lfsr); -- full-scale, exercises B_ACC end loop; wait until rising_edge(clk); report "PASS: " & integer'image(N_SAMPLES) & " samples, transposed = direct on every cycle" severity note; stop; end process; end architecture; ``` Two lines to run it, assuming both files are in the working directory: ```bash ghdl -a --std=08 fir_transposed.vhd tb_fir_transposed.vhd ghdl -r --std=08 tb_fir_transposed ``` A pass prints the `note`; a failure names the cycle it first diverged on and stops there. The `severity failure` matters: without it a mismatch is a printed line the run scrolls past, and a testbench you have to read the output of is a testbench that will eventually be believed while broken. The same harness works for the other three entities on this site -- the [CIC decimator](10-multirate/embedded.qmd#cic-on-an-fpga-registers-and-adders), the [moving average](08-smoothing/embedded.qmd#moving-average-on-an-fpga-delay-line-accumulator) and the [LFSR](../topics/noise-generation/embedded.qmd#the-lfsr-on-an-fpga-16-flip-flops-and-four-xor-gates) -- with the reference process swapped for whatever independent computation of the same quantity you trust. For the LFSR that reference is not a second implementation at all but the period itself: assert that the state returns to its seed after exactly $2^{16}-1$ clocks and never before. ::: {.callout-note title="Not simulated"} The entity and its testbench are faithful transcriptions of the structure checked in Python, but neither has been through a simulator here -- **including the testbench itself, which has not been run through `ghdl`** -- so both are verified by inspection and against the model above. Run the two commands above before trusting any of it on hardware. ::: <hr> ## Choosing a structure All of these realize the same $H(z)$. The differences are practical: | Structure | Delays | Robustness to coefficient rounding | Typical use | |-----------|:------:|------------------------------------|-------------| | **Direct Form I** | $M + N$ | Poor at high order; good per low-order section | Fixed-point biquads (separate delay lines resist overflow) | | **Direct Form II** | $\max(M,N)$ | Poor at high order; wide internal range | Rarely used directly at high order | | **Transposed DF-II** | $\max(M,N)$ | Good per section; small internal values | Default for floating-point (SciPy `sosfilt`, CMSIS `df2T`) | | **Cascade (SOS)** | $2\lceil N/2\rceil$ | Very good (low-order sections, exact zeros) | The standard production IIR form | | **Parallel** | $2\lceil N/2\rceil$ | Very good (independent sections) | Parallel hardware; sum-of-resonances designs | | **Lattice** | $M$ or $N$ | Very good; trivial $\lvert k_m\rvert < 1$ stability test | Speech LPC, adaptive filters | The headline advice is simple. Never implement a high-order IIR filter as a single $(b, a)$ pair. Factor it into **second-order sections** and run them in **transposed Direct Form II**, which is exactly what `scipy.signal.sosfilt` and CMSIS-DSP give you by default. Reach for **lattice** structures when you need a cheap stability guarantee or you are doing linear prediction. <hr> ## Exercises ::: {.callout-tip collapse="true" title="Exercise 1: Read a structure"} A first-order structure is described by the difference equation $$y[n] = b_0\, x[n] + b_1\, x[n-1] - a_1\, y[n-1].$$ Write its transfer function $H(z) = Y(z)/X(z)$. Is this Direct Form I or Direct Form II if the input and output histories use separate delay registers? ::: {.callout-note collapse="true" title="Answer"} Taking the $z$-transform of both sides and collecting terms, $$H(z) = \frac{b_0 + b_1 z^{-1}}{1 + a_1 z^{-1}}.$$ Separate delay registers for the input ($x[n-1]$) and the output ($y[n-1]$) make this **Direct Form I**. Sharing a single register between the feedforward and feedback paths would make it Direct Form II. ::: ::: ::: {.callout-tip collapse="true" title="Exercise 2: Biquad difference equation and resonance"} A biquad section has $$H(z) = \frac{1 + 0.1\,z^{-1} - 0.7199\,z^{-2}}{1 - 1.1786\,z^{-1} + 0.7246\,z^{-2}}.$$ (a) Write its Direct Form II difference equations. (b) With $f_s = 100$ kHz, find the resonance frequency from the pole angle. ::: {.callout-note collapse="true" title="Answer"} **(a)** With the shared state $w[n]$ and the denominator coefficients $a_1 = -1.1786$, $a_2 = 0.7246$: $$w[n] = x[n] + 1.1786\,w[n-1] - 0.7246\,w[n-2],$$ $$y[n] = w[n] + 0.1\,w[n-1] - 0.7199\,w[n-2].$$ **(b)** Write the denominator as $A(z) = 1 - 2R\cos\theta\,z^{-1} + R^2 z^{-2}$. Then $R = \sqrt{0.7246} = 0.851$ and $\cos\theta = 1.1786 / (2R) = 0.692$, so $\theta = 0.806$ rad. The resonance frequency is $$f_c = \frac{\theta}{2\pi}\, f_s = \frac{0.806}{2\pi}\,(100\ \text{kHz}) \approx 12.8\ \text{kHz}.$$ The poles have magnitude $R = 0.851 < 1$, so the section is stable. ::: ::: ::: {.callout-tip collapse="true" title="Exercise 3: Cascade decomposition in Python"} Take the 6th-order lowpass `b, a = butter(6, 0.2)`. Decompose it into second-order sections with `tf2sos`, and compare `sosfilt` against `lfilter` on a noise input. Which one do you trust at this order, and why? ::: {.callout-note collapse="true" title="Answer"} ```{python} b6, a6 = butter(6, 0.2) sos6 = tf2sos(b6, a6) xt = rng.standard_normal(4000) diff = np.max(np.abs(sosfilt(sos6, xt) - lfilter(b6, a6, xt))) print(f"{sos6.shape[0]} sections, max |sosfilt - lfilter| = {diff:.2e}") ``` The two agree here because the input is well-scaled and double precision is generous, but the single $(b, a)$ path is the fragile one: its 6th-order denominator has tightly clustered poles whose locations are hypersensitive to coefficient rounding. In fixed point, or at higher order, `lfilter` on the flat polynomial degrades while the three-section `sosfilt` cascade stays accurate. Trust the SOS form. ::: ::: ::: {.callout-tip collapse="true" title="Exercise 4: Transposition preserves the transfer function"} The transposed Direct Form II second-order section has the state equations $$y[n] = b_0\,x[n] + s_1[n-1], \quad s_1[n] = b_1\,x[n] - a_1\,y[n] + s_2[n-1], \quad s_2[n] = b_2\,x[n] - a_2\,y[n]$$ (this is the structure you derive by hand in Exercise 5). Show by direct calculation, taking $z$-transforms and eliminating the internal states $S_1(z)$ and $S_2(z)$, that it has the same transfer function as the direct form it was transposed from. ::: {.callout-note collapse="true" title="Answer"} Take the $z$-transform of each equation, where a one-sample delay becomes a factor $z^{-1}$: $$Y(z) = b_0 X(z) + z^{-1} S_1(z), \quad S_1(z) = b_1 X(z) - a_1 Y(z) + z^{-1} S_2(z), \quad S_2(z) = b_2 X(z) - a_2 Y(z).$$ Substitute $S_2(z)$ into the expression for $S_1(z)$: $$S_1(z) = (b_1 + b_2 z^{-1})\,X(z) - (a_1 + a_2 z^{-1})\,Y(z),$$ then substitute that into the equation for $Y(z)$: $$Y(z) = b_0 X(z) + z^{-1}\bigl[(b_1 + b_2 z^{-1})\,X(z) - (a_1 + a_2 z^{-1})\,Y(z)\bigr].$$ Collecting the $Y$ terms on the left, $$Y(z)\,\bigl(1 + a_1 z^{-1} + a_2 z^{-2}\bigr) = \bigl(b_0 + b_1 z^{-1} + b_2 z^{-2}\bigr)\,X(z),$$ so $$H(z) = \frac{Y(z)}{X(z)} = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}},$$ which is exactly the transfer function of the Direct Form II section it was transposed from. The transpose rearranges *where* the additions happen, not *what* is computed: same $H(z)$, different numerical behaviour. That is why the transposition theorem is "free." ::: ::: ::: {.callout-tip collapse="true" title="Exercise 5: Transpose Direct Form II by hand"} Starting from a second-order Direct Form II section, apply the three transposition rules (reverse every branch, swap branch nodes with adders, exchange input and output) and describe the resulting transposed Direct Form II. Where do the delays and the multiplies by $b_k$ and $-a_k$ end up? ::: {.callout-note collapse="true" title="Answer"} Reversing the graph turns the input fan-out node (where $x[n]$ splits into the $b_k$ taps) into the output summing node, and turns the output adder into the input fan-out. The two delay elements stay, but now each one sits *between* an adder and the next, accumulating a partial sum rather than holding a raw sample. The result is the structure $$y[n] = b_0\,x[n] + s_1[n-1], \quad s_1[n] = b_1\,x[n] - a_1\,y[n] + s_2[n-1], \quad s_2[n] = b_2\,x[n] - a_2\,y[n].$$ Each $b_k$ multiplies the input and each $-a_k$ multiplies the output, with the products injected into the delay chain. This is **transposed Direct Form II**, the form `sosfilt` and CMSIS `df2T` use, and it is preferred because the intermediate states $s_1, s_2$ stay smaller than the raw accumulator of plain DF-II. ::: ::: <hr> ## Summary - The same $H(z)$ has many **structures** (realizations). In exact arithmetic they are identical; in finite precision they differ in memory, operation count, and robustness. - **Direct forms** implement the difference equation literally. DF-II uses the minimum number of delays; the **transposed DF-II** has the best numerical behaviour and is the floating-point default. - **Transposition** turns any structure into a twin with the same $H(z)$ but different arithmetic, for free. - **Cascade** (product of biquads) and **parallel** (sum of sections) forms break a fragile high-order polynomial into robust low-order sections. Cascade SOS is the production standard. - **Lattice** structures realize the filter as a ladder of reflection-coefficient stages, with a trivial $|k_m| < 1$ stability test, favoured in speech and adaptive filtering. - The practical rule: factor into second-order sections and run them in transposed Direct Form II. You now have the full design toolchain: sampling and quantization, discrete-time systems, the z-domain and frequency domain, filter design, and the structures that turn a design into something a processor can run. The chapters that follow put it to work on real signals. ::: {.callout-note title="Going deeper"} The next chapters apply these structures directly: [biquad filters](09-biquad/index.qmd) builds cascaded second-order sections for audio equalizers and shows them on real microcontrollers, and [multirate systems](10-multirate/index.qmd) uses polyphase and CIC structures for efficient sample-rate conversion. In the [Topics](../topics/index.qmd) section, [adaptive filtering](../topics/adaptive-filtering/index.qmd) lets a filter's coefficients adjust themselves in real time. And the [capstone exercise](../exercises/capstone.qmd) ties the first six chapters together in one end-to-end pipeline; this chapter's structures are how its final design would actually be run. ::: ## References For the transposition theorem, the cascade and parallel forms, and the lattice (including the step-down recursion), see Oppenheim and Schafer [@oppenheim2010discrete, ch. 6], Proakis and Manolakis [@proakis2007digital], and Mitra [@mitra2006digital]. The cascade-of-biquads argument is developed further on the [biquad page](09-biquad/index.qmd), and the numerical-stability motivation appears in the [filter design chapter](06-filter-design.qmd#second-order-sections). ::: {#refs} :::