Designing filter coefficients using SciPy

Here is an example of designing a Chebyshev-II type low-pass filter from the passband/stopband specifications using the cheb2ord() and cheby2() functions from the scipy.signal library. Please refer to the Scipy documentations (https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.cheb2ord.html) for more detailed descriptions.

#!/usr/bin/env xmodelpy

import scipy.signal
import numpy as np

# design of low-pass chebyshev-II filter
wp = 2*np.pi * 0.9e6    # passband edge frequency (rad/s)
ws = 2*np.pi * 1.1e6    # stopband edge frequency (rad/s)
rp = 3                  # maximum attenuation within pass band in dB
rs = 40                 # minimum attenuation within stop band in dB

N, wn = scipy.signal.cheb2ord(wp, ws, rp, rs, analog=True)
z,p,k = scipy.signal.cheby2(N, rs, wn, 'lowpass', analog=True, output='zpk')

And here is an example of designing a Chebyshev-II type band-pass filter:

# design of bass-pass chebyshev-II filter
wp = [2*np.pi*1.1e6, 2*np.pi*2.9e6]     # passband edge frequency (rad/s)
ws = [2*np.pi*0.9e6, 2*np.pi*3.1e6]     # stopband edge frequency (rad/s)
rp = 3                                  # maximum attenuation within pass band in dB
rs = 40                                 # minimum attenuation within stop band in dB

N, wn = scipy.signal.cheb2ord(wp, ws, rp, rs, analog=True)
z,p,k = scipy.signal.cheby2(N, rs, wn, 'bandpass', analog=True, output='zpk')

Similarly, you can also design Chebyshev-I type filters using cheb1ord() and cheby1(), Butterworth-type filters using buttord() and butter(), and elliptic-type filters using ellipord() and ellip() functions. Wikipedia describes the differences between these filter types in details (https://en.wikipedia.org/wiki/Chebyshev_filter).

In the above examples, the z, p, k results store the zero frequencies, pole frequencies, and gain value of the designed linear filter, respectively. There are a few considerations when you use them as the 'poles', 'zeros' and 'gain' parameters for the XMODEL's 'filter' primitive:

1. The 'poles' and 'zeros' parameters list the pole/zero frequency values in Hertz whereas SciPy returns pole/zeros frequencies in rad/s.
2. The real part of the pole/zero frequencies listed in the 'poles' and 'zeros' parameters is positive for left-half plane (LHP) poles/zeros and negative for right-half plane (RHP) poles/zeros. For SciPy results, it is negative for LHP and positive for RHP poles/zeros.
3. The 'poles' and 'zeros' parameters list the real and imaginary parts of each pole/zero frequency as separate numbers as below. It is because SystemVerilog does not natively support complex numbers:

poles = '{real(fp1), imag(fp1), real(fp2), imag(fp2), ..., real(fpM), imag(fpM)}
zeros = '{real(fz1), imag(fz1), real(fz2), imag(fz2), ..., real(fzN), imag(fzN)}

4. The complex poles/zeros listed in the 'poles' and 'zeros' parameters must be in exact conjugate pairs.
5. To express a filter transfer function without any zeros, the 'zeros' parameter should be defined as '{0.0}.
6. For a strictly proper system, the number of poles of a 'filter' primitive must be greater than the number of zeros. To satisfy this requirement, you may need to add a dummy pole (e.g. a pole at 1THz) that does not affect the filter characteristics.

For example, for z, p, k values shown below:

z = [ -3.72529030e-09+16418754.44763247j  -5.23868948e-09-16418754.44763249j
        1.92618690e-09 +6800870.76948295j   7.78721962e-11 -6800870.76948295j]
p = [ -1075430.76302337+2991438.64240786j -3170099.69183802+1512931.2492777j
      -3170099.69183801-1512931.2492777j  -1075430.76302337-2991438.64240787j]
k = (0.01+0j)

The corresponding filter instance may be:

filter #(.gain(0.01),
    .poles('{171160, 476102, 504537, 240790, 504537, -240790, 171160, -476102, 1e12, 0.0}),
    .zeros('{0.0, 2613125.9, 0.0, -2613125.9, 0.0, 1082392.2, 0.0, -1082392.2})) 
    filter_inst (.in(in), .out(out));

Note that the pole/zero frequencies are reordered to arrange the complex poles/zeros in conjugate pairs. Also, a dummy pole at a frequency well above the other pole/zero frequencies (1THz) is added in order to keep the filter strictly proper.

Here is one way to convert z, p, k values to gain, poles, zeros parameters:

    def sort_pz(p):
        res = []
        for w in sorted(p, key=np.real):
            f = w / (2*np.pi)
            if f.imag == 0.0:
                res += [f.real, 0.0]
            elif f.imag > 0.0:
                res += [f.real, f.imag, f.real, -f.imag]
        return res

    gain = (k / np.prod(-p) * np.prod(-z)).real
    poles = sort_pz(-p)
    zeros = sort_pz(-z)
    if len(p) <= len(z): poles += [1e12, 0.0] * (len(z)-len(p)+1)
    if len(z) == 0: zeros += [0.0]

    print("gain  : %g" % gain)
    print("poles : %s" % ', '.join(['%g' % v for v in poles]))
    print("zeros : %s" % ', '.join(['%g' % v for v in zeros]))