MZTi Matched Z-Transform improved lowpass example [by karrikuh, Matlab→Python translation]

readme.md

Translated from karrikuh's Matlab code: Struggling with MZTi lowpass

mzti.m.py

from numpy import pi, exp, tan, sqrt, log10
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt


def lp2_mzti(w0=None, K=None):
    """
    Digital 2nd order lowpass filter design using
    improved Matched Z transform (MZTi)

    Analog prototype is

                     1
    H(s) = -----------------------
           (s/w0)^2 + K*(s/w0) + 1

    Returns coeffs of biquadratic digital filter
    """

    # Pole mapping via conventional Matched Z transform
    aa = [1 / w0**2, K / w0, 1]
    p = np.roots(aa)  # analog poles
    a = signal.convolve([1, -exp(p[0])], [1, -exp(p[1])])

    # FIR correction filter
    w = pi * np.arange(3) / 3
    H = abs(signal.freqs([1], aa, w)[1] /
            signal.freqz([1], a, w)[1])  # relative magnitude error

    b1 = 0.5 * (H[0] - sqrt(H[0] ** 2 - 2 * H[1] ** 2 + 2 * H[2] ** 2))
    b2 = (3 * (H[0] - b1) -
          sqrt(-3 * H[0] ** 2 + 12 * H[1] ** 2 - 6 * H[0] * b1 - 3 * b1 ** 2)
          ) / 6
    b0 = H[0] - b1 - b2
    b = np.array([b0, b1, b2])

    return b, a


if __name__ == '__main__':
    w0 = 0.4 * pi
    Q = 4

    w = np.logspace(-1, 0, 1024) * pi

    # analog prototype
    Hs = signal.freqs([1], [1 / w0 ** 2, 1 / (w0 * Q), 1], w)[1]

    # bilinear Z transform (BLT)
    z = exp(-1j * w)
    s = 2 * (1 - z) / (1 + z) / (2 * tan(w0 / 2))
    Hblt = 1 / (s ** 2 + s / Q + 1)

    [b, a] = lp2_mzti(w0, 1 / Q)
    Hmzti = signal.freqz(b, a, w)[1]

    plt.semilogx(w / pi, 20*log10(abs(Hs)), label='Analog')
    plt.semilogx(w / pi, 20*log10(abs(Hblt)), 'g', label='BLT')
    plt.semilogx(w / pi, 20*log10(abs(Hmzti)), 'r', label='MZTi')
    plt.grid(True, which='both')
    plt.ylim([-40, 20])
    plt.legend()