test_long_ZPK_stability

wield.control.SISO.test.test_long_ZPK_stability

This is a pytest module needing documentation

pytest-html report

Functions

chebcompanion2(c[, scale])

Return the scaled companion matrix of c.

gen_filt()

This is a pytest needing documentation

test_long_ZPK_stability_plot()

This is a test of the numerical stability of a very large ZPK filter that spans 18 orders of magnitude

test_long_cascade()

This is a demonstration showing the construction of a statespace using a Chechen companion matrix.

test_long_cheby()

This is a demonstration showing the construction of a statespace using a Chebychev companion matrix.

Details

chebcompanion2(c, scale=None)[source][github]

Return the scaled companion matrix of c.

The basis polynomials are scaled so that the companion matrix is symmetric when c is a Chebyshev basis polynomial. This provides better eigenvalue estimates than the unscaled case and for basis polynomials the eigenvalues are guaranteed to be real if numpy.linalg.eigvalsh is used to obtain them.

Parameters:

c (array_like) – 1-D array of Chebyshev series coefficients ordered from low to high degree.

Returns:

mat – Scaled companion matrix of dimensions (deg, deg).

Return type:

ndarray

Notes

Added in version 1.7.0.

code
docstring
"""Return the scaled companion matrix of c.

The basis polynomials are scaled so that the companion matrix is
symmetric when `c` is a Chebyshev basis polynomial. This provides
better eigenvalue estimates than the unscaled case and for basis
polynomials the eigenvalues are guaranteed to be real if
`numpy.linalg.eigvalsh` is used to obtain them.

Parameters
----------
c : array_like
    1-D array of Chebyshev series coefficients ordered from low to high
    degree.

Returns
-------
mat : ndarray
    Scaled companion matrix of dimensions (deg, deg).

Notes
-----

.. versionadded:: 1.7.0

"""

 1def chebcompanion2(c, scale = None):
 2
 3    # c is a trimmed copy
 4    if scale is None:
 5        scale = c[-1]
 6    if len(c) < 2:
 7        raise ValueError('Series must have maximum degree of at least 1.')
 8    if len(c) == 2:
 9        return np.array([[-c[0]/c[1]]])
10
11    n = len(c) - 1
12    mat = np.zeros((n, n), dtype=c.dtype)
13    scl = np.array([1.] + [np.sqrt(.5)]*(n-1))
14    top = mat.reshape(-1)[1::n+1]
15    bot = mat.reshape(-1)[n::n+1]
16    top[0] = np.sqrt(.5)
17    top[1:] = 1/2
18    bot[...] = top
19    mat[:, -1] -= (c[:-1]/scale)*(scl/scl[-1])*.5
20    return mat
pytest information

This code is wrapped in a pytest function using conventions detailed in Pytest Conventions. The full name of this test, as known by the documentation, is:

wield.control.SISO.test.test_long_ZPK_stability.chebcompanion2

The full name is useful when building documentation, to link a reference to this page using :func:`name`, or directly include it with an autofunction directive. The collapse nodes below show every instance of the test run. There may only be one, but if the test was run multiple times through pytest parametrizations, the list can be longer.

gen_filt()[source][github]

This is a pytest needing documentation

code
 1def gen_filt():
 2    F_Fq_lo = 10 * 2 * np.pi
 3    F_Fq_hi = 120 * 2 * np.pi
 4    hp_order = 1 #lo_ord
 5    lp_order = 2
 6    F_p = []
 7    F_z = []
 8    for ord in range(hp_order):
 9        F_p.append(-F_Fq_lo + 0*1j*F_Fq_lo)
10        F_p.append(-F_Fq_lo - 0*1j*F_Fq_lo)
11        F_z.append(0)
12        F_z.append(0)
13
14    for ord in range(lp_order):
15        F_p.append(-F_Fq_hi)
16    F_k = 1
17    c_zpk = cheby2(8, 160, 1.25*2*np.pi, btype='high', analog=True, output='zpk')
18    F_z.extend(c_zpk[0])
19    F_p.extend(c_zpk[1])
20    F_k *= c_zpk[2]
21
22    ian_lfq = -0.7
23    ian_lfq2 = -2
24    ian_lfq3 = -0.5
25    F_z.append(0)
26    F_p.append(ian_lfq)
27    F_z.append(0)
28    F_p.append(ian_lfq)
29
30    # F_z.append(ian_lfq3)
31    # F_z.append(ian_lfq3)
32    # F_p.append(ian_lfq2)
33    # F_p.append(ian_lfq2)
34
35    # make it more like a power spectrum
36    F_p = F_p
37    F_z = F_z
38    # F_z = [0] *  13
39    print()
40    print("Poles list: ", "[" + "; ".join(["{} + {}j".format(p.real, p.imag) for p in F_p]), "]")
41    print("Zeros list: ", "[" + "; ".join(["{} + {}j".format(z.real, z.imag) for z in F_z]), "]")
42
43    filt = SISO.zpk(F_z, F_p, F_k, angular=True, fiducial_rtol=1e-5, fiducial_atol=1e-10)
44    filtss = filt.asSS
45    filt = filt / abs(filtss.Linf_norm()[0]) 
46    return filt
pytest information

This code is wrapped in a pytest function using conventions detailed in Pytest Conventions. The full name of this test, as known by the documentation, is:

wield.control.SISO.test.test_long_ZPK_stability.gen_filt

The full name is useful when building documentation, to link a reference to this page using :func:`name`, or directly include it with an autofunction directive. The collapse nodes below show every instance of the test run. There may only be one, but if the test was run multiple times through pytest parametrizations, the list can be longer.

test_long_ZPK_stability_plot()[source][github]

This is a test of the numerical stability of a very large ZPK filter that spans 18 orders of magnitude

code
docstring
"""
This is a test of the numerical stability of a very large ZPK filter that spans 18 orders of magnitude

"""

 1def test_long_ZPK_stability_plot():
 2
 3    filt = gen_filt()
 4    filtss = filt.asSS
 5    filt = filt * filt
 6    filtss_bal2 = filtss.balance_and_truncate()
 7    filtss = filtss.balance()
 8    filtss = filtss.schur_form()
 9
10    filtss_bal3 = filtss.ss.transpose() @ filtss.ss.transpose()
11    filtss_bal3.print_nonzero()
12
13    filtss = filtss * filtss
14    filtss_bal2 = filtss_bal2 * filtss_bal2
15    # print(filtss.A.shape, filtss.E.shape)
16
17    filt = filt / abs(filtss.Linf_norm()[0]) 
18    filtss = filtss / abs(filtss.Linf_norm()[0]) 
19    # print(filtss.A.shape, filtss.E.shape)
20
21    filtss_bal = filtss.balance_and_truncate()
22    filtss_bal2 = filtss_bal2.balance_and_truncate()
23
24
25    filtss.print_nonzero()
26    filtss.print_nonzero()
27    filtss_bal.print_nonzero()
28
29    F_Hz = np.geomspace(1e-3, 1e3, 1000)
30
31    axB = mplfigB(Nrows=2)
32
33    mag, phase, omega = control.bode(control.ss(
34        filtss.A,
35        filtss.B,
36        filtss.C,
37        filtss.D,
38    ), Hz=False, plot=False, omega=F_Hz*2*np.pi, wrap_phase=True, label = "python-control")
39    axB.ax0.loglog(omega/(2*np.pi), mag, label="python-control")
40    axB.ax1.semilogx(omega/(2*np.pi), phase * 180 / np.pi, label="python-control")
41
42    xfer = filt.fresponse(f=F_Hz)
43    axB.ax0.loglog(*xfer.fplot_mag, label="ZPK, wield")
44    axB.ax1.semilogx(*xfer.fplot_deg180, label="ZPK, wield")
45
46    xfer = filtss.fresponse(f=F_Hz)
47    axB.ax0.loglog(*xfer.fplot_mag, label="state space, wield")
48    axB.ax1.semilogx(*xfer.fplot_deg180, label="state space, wield")
49
50    xfer = filtss_bal.fresponse(f=F_Hz)
51    axB.ax0.loglog(*xfer.fplot_mag, label="state space, wield bal")
52    axB.ax1.semilogx(*xfer.fplot_deg180, label="state space, wield bal")
53
54    xfer = filtss_bal2.fresponse(f=F_Hz)
55    axB.ax0.loglog(*xfer.fplot_mag, label="state space, wield bal2")
56    axB.ax1.semilogx(*xfer.fplot_deg180, label="state space, wield bal2")
57
58    axB.ax0.set_ylim(1e-25, 10)
59    axB.ax0.legend()
60
61    axB.save(tjoin('long_zpk'))
62    return
pytest information

This code is wrapped in a pytest function using conventions detailed in Pytest Conventions. The full name of this test, as known by the documentation, is:

wield.control.SISO.test.test_long_ZPK_stability.test_long_ZPK_stability_plot

The full name is useful when building documentation, to link a reference to this page using :func:`name`, or directly include it with an autofunction directive. The collapse nodes below show every instance of the test run. There may only be one, but if the test was run multiple times through pytest parametrizations, the list can be longer.

test_long_ZPK_stability_plot
output
Poles list:  [-62.83185307179586 + 0.0j; -62.83185307179586 + 0.0j; -753.9822368615503 + 0.0j; -753.9822368615503 + 0.0j; -8.284321652522745 + 42.35447200411846j; -23.59175208353443 + 35.90638758573623j; -35.30755211959318 + 23.991881149639656j; -41.64809740920085 + 8.42482829545465j; -41.64809740920085 + -8.42482829545465j; -35.30755211959318 + -23.991881149639656j; -23.59175208353443 + -35.90638758573623j; -8.284321652522745 + -42.35447200411846j; -0.7 + 0.0j; -0.7 + 0.0j ]
Zeros list:  [0 + 0j; 0 + 0j; -0.0 + 7.703069579159485j; -0.0 + 6.530347064232074j; -0.0 + 4.363438406518879j; -0.0 + 1.532235806080839j; 0.0 + -1.532235806080839j; 0.0 + -4.363438406518879j; 0.0 + -6.530347064232074j; 0.0 + -7.703069579159485j; 0 + 0j; 0 + 0j ]
   |   ABCDEFGHIJKLMNOPQRSTUVWXYZab   |       |                                 
0  | [[33333333............44545444]  | [[.]  | [[1...........................] 
1  |  [m3222222............44545444]  |  [.]  |  [.1..........................] 
2  |  [..233333............33445443]  |  [.]  |  [..1.........................] 
3  |  [..221111............33445443]  |  [.]  |  [...1........................] 
4  |  [....3333............33445443]  |  [.]  |  [....1.......................] 
5  |  [....1333............33445443]  |  [.]  |  [.....1......................] 
6  |  [......4m1111aa......44545444]  |  [.]  |  [......1.....................] 
7  |  [......441111aa......44545444]  |  [.]  |  [.......1....................] 
8  |  [........232211..............]  |  [.]  |  [........1...................] 
9  |  [........322211..............]  |  [.]  |  [.........1..................] 
10 |  [..........3321..............]  |  [.]  |  [..........1.................] 
11 |  [..........2311..............]  |  [.]  |  [...........1................] 
12 |  [............11..............]  |  [.]  |  [............1...............] 
13 |  [.............1..............]  |  [.]  |  [.............1..............] 
14 |  [..............33333333......]  |  [1]  |  [..............1.............] 
15 |  [..............m3222222......]  |  [a]  |  [...............1............] 
16 |  [................233333......]  |  [a]  |  [................1...........] 
17 |  [................221111......]  |  [a]  |  [.................1..........] 
18 |  [..................3333......]  |  [a]  |  [..................1.........] 
19 |  [..................1333......]  |  [a]  |  [...................1........] 
20 |  [....................4m1111aa]  |  [1]  |  [....................1.......] 
21 |  [....................441111aa]  |  [1]  |  [.....................1......] 
22 |  [......................232211]  |  [.]  |  [......................1.....] 
23 |  [......................322211]  |  [.]  |  [.......................1....] 
24 |  [........................3321]  |  [.]  |  [........................1...] 
25 |  [........................2311]  |  [.]  |  [.........................1..] 
26 |  [..........................11]  |  [.]  |  [..........................1.] 
27 |  [...........................1]] |  [.]] |  [...........................1]]
   |                                  |       |                                 
   | [[......44555554..............]] | [[.]] |                                 
   |   ABCDEFGHIJKLMNOPQRSTUVWXYZab   |       |                                 
0  | [[33333333....................]  | [[.]  | [[1...........................] 
1  |  [m3222222....................]  |  [.]  |  [.1..........................] 
2  |  [..233333....................]  |  [.]  |  [..1.........................] 
3  |  [..221111....................]  |  [.]  |  [...1........................] 
4  |  [....3333....................]  |  [.]  |  [....1.......................] 
5  |  [....1333....................]  |  [.]  |  [.....1......................] 
6  |  [......4m1111aa44333344......]  |  [.]  |  [......1.....................] 
7  |  [......441111aa44333344......]  |  [.]  |  [.......1....................] 
8  |  [........23221155444455......]  |  [.]  |  [........1...................] 
9  |  [........32221144444444......]  |  [.]  |  [.........1..................] 
10 |  [..........332155555555......]  |  [.]  |  [..........1.................] 
11 |  [..........231144444444......]  |  [.]  |  [...........1................] 
12 |  [............1144444444......]  |  [.]  |  [............1...............] 
13 |  [.............144333344......]  |  [.]  |  [.............1..............] 
14 |  [..............33333333......]  |  [.]  |  [..............1.............] 
15 |  [..............m3222222......]  |  [.]  |  [...............1............] 
16 |  [................233333......]  |  [.]  |  [................1...........] 
17 |  [................221111......]  |  [.]  |  [.................1..........] 
18 |  [..................3333......]  |  [.]  |  [..................1.........] 
19 |  [..................1333......]  |  [.]  |  [...................1........] 
20 |  [....................4m1111aa]  |  [4]  |  [....................1.......] 
21 |  [....................441111aa]  |  [4]  |  [.....................1......] 
22 |  [......................232211]  |  [5]  |  [......................1.....] 
23 |  [......................322211]  |  [5]  |  [.......................1....] 
24 |  [........................3321]  |  [5]  |  [........................1...] 
25 |  [........................2311]  |  [5]  |  [.........................1..] 
26 |  [..........................11]  |  [5]  |  [..........................1.] 
27 |  [...........................1]] |  [4]] |  [...........................1]]
   |                                  |       |                                 
   | [[1aaaaa11....................]] | [[.]] |                                 
   |   ABCDEFGHIJKLMNOPQRSTUVWXYZab   |       |                                 
0  | [[33333333....................]  | [[.]  | [[1...........................] 
1  |  [m3222222....................]  |  [.]  |  [.1..........................] 
2  |  [..233333....................]  |  [.]  |  [..1.........................] 
3  |  [..221111....................]  |  [.]  |  [...1........................] 
4  |  [....3333....................]  |  [.]  |  [....1.......................] 
5  |  [....1333....................]  |  [.]  |  [.....1......................] 
6  |  [......4m1111aa44333344......]  |  [.]  |  [......1.....................] 
7  |  [......441111aa44333344......]  |  [.]  |  [.......1....................] 
8  |  [........23221155444455......]  |  [.]  |  [........1...................] 
9  |  [........32221144444444......]  |  [.]  |  [.........1..................] 
10 |  [..........332155555555......]  |  [.]  |  [..........1.................] 
11 |  [..........231144444444......]  |  [.]  |  [...........1................] 
12 |  [............1144444444......]  |  [.]  |  [............1...............] 
13 |  [.............144333344......]  |  [.]  |  [.............1..............] 
14 |  [..............33333333......]  |  [.]  |  [..............1.............] 
15 |  [..............m3222222......]  |  [.]  |  [...............1............] 
16 |  [................233333......]  |  [.]  |  [................1...........] 
17 |  [................221111......]  |  [.]  |  [.................1..........] 
18 |  [..................3333......]  |  [.]  |  [..................1.........] 
19 |  [..................1333......]  |  [.]  |  [...................1........] 
20 |  [....................4m1111aa]  |  [4]  |  [....................1.......] 
21 |  [....................441111aa]  |  [4]  |  [.....................1......] 
22 |  [......................232211]  |  [5]  |  [......................1.....] 
23 |  [......................322211]  |  [5]  |  [.......................1....] 
24 |  [........................3321]  |  [5]  |  [........................1...] 
25 |  [........................2311]  |  [5]  |  [.........................1..] 
26 |  [..........................11]  |  [5]  |  [..........................1.] 
27 |  [...........................1]] |  [4]] |  [...........................1]]
   |                                  |       |                                 
   | [[1aaaaa11....................]] | [[.]] |                                 
   |   ABCDEFGHIJKLMNOPQRSTUVWX   |       |                             
0  | [[233122321321a21aaabbccde]  | [[4]  | [[1.......................] 
1  |  [312122311221a21aabbccdee]  |  [4]  |  [.1......................] 
2  |  [32332332233213211aabbcde]  |  [4]  |  [..1.....................] 
3  |  [113b111aa11ab1abbccddeef]  |  [3]  |  [...1....................] 
4  |  [222123312321a211aabbccde]  |  [4]  |  [....1...................] 
5  |  [22313333233213211aabbcdd]  |  [4]  |  [.....1..................] 
6  |  [333133332432132211aabbcd]  |  [4]  |  [......1.................] 
7  |  [212a13312221a21aaabbcdde]  |  [3]  |  [.......1................] 
8  |  [112a22222332a211aaabccde]  |  [3]  |  [........1...............] 
9  |  [3231334234432432211aabcd]  |  [4]  |  [.........1..............] 
10 |  [2231233234331422211aabcd]  |  [4]  |  [..........1.............] 
11 |  [112a12212332132211aabbcd]  |  [3]  |  [...........1............] 
12 |  [aa1ba11aa211a311aabbccde]  |  [2]  |  [............1...........] 
13 |  [223123322443343332211abc]  |  [4]  |  [.............1..........] 
14 |  [112a1221132213322211abbc]  |  [3]  |  [..............1.........] 
15 |  [aa1b112a122213222211abbc]  |  [2]  |  [...............1........] 
16 |  [aa1ba11aa221a32222211abc]  |  [2]  |  [................1.......] 
17 |  [abacaa1aa111a22222221aab]  |  [1]  |  [.................1......] 
18 |  [bbacbaaba11ab211222211ab]  |  [1]  |  [..................1.....] 
19 |  [bcbdbbabbaaab1111222211a]  |  [a]  |  [...................1....] 
20 |  [ccbdcbbccaabc1aa1112221a]  |  [a]  |  [....................1...] 
21 |  [cdceccbdcbbbcabbaa112221]  |  [b]  |  [.....................1..] 
22 |  [dededdcddcccdbbbbaa11222]  |  [c]  |  [......................1.] 
23 |  [eeefeddeedddeccccbbaa122]] |  [d]] |  [.......................1]]
   |                              |       |                             
   | [[bbacbaabcabcdbcdddeefggh]] | [[.]] |                             
captured errors:
def test_long_ZPK_stability_plot():
        """
        This is a test of the numerical stability of a very large ZPK filter that spans 18 orders of magnitude

        """
        filt = gen_filt()
        filtss = filt.asSS
        filt = filt * filt
        filtss_bal2 = filtss.balance_and_truncate()
        filtss = filtss.balance()
        filtss = filtss.schur_form()

        filtss_bal3 = filtss.ss.transpose() @ filtss.ss.transpose()
        filtss_bal3.print_nonzero()

        filtss = filtss * filtss
        filtss_bal2 = filtss_bal2 * filtss_bal2
        # print(filtss.A.shape, filtss.E.shape)

        filt = filt / abs(filtss.Linf_norm()[0])
        filtss = filtss / abs(filtss.Linf_norm()[0])
        # print(filtss.A.shape, filtss.E.shape)

        filtss_bal = filtss.balance_and_truncate()
        filtss_bal2 = filtss_bal2.balance_and_truncate()


        filtss.print_nonzero()
        filtss.print_nonzero()
        filtss_bal.print_nonzero()

        F_Hz = np.geomspace(1e-3, 1e3, 1000)

        axB = mplfigB(Nrows=2)

>       mag, phase, omega = control.bode(control.ss(
            filtss.A,
            filtss.B,
            filtss.C,
            filtss.D,
        ), Hz=False, plot=False, omega=F_Hz*2*np.pi, wrap_phase=True, label = "python-control")
E       AttributeError: module 'control' has no attribute 'bode'

../../src/wield/control/SISO/test/test_long_ZPK_stability.py:94: AttributeError
test_long_cascade()[source][github]

This is a demonstration showing the construction of a statespace using a Chechen companion matrix.

code
docstring
"""
This is a demonstration showing the construction of a statespace using a
Chechen companion matrix. 
"""

 1def test_long_cascade():
 2
 3    filt = SISO.zpk(
 4        [-3000] * 8,
 5        [-1] * 8,
 6        1,
 7        angular=True,
 8        fiducial_rtol=1e-5,
 9        fiducial_atol=1e-10
10    )
11    filt = gen_filt()
12    filtss = filt.asSS
13    norm = abs(filtss.Linf_norm()[0])
14    # filt = filt.inv()
15    # filtss = filtss.inv()
16
17    print(filtss.A)
18
19    axB = mplfigB(Nrows=2)
20
21    F_Hz = np.geomspace(1e-3, 1e5, 1000)
22
23    xfer = filt.fresponse(f=F_Hz)
24    axB.ax0.loglog(*xfer.fplot_mag, label="ZPK, wield")
25    axB.ax1.semilogx(*xfer.fplot_deg180, label="ZPK, wield")
26    axB.ax0.axhline(1e-16, ls='--', color='black', lw=1)
27
28    xfer = filtss.fresponse(f=F_Hz)
29    axB.ax0.loglog(*xfer.fplot_mag, label="filt2")
30    axB.ax1.semilogx(*xfer.fplot_deg180, label="filt2")
31    #axB.ax0.set_ylim(1e-25, 10)
32    axB.ax0.legend()
33
34    axB.save(tjoin('long_zpk'))
35    return
pytest information

This code is wrapped in a pytest function using conventions detailed in Pytest Conventions. The full name of this test, as known by the documentation, is:

wield.control.SISO.test.test_long_ZPK_stability.test_long_cascade

The full name is useful when building documentation, to link a reference to this page using :func:`name`, or directly include it with an autofunction directive. The collapse nodes below show every instance of the test run. There may only be one, but if the test was run multiple times through pytest parametrizations, the list can be longer.

test_long_cascade
output
Poles list:  [-62.83185307179586 + 0.0j; -62.83185307179586 + 0.0j; -753.9822368615503 + 0.0j; -753.9822368615503 + 0.0j; -8.284321652522745 + 42.35447200411846j; -23.59175208353443 + 35.90638758573623j; -35.30755211959318 + 23.991881149639656j; -41.64809740920085 + 8.42482829545465j; -41.64809740920085 + -8.42482829545465j; -35.30755211959318 + -23.991881149639656j; -23.59175208353443 + -35.90638758573623j; -8.284321652522745 + -42.35447200411846j; -0.7 + 0.0j; -0.7 + 0.0j ]
Zeros list:  [0 + 0j; 0 + 0j; -0.0 + 7.703069579159485j; -0.0 + 6.530347064232074j; -0.0 + 4.363438406518879j; -0.0 + 1.532235806080839j; 0.0 + -1.532235806080839j; 0.0 + -4.363438406518879j; 0.0 + -6.530347064232074j; 0.0 + -7.703069579159485j; 0 + 0j; 0 + 0j ]
[[-1.25663706e+02  9.05096680e+01 -6.28318531e+01  0.00000000e+00
  -6.28318531e+01  0.00000000e+00 -1.25663706e+02  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00]
 [-4.36179012e+01  0.00000000e+00 -2.14811554e+01  0.00000000e+00
  -2.14811554e+01  0.00000000e+00 -4.29623107e+01  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00 -4.71835042e+01  4.52548340e+01
  -4.71835042e+01  0.00000000e+00 -9.43670083e+01  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00 -4.07876744e+01  0.00000000e+00
  -4.03669549e+01  0.00000000e+00 -8.07339097e+01  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  -8.32961948e+01  4.52548340e+01 -1.66592390e+02  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  -3.98972130e+01  0.00000000e+00 -7.97944259e+01  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00 -1.50796447e+03  7.24077344e+02
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00 -7.85122222e+02  0.00000000e+00
   5.65685425e+00  0.00000000e+00  5.65685425e+00  0.00000000e+00
   7.07106781e-01  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  -1.65686433e+01  4.52548340e+01 -1.65686433e+01  0.00000000e+00
  -2.07108041e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  -4.11565157e+01  0.00000000e+00 -4.02141758e+01  0.00000000e+00
  -5.02677198e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00 -7.06151042e+01  4.52548340e+01
  -8.82688803e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00 -4.02660542e+01  0.00000000e+00
  -5.02677198e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  -1.40000000e+00  1.41421356e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  -3.46482323e-01  0.00000000e+00]]
test_long_cheby()[source][github]

This is a demonstration showing the construction of a statespace using a Chebychev companion matrix.

code
docstring
"""
This is a demonstration showing the construction of a statespace using a
Chebychev companion matrix. 
"""

 1def test_long_cheby():
 2
 3    filt = gen_filt()
 4    filt = filt * filt
 5    p = filt.p
 6    z = filt.z
 7    norm = max(abs(p))**0.5
 8    norm = 1
 9    p = p / norm
10    z = z / norm
11    from numpy.polynomial.chebyshev import (
12        chebcompanion,
13        chebfromroots,
14        chebroots,
15        chebdiv,
16    )
17    cp = chebfromroots(p).real
18    cz = chebfromroots(z).real
19    nd = len(cp) - len(cz)
20    if nd == 0:
21        qz, cz = chebdiv(cz, cp)
22        # surprisingly, the nd needs to stay 0 for the normalization call below,
23        # rather than being reset
24        # nd = len(cp) - len(cz)
25    else:
26        qz = 0
27    print("nd: ", nd)
28    czl = np.concatenate([cz, np.zeros(len(cp) - len(cz))])
29    print("cp: ", cp)
30    print("cz: ", cz)
31    cp_alt = chebfromroots(p * 1j)
32    print("cp alt!: ", cp_alt)
33    A = chebcompanion(cp)
34    Az = chebcompanion2(czl, scale = cp[-1])
35
36    n = len(cp) - 1
37    scl = np.array([1.] + [np.sqrt(.5)]*(n-1))
38    # this is basically the last line of the companion matrix, but not subtracting
39    # the companion coupling alpha beta gamma contribution.
40    czs = (czl[:-1]/cp[-1])*(scl/scl[-1])*.5
41
42    print(A[:, -1])
43    print(Az[:, -1])
44    print("poles1: ", p * norm)
45    print("poles1: ", chebroots(cp) * norm)
46    p2 = chebroots(cp) * norm
47
48    print("poles2: ", np.linalg.eigvals(A) * norm)
49    z2 = z * norm
50
51    filt2 = SISO.zpk(z2, p2, filt.k, angular=True, fiducial_rtol=1e-5, fiducial_atol=1e-10)
52
53    C = np.zeros(len(cp) - 1)
54    C[-1] = filt.k
55    B = czs
56    # B = np.zeros(len(cp) - 1)
57    # B[0] = 1 / np.sqrt(.5)*(n-1)
58    D = np.asarray([qz]).reshape(1, 1) * filt.k
59    filt3 = SISO.SISOStateSpace(
60        A=A * norm,
61        B=B.reshape(-1, 1) / norm**(nd - 1),
62        C=C.reshape(1, -1),
63        # C=B.reshape(1, -1),
64        # B=C.reshape(-1, 1),
65        D = D,
66    )
67
68    axB = mplfigB(Nrows=2)
69
70    F_Hz = np.geomspace(1e-3, 1e5, 1000)
71
72    #filt = filt * filt
73    xfer = filt.fresponse(f=F_Hz)
74    axB.ax0.loglog(*xfer.fplot_mag, label="ZPK, wield")
75    axB.ax1.semilogx(*xfer.fplot_deg180, label="ZPK, wield")
76    axB.ax0.axhline(1e-16, ls='--', color='black', lw=1)
77
78    xfer = filt.asSS.fresponse(f=F_Hz)
79    axB.ax0.loglog(*xfer.fplot_mag, label="SS, wield")
80    axB.ax1.semilogx(*xfer.fplot_deg180, label="SS, wield")
81    axB.ax0.axhline(1e-16, ls='--', color='black', lw=1)
82
83    #filt2 = filt2 * filt2
84    xfer = filt2.fresponse(f=F_Hz)
85    axB.ax0.loglog(*xfer.fplot_mag, label="filt2")
86    axB.ax1.semilogx(*xfer.fplot_deg180, label="filt2")
87    #axB.ax0.set_ylim(1e-25, 10)
88    axB.ax0.legend()
89
90    #filt3 = filt3 * filt3
91    xfer = filt3.fresponse(f=F_Hz)
92    axB.ax0.loglog(*xfer.fplot_mag, label="filt3")
93    axB.ax1.semilogx(*xfer.fplot_deg180, label="filt3")
94    axB.ax0.set_ylim(1e-25, 10)
95    axB.ax0.legend()
96
97    axB.save(tjoin('long_zpk'))
98    return
pytest information

This code is wrapped in a pytest function using conventions detailed in Pytest Conventions. The full name of this test, as known by the documentation, is:

wield.control.SISO.test.test_long_ZPK_stability.test_long_cheby

The full name is useful when building documentation, to link a reference to this page using :func:`name`, or directly include it with an autofunction directive. The collapse nodes below show every instance of the test run. There may only be one, but if the test was run multiple times through pytest parametrizations, the list can be longer.

test_long_cheby
output
Poles list:  [-62.83185307179586 + 0.0j; -62.83185307179586 + 0.0j; -753.9822368615503 + 0.0j; -753.9822368615503 + 0.0j; -8.284321652522745 + 42.35447200411846j; -23.59175208353443 + 35.90638758573623j; -35.30755211959318 + 23.991881149639656j; -41.64809740920085 + 8.42482829545465j; -41.64809740920085 + -8.42482829545465j; -35.30755211959318 + -23.991881149639656j; -23.59175208353443 + -35.90638758573623j; -8.284321652522745 + -42.35447200411846j; -0.7 + 0.0j; -0.7 + 0.0j ]
Zeros list:  [0 + 0j; 0 + 0j; -0.0 + 7.703069579159485j; -0.0 + 6.530347064232074j; -0.0 + 4.363438406518879j; -0.0 + 1.532235806080839j; 0.0 + -1.532235806080839j; 0.0 + -4.363438406518879j; 0.0 + -6.530347064232074j; 0.0 + -7.703069579159485j; 0 + 0j; 0 + 0j ]
nd:  4
cp:  [1.73871847e+45 2.93833543e+45 1.75173170e+45 6.99802619e+44
 1.68985598e+44 1.86505237e+43 1.24451973e+42 5.80798455e+40
 2.05363601e+39 5.77444604e+37 1.33201211e+36 2.57336770e+34
 4.22156066e+32 5.93264302e+30 7.17756842e+28 7.48841715e+26
 6.72819927e+24 5.18298064e+22 3.39651603e+20 1.87118441e+18
 8.52104331e+15 3.13290481e+13 9.00227663e+10 1.93287568e+08
 2.91727598e+05 2.86671187e+02 1.70478558e-01 5.52145131e-05
 7.45058060e-09]
cz:  [7.91503198e+09 0.00000000e+00 1.30181133e+10 0.00000000e+00
 7.15492826e+09 0.00000000e+00 2.52858686e+09 0.00000000e+00
 5.33062363e+08 0.00000000e+00 5.95856036e+07 0.00000000e+00
 3.35321350e+06 0.00000000e+00 9.13324308e+04 0.00000000e+00
 1.32555695e+03 0.00000000e+00 1.08198982e+01 0.00000000e+00
 4.97953370e-02 0.00000000e+00 1.20515876e-04 0.00000000e+00
 1.19209290e-07]
cp alt!:  [-4.59758929e+44-2.37684488e+29j -2.97105609e+29+7.09862211e+44j
 -4.59027297e+44-3.16912650e+29j -2.27780967e+29+5.17337747e+44j
  1.54247633e+44+5.94211219e+28j  5.26124517e+27-1.78446573e+43j
 -1.21187413e+42-2.32113757e+26j -1.20892582e+25+5.70455301e+40j
  2.02709692e+39+4.15568250e+23j  9.44473297e+21-5.71800141e+37j
 -1.32190446e+36-2.21360929e+20j -2.88230376e+18+2.55797201e+34j
  4.20149359e+32+5.40431955e+16j  7.74056186e+14-5.91020419e+30j
 -7.15605774e+28-8.35800636e+12j -7.15380490e+10+7.47080710e+26j
  6.71597795e+24+5.81435392e+08j  5.08723200e+06-5.17587264e+22j
 -3.39310841e+20-3.42400000e+04j -1.93500000e+02+1.86986878e+18j
  8.51708261e+15+8.12500000e-01j  2.68554688e-03-3.13201572e+13j
 -9.00087636e+10-5.66244125e-06j -8.49831849e-09+1.93273235e+08j
  2.91718733e+05+6.50857146e-12j  2.74780199e-15-2.86668206e+02j
 -1.70478141e-01-5.14996032e-19j  0.00000000e+00+5.52145131e-05j
  7.45058060e-09+0.00000000e+00j]
[-1.65015277e+53 -1.97188353e+53 -1.17556724e+53 -4.69629588e+52
 -1.13404315e+52 -1.25161546e+51 -8.35183053e+49 -3.89767245e+48
 -1.37817179e+47 -3.87516514e+45 -8.93898196e+43 -1.72695783e+42
 -2.83304140e+40 -3.98132934e+38 -4.81678463e+36 -5.02539168e+34
 -4.51521810e+32 -3.47823943e+30 -2.27936333e+28 -1.25573060e+26
 -5.71837536e+23 -2.10245683e+21 -6.04132558e+18 -1.29713091e+16
 -1.95775077e+13 -1.92381777e+10 -1.14406218e+07 -3.70538325e+03]
[-7.51186128e+17  0.00000000e+00 -8.73630797e+17  0.00000000e+00
 -4.80159107e+17  0.00000000e+00 -1.69690592e+17  0.00000000e+00
 -3.57732096e+16  0.00000000e+00 -3.99872217e+15  0.00000000e+00
 -2.25030349e+14  0.00000000e+00 -6.12921568e+12  0.00000000e+00
 -8.89566212e+10  0.00000000e+00 -7.26111075e+08  0.00000000e+00
 -3.34170850e+06  0.00000000e+00 -8.08768352e+03  0.00000000e+00
 -8.00000000e+00  0.00000000e+00  5.00000000e-01  0.00000000e+00]
poles1:  [-6.28318531e+01 +0.j         -6.28318531e+01 +0.j
 -7.53982237e+02 +0.j         -7.53982237e+02 +0.j
 -7.00000000e-01 +0.j         -7.00000000e-01 +0.j
 -6.28318531e+01 +0.j         -6.28318531e+01 +0.j
 -7.53982237e+02 +0.j         -7.53982237e+02 +0.j
 -7.00000000e-01 +0.j         -7.00000000e-01 +0.j
 -8.28432165e+00+42.354472j   -8.28432165e+00-42.354472j
 -2.35917521e+01+35.90638759j -2.35917521e+01-35.90638759j
 -3.53075521e+01+23.99188115j -3.53075521e+01-23.99188115j
 -4.16480974e+01 +8.4248283j  -4.16480974e+01 -8.4248283j
 -8.28432165e+00+42.354472j   -8.28432165e+00-42.354472j
 -2.35917521e+01+35.90638759j -2.35917521e+01-35.90638759j
 -3.53075521e+01+23.99188115j -3.53075521e+01-23.99188115j
 -4.16480974e+01 +8.4248283j  -4.16480974e+01 -8.4248283j ]
poles1:  [-7.54139995e+02-1.57768934e-01j -7.54139995e+02+1.57768934e-01j
 -7.53824479e+02-1.57746950e-01j -7.53824479e+02+1.57746950e-01j
 -6.31121508e+01-2.79649662e-01j -6.31121508e+01+2.79649662e-01j
 -6.25515382e+01-2.80512935e-01j -6.25515382e+01+2.80512935e-01j
 -4.16584727e+01-8.41704179e+00j -4.16584727e+01+8.41704179e+00j
 -4.16377405e+01-8.43258575e+00j -4.16377405e+01+8.43258575e+00j
 -3.53102836e+01-2.39944354e+01j -3.53102836e+01+2.39944354e+01j
 -3.53048195e+01-2.39893297e+01j -3.53048195e+01+2.39893297e+01j
 -2.35918233e+01-3.59058079e+01j -2.35918233e+01+3.59058079e+01j
 -2.35916808e+01-3.59069672e+01j -2.35916808e+01+3.59069672e+01j
 -8.28433625e+00-4.23544616e+01j -8.28433625e+00+4.23544616e+01j
 -8.28430705e+00-4.23544824e+01j -8.28430705e+00+4.23544824e+01j
 -7.00379295e-01+0.00000000e+00j -7.00000023e-01-3.79318221e-04j
 -7.00000023e-01+3.79318221e-04j -6.99620659e-01+0.00000000e+00j]
poles2:  [-8.28418199e+00+4.23544686e+01j -8.28418199e+00-4.23544686e+01j
 -8.28446132e+00+4.23544754e+01j -8.28446132e+00-4.23544754e+01j
 -6.99449440e-01+5.50482067e-04j -6.99449440e-01-5.50482067e-04j
 -7.00550560e-01+5.50639497e-04j -7.00550560e-01-5.50639497e-04j
 -2.35907098e+01+3.59069203e+01j -2.35907098e+01-3.59069203e+01j
 -2.35927946e+01+3.59058547e+01j -2.35927946e+01-3.59058547e+01j
 -3.53050864e+01+2.39969196e+01j -3.53050864e+01-2.39969196e+01j
 -3.53100174e+01+2.39868363e+01j -3.53100174e+01-2.39868363e+01j
 -4.16337614e+01+8.41256983e+00j -4.16337614e+01-8.41256983e+00j
 -4.16624361e+01+8.43715375e+00j -4.16624361e+01-8.43715375e+00j
 -6.21115754e+01+0.00000000e+00j -6.28418801e+01+7.09791433e-01j
 -6.28418801e+01-7.09791433e-01j -6.35320717e+01+0.00000000e+00j
 -7.53816469e+02+1.65585782e-01j -7.53816469e+02-1.65585782e-01j
 -7.54148005e+02+1.65948995e-01j -7.54148005e+02-1.65948995e-01j]