# Author: Jeffrey Armstrong # April 4, 2011 from __future__ import division, print_function, absolute_import import warnings import numpy as np from numpy.testing import (TestCase, run_module_suite, assert_equal, assert_array_almost_equal, assert_array_equal, assert_allclose, assert_, assert_raises, assert_almost_equal) from scipy.signal import (dlsim, dstep, dimpulse, tf2zpk, lti, dlti, StateSpace, TransferFunction, ZerosPolesGain, dfreqresp, dbode) class TestDLTI(TestCase): def test_dlsim(self): a = np.asarray([[0.9, 0.1], [-0.2, 0.9]]) b = np.asarray([[0.4, 0.1, -0.1], [0.0, 0.05, 0.0]]) c = np.asarray([[0.1, 0.3]]) d = np.asarray([[0.0, -0.1, 0.0]]) dt = 0.5 # Create an input matrix with inputs down the columns (3 cols) and its # respective time input vector u = np.hstack((np.asmatrix(np.linspace(0, 4.0, num=5)).transpose(), 0.01 * np.ones((5, 1)), -0.002 * np.ones((5, 1)))) t_in = np.linspace(0, 2.0, num=5) # Define the known result yout_truth = np.asmatrix([-0.001, -0.00073, 0.039446, 0.0915387, 0.13195948]).transpose() xout_truth = np.asarray([[0, 0], [0.0012, 0.0005], [0.40233, 0.00071], [1.163368, -0.079327], [2.2402985, -0.3035679]]) tout, yout, xout = dlsim((a, b, c, d, dt), u, t_in) assert_array_almost_equal(yout_truth, yout) assert_array_almost_equal(xout_truth, xout) assert_array_almost_equal(t_in, tout) # Make sure input with single-dimension doesn't raise error dlsim((1, 2, 3), 4) # Interpolated control - inputs should have different time steps # than the discrete model uses internally u_sparse = u[[0, 4], :] t_sparse = np.asarray([0.0, 2.0]) tout, yout, xout = dlsim((a, b, c, d, dt), u_sparse, t_sparse) assert_array_almost_equal(yout_truth, yout) assert_array_almost_equal(xout_truth, xout) assert_equal(len(tout), yout.shape[0]) # Transfer functions (assume dt = 0.5) num = np.asarray([1.0, -0.1]) den = np.asarray([0.3, 1.0, 0.2]) yout_truth = np.asmatrix([0.0, 0.0, 3.33333333333333, -4.77777777777778, 23.0370370370370]).transpose() # Assume use of the first column of the control input built earlier tout, yout = dlsim((num, den, 0.5), u[:, 0], t_in) assert_array_almost_equal(yout, yout_truth) assert_array_almost_equal(t_in, tout) # Retest the same with a 1-D input vector uflat = np.asarray(u[:, 0]) uflat = uflat.reshape((5,)) tout, yout = dlsim((num, den, 0.5), uflat, t_in) assert_array_almost_equal(yout, yout_truth) assert_array_almost_equal(t_in, tout) # zeros-poles-gain representation zd = np.array([0.5, -0.5]) pd = np.array([1.j / np.sqrt(2), -1.j / np.sqrt(2)]) k = 1.0 yout_truth = np.asmatrix([0.0, 1.0, 2.0, 2.25, 2.5]).transpose() tout, yout = dlsim((zd, pd, k, 0.5), u[:, 0], t_in) assert_array_almost_equal(yout, yout_truth) assert_array_almost_equal(t_in, tout) # Raise an error for continuous-time systems system = lti([1], [1, 1]) assert_raises(AttributeError, dlsim, system, u) def test_dstep(self): a = np.asarray([[0.9, 0.1], [-0.2, 0.9]]) b = np.asarray([[0.4, 0.1, -0.1], [0.0, 0.05, 0.0]]) c = np.asarray([[0.1, 0.3]]) d = np.asarray([[0.0, -0.1, 0.0]]) dt = 0.5 # Because b.shape[1] == 3, dstep should result in a tuple of three # result vectors yout_step_truth = (np.asarray([0.0, 0.04, 0.052, 0.0404, 0.00956, -0.036324, -0.093318, -0.15782348, -0.226628324, -0.2969374948]), np.asarray([-0.1, -0.075, -0.058, -0.04815, -0.04453, -0.0461895, -0.0521812, -0.061588875, -0.073549579, -0.08727047595]), np.asarray([0.0, -0.01, -0.013, -0.0101, -0.00239, 0.009081, 0.0233295, 0.03945587, 0.056657081, 0.0742343737])) tout, yout = dstep((a, b, c, d, dt), n=10) assert_equal(len(yout), 3) for i in range(0, len(yout)): assert_equal(yout[i].shape[0], 10) assert_array_almost_equal(yout[i].flatten(), yout_step_truth[i]) # Check that the other two inputs (tf, zpk) will work as well tfin = ([1.0], [1.0, 1.0], 0.5) yout_tfstep = np.asarray([0.0, 1.0, 0.0]) tout, yout = dstep(tfin, n=3) assert_equal(len(yout), 1) assert_array_almost_equal(yout[0].flatten(), yout_tfstep) zpkin = tf2zpk(tfin[0], tfin[1]) + (0.5,) tout, yout = dstep(zpkin, n=3) assert_equal(len(yout), 1) assert_array_almost_equal(yout[0].flatten(), yout_tfstep) # Raise an error for continuous-time systems system = lti([1], [1, 1]) assert_raises(AttributeError, dstep, system) def test_dimpulse(self): a = np.asarray([[0.9, 0.1], [-0.2, 0.9]]) b = np.asarray([[0.4, 0.1, -0.1], [0.0, 0.05, 0.0]]) c = np.asarray([[0.1, 0.3]]) d = np.asarray([[0.0, -0.1, 0.0]]) dt = 0.5 # Because b.shape[1] == 3, dimpulse should result in a tuple of three # result vectors yout_imp_truth = (np.asarray([0.0, 0.04, 0.012, -0.0116, -0.03084, -0.045884, -0.056994, -0.06450548, -0.068804844, -0.0703091708]), np.asarray([-0.1, 0.025, 0.017, 0.00985, 0.00362, -0.0016595, -0.0059917, -0.009407675, -0.011960704, -0.01372089695]), np.asarray([0.0, -0.01, -0.003, 0.0029, 0.00771, 0.011471, 0.0142485, 0.01612637, 0.017201211, 0.0175772927])) tout, yout = dimpulse((a, b, c, d, dt), n=10) assert_equal(len(yout), 3) for i in range(0, len(yout)): assert_equal(yout[i].shape[0], 10) assert_array_almost_equal(yout[i].flatten(), yout_imp_truth[i]) # Check that the other two inputs (tf, zpk) will work as well tfin = ([1.0], [1.0, 1.0], 0.5) yout_tfimpulse = np.asarray([0.0, 1.0, -1.0]) tout, yout = dimpulse(tfin, n=3) assert_equal(len(yout), 1) assert_array_almost_equal(yout[0].flatten(), yout_tfimpulse) zpkin = tf2zpk(tfin[0], tfin[1]) + (0.5,) tout, yout = dimpulse(zpkin, n=3) assert_equal(len(yout), 1) assert_array_almost_equal(yout[0].flatten(), yout_tfimpulse) # Raise an error for continuous-time systems system = lti([1], [1, 1]) assert_raises(AttributeError, dimpulse, system) def test_dlsim_trivial(self): a = np.array([[0.0]]) b = np.array([[0.0]]) c = np.array([[0.0]]) d = np.array([[0.0]]) n = 5 u = np.zeros(n).reshape(-1, 1) tout, yout, xout = dlsim((a, b, c, d, 1), u) assert_array_equal(tout, np.arange(float(n))) assert_array_equal(yout, np.zeros((n, 1))) assert_array_equal(xout, np.zeros((n, 1))) def test_dlsim_simple1d(self): a = np.array([[0.5]]) b = np.array([[0.0]]) c = np.array([[1.0]]) d = np.array([[0.0]]) n = 5 u = np.zeros(n).reshape(-1, 1) tout, yout, xout = dlsim((a, b, c, d, 1), u, x0=1) assert_array_equal(tout, np.arange(float(n))) expected = (0.5 ** np.arange(float(n))).reshape(-1, 1) assert_array_equal(yout, expected) assert_array_equal(xout, expected) def test_dlsim_simple2d(self): lambda1 = 0.5 lambda2 = 0.25 a = np.array([[lambda1, 0.0], [0.0, lambda2]]) b = np.array([[0.0], [0.0]]) c = np.array([[1.0, 0.0], [0.0, 1.0]]) d = np.array([[0.0], [0.0]]) n = 5 u = np.zeros(n).reshape(-1, 1) tout, yout, xout = dlsim((a, b, c, d, 1), u, x0=1) assert_array_equal(tout, np.arange(float(n))) # The analytical solution: expected = (np.array([lambda1, lambda2]) ** np.arange(float(n)).reshape(-1, 1)) assert_array_equal(yout, expected) assert_array_equal(xout, expected) def test_more_step_and_impulse(self): lambda1 = 0.5 lambda2 = 0.75 a = np.array([[lambda1, 0.0], [0.0, lambda2]]) b = np.array([[1.0, 0.0], [0.0, 1.0]]) c = np.array([[1.0, 1.0]]) d = np.array([[0.0, 0.0]]) n = 10 # Check a step response. ts, ys = dstep((a, b, c, d, 1), n=n) # Create the exact step response. stp0 = (1.0 / (1 - lambda1)) * (1.0 - lambda1 ** np.arange(n)) stp1 = (1.0 / (1 - lambda2)) * (1.0 - lambda2 ** np.arange(n)) assert_allclose(ys[0][:, 0], stp0) assert_allclose(ys[1][:, 0], stp1) # Check an impulse response with an initial condition. x0 = np.array([1.0, 1.0]) ti, yi = dimpulse((a, b, c, d, 1), n=n, x0=x0) # Create the exact impulse response. imp = (np.array([lambda1, lambda2]) ** np.arange(-1, n + 1).reshape(-1, 1)) imp[0, :] = 0.0 # Analytical solution to impulse response y0 = imp[:n, 0] + np.dot(imp[1:n + 1, :], x0) y1 = imp[:n, 1] + np.dot(imp[1:n + 1, :], x0) assert_allclose(yi[0][:, 0], y0) assert_allclose(yi[1][:, 0], y1) # Check that dt=0.1, n=3 gives 3 time values. system = ([1.0], [1.0, -0.5], 0.1) t, (y,) = dstep(system, n=3) assert_allclose(t, [0, 0.1, 0.2]) assert_array_equal(y.T, [[0, 1.0, 1.5]]) t, (y,) = dimpulse(system, n=3) assert_allclose(t, [0, 0.1, 0.2]) assert_array_equal(y.T, [[0, 1, 0.5]]) class TestDlti(object): def test_dlti_instantiation(self): # Test that lti can be instantiated. dt = 0.05 # TransferFunction s = dlti([1], [-1], dt=dt) assert_(isinstance(s, TransferFunction)) assert_(isinstance(s, dlti)) assert_(not isinstance(s, lti)) assert_equal(s.dt, dt) # ZerosPolesGain s = dlti(np.array([]), np.array([-1]), 1, dt=dt) assert_(isinstance(s, ZerosPolesGain)) assert_(isinstance(s, dlti)) assert_(not isinstance(s, lti)) assert_equal(s.dt, dt) # StateSpace s = dlti([1], [-1], 1, 3, dt=dt) assert_(isinstance(s, StateSpace)) assert_(isinstance(s, dlti)) assert_(not isinstance(s, lti)) assert_equal(s.dt, dt) # Number of inputs assert_raises(ValueError, dlti, 1) assert_raises(ValueError, dlti, 1, 1, 1, 1, 1) class TestStateSpaceDisc(object): def test_initialization(self): # Check that all initializations work dt = 0.05 s = StateSpace(1, 1, 1, 1, dt=dt) s = StateSpace([1], [2], [3], [4], dt=dt) s = StateSpace(np.array([[1, 2], [3, 4]]), np.array([[1], [2]]), np.array([[1, 0]]), np.array([[0]]), dt=dt) s = StateSpace(1, 1, 1, 1, dt=True) def _compare_systems(self, sys1, sys2): # Compare the contents of two systems assert_equal(sys1.A, sys2.A) assert_equal(sys1.B, sys2.B) assert_equal(sys1.C, sys2.C) assert_equal(sys1.D, sys2.D) def test_conversion(self): # Check the conversion functions s = StateSpace(1, 2, 3, 4, dt=0.05) assert_(isinstance(s.to_ss(), StateSpace)) assert_(isinstance(s.to_tf(), TransferFunction)) assert_(isinstance(s.to_zpk(), ZerosPolesGain)) # Make sure copies work assert_(StateSpace(s) is not s) assert_(s.to_ss() is not s) def test_properties(self): # Test setters/getters for cross class properties. # This implicitly tests to_tf() and to_zpk() # Getters s = StateSpace(1, 1, 1, 1, dt=0.05) assert_equal(s.poles, [1]) assert_equal(s.zeros, [0]) with warnings.catch_warnings(): warnings.simplefilter("ignore", DeprecationWarning) assert_equal(s.gain, 1) assert_equal(s.num, [1, 0]) assert_equal(s.den, [1, -1]) # transfer function setters s2 = StateSpace(2, 2, 2, 2, dt=0.05) with warnings.catch_warnings(): warnings.simplefilter("ignore", DeprecationWarning) s2.num = [1, 0] s2.den = [1, -1] self._compare_systems(s, s2) # zpk setters s2 = StateSpace(2, 2, 2, 2, dt=0.05) with warnings.catch_warnings(): warnings.simplefilter("ignore", DeprecationWarning) s2.poles = 1 s2.zeros = 0 s2.gain = 1 self._compare_systems(s, s2) class TestTransferFunction(object): def test_initialization(self): # Check that all initializations work dt = 0.05 s = TransferFunction(1, 1, dt=dt) s = TransferFunction([1], [2], dt=dt) s = TransferFunction(np.array([1]), np.array([2]), dt=dt) s = TransferFunction(1, 1, dt=True) def _compare_systems(self, sys1, sys2): # Compare the contents of two systems assert_equal(sys1.num, sys2.num) assert_equal(sys1.den, sys2.den) def test_conversion(self): # Check the conversion functions s = TransferFunction([1, 0], [1, -1], dt=0.05) assert_(isinstance(s.to_ss(), StateSpace)) assert_(isinstance(s.to_tf(), TransferFunction)) assert_(isinstance(s.to_zpk(), ZerosPolesGain)) # Make sure copies work assert_(TransferFunction(s) is not s) assert_(s.to_tf() is not s) def test_properties(self): # Test setters/getters for cross class properties. # This implicitly tests to_ss() and to_zpk() # Getters s = TransferFunction([1, 0], [1, -1], dt=0.05) assert_equal(s.poles, [1]) assert_equal(s.zeros, [0]) with warnings.catch_warnings(): warnings.simplefilter("ignore", DeprecationWarning) assert_equal(s.gain, 1) assert_equal(s.A, 1) assert_equal(s.B, 1) assert_equal(s.C, 1) assert_equal(s.D, 1) # state space setters s2 = TransferFunction([2, 3], [4, 5], dt=0.05) with warnings.catch_warnings(): warnings.simplefilter("ignore", DeprecationWarning) s2.A = 1 s2.B = 1 s2.C = 1 s2.D = 1 self._compare_systems(s, s2) # zpk setters s2 = TransferFunction([2, 3], [4, 5], dt=0.05) with warnings.catch_warnings(): warnings.simplefilter("ignore", DeprecationWarning) s2.poles = 1 s2.zeros = 0 s2.gain = 1 self._compare_systems(s, s2) class TestZerosPolesGain(object): def test_initialization(self): # Check that all initializations work dt = 0.05 s = ZerosPolesGain(1, 1, 1, dt=dt) s = ZerosPolesGain([1], [2], 1, dt=dt) s = ZerosPolesGain(np.array([1]), np.array([2]), 1, dt=dt) s = ZerosPolesGain(1, 1, 1, dt=True) def _compare_systems(self, sys1, sys2): # Compare the contents of two systems assert_equal(sys1.poles, sys2.poles) assert_equal(sys1.zeros, sys2.zeros) assert_equal(sys1.gain, sys2.gain) def test_conversion(self): # Check the conversion functions s = ZerosPolesGain(1, 2, 3, dt=0.05) assert_(isinstance(s.to_ss(), StateSpace)) assert_(isinstance(s.to_tf(), TransferFunction)) assert_(isinstance(s.to_zpk(), ZerosPolesGain)) # Make sure copies work assert_(ZerosPolesGain(s) is not s) assert_(s.to_zpk() is not s) def test_properties(self): # Test setters/getters for cross class properties. # This implicitly tests to_ss() and to_tf() # Getters s = ZerosPolesGain(0, 1, 1, dt=0.05) with warnings.catch_warnings(): warnings.simplefilter("ignore", DeprecationWarning) assert_equal(s.num, [1, 0]) assert_equal(s.den, [1, -1]) assert_equal(s.A, 1) assert_equal(s.B, 1) assert_equal(s.C, 1) assert_equal(s.D, 1) # state space setters s2 = ZerosPolesGain([2], [6], 3, dt=0.05) with warnings.catch_warnings(): warnings.simplefilter("ignore", DeprecationWarning) s2.A = 1 s2.B = 1 s2.C = 1 s2.D = 1 self._compare_systems(s, s2) # tf setters s2 = ZerosPolesGain([2], [5], 3, dt=0.05) with warnings.catch_warnings(): warnings.simplefilter("ignore", DeprecationWarning) s2.num = [1, 0] s2.den = [1, -1] self._compare_systems(s, s2) class Test_dfreqresp(TestCase): def test_manual(self): # Test dfreqresp() real part calculation (manual sanity check). # 1st order low-pass filter: H(z) = 1 / (z - 0.2), system = TransferFunction(1, [1, -0.2], dt=0.1) w = [0.1, 1, 10] w, H = dfreqresp(system, w=w) # test real expected_re = [1.2383, 0.4130, -0.7553] assert_almost_equal(H.real, expected_re, decimal=4) # test imag expected_im = [-0.1555, -1.0214, 0.3955] assert_almost_equal(H.imag, expected_im, decimal=4) def test_auto(self): # Test dfreqresp() real part calculation. # 1st order low-pass filter: H(z) = 1 / (z - 0.2), system = TransferFunction(1, [1, -0.2], dt=0.1) w = [0.1, 1, 10, 100] w, H = dfreqresp(system, w=w) jw = np.exp(w * 1j) y = np.polyval(system.num, jw) / np.polyval(system.den, jw) # test real expected_re = y.real assert_almost_equal(H.real, expected_re) # test imag expected_im = y.imag assert_almost_equal(H.imag, expected_im) def test_freq_range(self): # Test that freqresp() finds a reasonable frequency range. # 1st order low-pass filter: H(z) = 1 / (z - 0.2), # Expected range is from 0.01 to 10. system = TransferFunction(1, [1, -0.2], dt=0.1) n = 10 expected_w = np.linspace(0, np.pi, 10, endpoint=False) w, H = dfreqresp(system, n=n) assert_almost_equal(w, expected_w) def test_pole_one(self): # Test that freqresp() doesn't fail on a system with a pole at 0. # integrator, pole at zero: H(s) = 1 / s system = TransferFunction([1], [1, -1], dt=0.1) with warnings.catch_warnings(): warnings.simplefilter("ignore", RuntimeWarning) w, H = dfreqresp(system, n=2) assert_equal(w[0], 0.) # a fail would give not-a-number def test_error(self): # Raise an error for continuous-time systems system = lti([1], [1, 1]) assert_raises(AttributeError, dfreqresp, system) class Test_bode(object): def test_manual(self): # Test bode() magnitude calculation (manual sanity check). # 1st order low-pass filter: H(s) = 0.3 / (z - 0.2), dt = 0.1 system = TransferFunction(0.3, [1, -0.2], dt=dt) w = [0.1, 0.5, 1, np.pi] w2, mag, phase = dbode(system, w=w) # Test mag expected_mag = [-8.5329, -8.8396, -9.6162, -12.0412] assert_almost_equal(mag, expected_mag, decimal=4) # Test phase expected_phase = [-7.1575, -35.2814, -67.9809, -180.0000] assert_almost_equal(phase, expected_phase, decimal=4) # Test frequency assert_equal(np.array(w) / dt, w2) def test_auto(self): # Test bode() magnitude calculation. # 1st order low-pass filter: H(s) = 0.3 / (z - 0.2), system = TransferFunction(0.3, [1, -0.2], dt=0.1) w = np.array([0.1, 0.5, 1, np.pi]) w2, mag, phase = dbode(system, w=w) jw = np.exp(w * 1j) y = np.polyval(system.num, jw) / np.polyval(system.den, jw) # Test mag expected_mag = 20.0 * np.log10(abs(y)) assert_almost_equal(mag, expected_mag) # Test phase expected_phase = np.rad2deg(np.angle(y)) assert_almost_equal(phase, expected_phase) def test_range(self): # Test that bode() finds a reasonable frequency range. # 1st order low-pass filter: H(s) = 0.3 / (z - 0.2), dt = 0.1 system = TransferFunction(0.3, [1, -0.2], dt=0.1) n = 10 # Expected range is from 0.01 to 10. expected_w = np.linspace(0, np.pi, n, endpoint=False) / dt w, mag, phase = dbode(system, n=n) assert_almost_equal(w, expected_w) def test_pole_one(self): # Test that freqresp() doesn't fail on a system with a pole at 0. # integrator, pole at zero: H(s) = 1 / s system = TransferFunction([1], [1, -1], dt=0.1) with warnings.catch_warnings(): warnings.simplefilter("ignore", RuntimeWarning) w, mag, phase = dbode(system, n=2) assert_equal(w[0], 0.) # a fail would give not-a-number def test_imaginary(self): # bode() should not fail on a system with pure imaginary poles. # The test passes if bode doesn't raise an exception. system = TransferFunction([1], [1, 0, 100], dt=0.1) dbode(system, n=2) def test_error(self): # Raise an error for continuous-time systems system = lti([1], [1, 1]) assert_raises(AttributeError, dbode, system) class TestTransferFunctionZConversion(object): """Test private conversions between 'z' and 'z**-1' polynomials.""" def test_full(self): # Numerator and denominator same order num = [2, 3, 4] den = [5, 6, 7] num2, den2 = TransferFunction._z_to_zinv(num, den) assert_equal(num, num2) assert_equal(den, den2) num2, den2 = TransferFunction._zinv_to_z(num, den) assert_equal(num, num2) assert_equal(den, den2) def test_numerator(self): # Numerator lower order than denominator num = [2, 3] den = [5, 6, 7] num2, den2 = TransferFunction._z_to_zinv(num, den) assert_equal([0, 2, 3], num2) assert_equal(den, den2) num2, den2 = TransferFunction._zinv_to_z(num, den) assert_equal([2, 3, 0], num2) assert_equal(den, den2) def test_denominator(self): # Numerator higher order than denominator num = [2, 3, 4] den = [5, 6] num2, den2 = TransferFunction._z_to_zinv(num, den) assert_equal(num, num2) assert_equal([0, 5, 6], den2) num2, den2 = TransferFunction._zinv_to_z(num, den) assert_equal(num, num2) assert_equal([5, 6, 0], den2) if __name__ == "__main__": run_module_suite()