def simulation(self, ts_length=90, random_state=None):
        """
        Compute a simulated sample path assuming Gaussian shocks.

        Parameters
        ----------
        ts_length : scalar(int), optional(default=90)
            Number of periods to simulate for

        random_state : int or np.random.RandomState, optional
            Random seed (integer) or np.random.RandomState instance to set
            the initial state of the random number generator for
            reproducibility. If None, a randomly initialized RandomState is
            used.

        Returns
        -------
        vals : array_like(float)
            A simulation of the model that corresponds to this class

        """
        random_state = check_random_state(random_state)

        sys = self.ma_poly, self.ar_poly, 1
        u = random_state.randn(ts_length, 1) * self.sigma
        vals = dlsim(sys, u)[1]

        return vals.flatten()

 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 main():
    dirac = zeros((100, 1))
    dirac[0] = 1

    for i in [
        [0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9],
    ]:
        reflectance = check_surface_filters(i, False)
        plt.figure()
        plt.plot(*dlsim((reflectance[0], reflectance[1], 1), dirac))

        plt.show()

    def simulation(self, ts_length=90):
        """
        Compute a simulated sample path assuming Gaussian shocks.
        Parameters
        ----------
        ts_length : scalar(int), optional(default=90)
            Number of periods to simulate for
        Returns
        -------
        vals : array_like(float)
            A simulation of the model that corresponds to this class
        """
        sys = self.ma_poly, self.ar_poly, 1
        u = np.random.randn(ts_length, 1) * self.sigma
        vals = dlsim(sys, u)[1]

        return vals.flatten()

 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_discrete_approx(self):
        """
        Test that the solution to the discrete approximation of a continuous
        system actually approximates the solution to the continuous sytem.
        This is an indirect test of the correctness of the implementation
        of cont2discrete.
        """

        def u(t):
            return np.sin(2.5 * t)

        a = np.array([[-0.01]])
        b = np.array([[1.0]])
        c = np.array([[1.0]])
        d = np.array([[0.2]])
        x0 = 1.0

        t = np.linspace(0, 10.0, 101)
        dt = t[1] - t[0]
        u1 = u(t)

        # Use lsim2 to compute the solution to the continuous system.
        t, yout, xout = lsim2((a, b, c, d), T=t, U=u1, X0=x0,
                              rtol=1e-9, atol=1e-11)

        # Convert the continuous system to a discrete approximation.
        dsys = c2d((a, b, c, d), dt, method='bilinear')

        # Use dlsim with the pairwise averaged input to compute the output
        # of the discrete system.
        u2 = 0.5 * (u1[:-1] + u1[1:])
        t2 = t[:-1]
        td2, yd2, xd2 = dlsim(dsys, u=u2.reshape(-1, 1), t=t2, x0=x0)

        # ymid is the average of consecutive terms of the "exact" output
        # computed by lsim2.  This is what the discrete approximation
        # actually approximates.
        ymid = 0.5 * (yout[:-1] + yout[1:])

        assert_allclose(yd2.ravel(), ymid, rtol=1e-4)

def runLinearAero(E,F,G,C,D,delS,nT,u,x0 = None):
    """@details run time-domain simulation of linear aerodynamics.
    @param E Discrete-time state-space matrix.
    @param F Discrete-time state-space matrix.
    @param G Discrete-time state-space matrix.
    @param C Discrete-time state-space matrix.
    @param D Discrete-time state-space matrix.
    @param delS Non-dimensional time step of model.
    @param nT Number of time-steps to run simulation.
    @param u Inputs, an nT x m array.
    @return x State history, an nT x n array.
    @return y Output history, an nT x l array.
    """
    
    invE = np.linalg.inv(E)    
    # run simulation
    tOut, yOut, xOut = dlsim((np.dot(invE,F),np.dot(invE,G),C,D,delS),
                             u,
                             None,
                             x0)
    
    return tOut, yOut, xOut

    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 simulation(self, ts_length=90) :
     " Compute a simulated sample path. "        
     sys = self.ma_poly, self.ar_poly, 1
     u = np.random.randn(ts_length, 1)
     vals = dlsim(sys, u)[1]
     return vals.flatten()

 def simulation(self, ts_length=90) :
     " Compute a simulated sample path. "        
     sys = self.num, self.den, 1
     u = np.random.randn(ts_length, 1)
     return dlsim(sys, u)[1]

 def RunFilter(self,x_vec,t_vec):
     den = signal.convolve(np.array([self.Ts1,1]),np.array([self.Ts2,1]));
     num = self.dc_gain;
     filter_obj = signal.cont2discrete((num,den),self.Ts,method='zoh');
     tout, x_filt= signal.dlsim(filter_obj,x_vec,t=t_vec);
     return x_filt;

#.........这里部分代码省略.........
        # Run simulation
        if 'writeDict' in kwords or 'plotDict' in kwords \
        or 'mpcCont' in kwords or 'gust' in  kwords:
            # Determine whether to write and/or plot
            if 'writeDict' in kwords and Settings.WriteOut == True:
                write = True
            
            if 'plotDict' in kwords and Settings.PlotOut == True:
                #plot = True
                raise NotImplementedError()
            
            # Check function inputs for discrete time system plotting/writing
            if t is None:
                out_samples = U.shape[0]
                stoptime = (out_samples - 1) * disSys._Ts
            else:
                stoptime = t[-1]
                out_samples = int(np.floor(stoptime / disSys._Ts)) + 1
        
            # Pre-build output arrays
            xout = np.zeros((out_samples, disSys.A.shape[0]))
            yout = np.zeros((out_samples, disSys.C.shape[0]))
            tout = np.linspace(0.0, stoptime, num=out_samples)
        
            # Check initial condition
            if args[3] is None:
                xout[0,:] = np.zeros((disSys.A.shape[1],))
            else:
                xout[0,:] = np.asarray(args[3])
        
            # Pre-interpolate inputs into the desired time steps
            if t is None and U[0,0] is not None:
                u_dt = U
            elif U[0,0] is None and 'mpcCont' in kwords:
                u_dt = np.zeros((t.shape[0],disSys.nU))
            else:
                if len(U.shape) == 1:
                    U = U[:, np.newaxis]
        
                u_dt_interp = interp1d(t, U.transpose(),
                                       copy=False,
                                       bounds_error=True)
                u_dt = u_dt_interp(tout).transpose()
                
            if write == True:
                # Write output file header
                outputIndices = list(writeDict.values())
                ofile = Settings.OutputDir + \
                        Settings.OutputFileRoot + \
                        '_SOL302_out.dat'
                fp = open(ofile,'w')
                fp.write("{:<14}".format("Time"))
                for output in writeDict.keys():
                    fp.write("{:<14}".format(output))
                fp.write("\n")
                fp.flush()
            # END if write
                
            # Simulate the system
            for i in range(0, out_samples - 1):
                
                # get optimal control action
                if 'mpcCont' in kwords:             
                    u_dt[i,:kwords['mpcCont'].mpcU] = kwords['mpcCont'].getUopt(xout[i],True)
                    
                # get gust velocity at current time step
                if 'gust' in kwords:
                    raise NotImplementedError()
                
                xout[i+1,:] = np.dot(disSys.A, xout[i,:]) + np.dot(disSys.B, u_dt[i,:])
                yout[i,:] = np.dot(disSys.C, xout[i,:]) + np.dot(disSys.D, u_dt[i,:])
                
                if write == True:
                    fp.write("{:<14,e}".format(tout[i]))
                    for j in yout[i,outputIndices]:
                        fp.write("{:<14,e}".format(j))
                    fp.write("\n")
                    fp.flush()
        
            # Last point
            yout[out_samples-1,:] = np.dot(disSys.C, xout[out_samples-1,:]) + \
                                    np.dot(disSys.D, u_dt[out_samples-1,:])
            # Write final output and close
            fp.write("{:<14,e}".format(tout[out_samples-1]))
            for j in yout[out_samples-1,outputIndices]:
                fp.write("{:<14,e}".format(j))      
            fp.close()
                                    
            return tout, yout, xout
        else:
            # Run discrete time sim using scipy solver
            tout, yout, xout = dlsim((disSys.A,disSys.B,
                                      disSys.C,disSys.D,
                                      disSys._Ts),
                                      U,
                                      t = args[2],
                                      x0 = args[3])
        return tout, yout, xout
    else:
        raise ValueError("Needs 4 positional arguments")

import scipy.io.wavfile as wavfile


D = 500

f, data = wavfile.read('sp04.wav')
t = np.arange(0, len(data) / f, 1 / f)


echoes = np.concatenate((data[:D], data[D:] + 0.5 * data[:-D]))
wavfile.write('echoes.wav', f, echoes)

num = np.zeros(1)
num[0] = 1.0
den = np.zeros(D + 1)
den[0] = 1.0

for a in [0.5, 0.9, 0.25]:
    den[-1] = -a

    tf = (num, den, 1 / f)
    _, y = signal.dlsim(tf, echoes, t=t)
    wavfile.write('q2_minus_{0:.2f}.wav'.format(a), f, y)

for a in [0.5, 0.9, 0.25]:
    den[-1] = a

    tf = (num, den, 1 / f)
    _, y = signal.dlsim(tf, echoes, t=t)
    wavfile.write('q2_plus_{0:.2f}.wav'.format(a), f, y)

#----------------------------------------------------------------------#

# Take parameters as given
uar_alpha1 = .5
uar_alpha2 = 0.3
uar_alpha3 = -0.5
uar_alpha4 = -0.3

# Define the initial condition, numerator and denominator
uar_x0 = 10.
uar_num = [1.]
uar_den = [1., -uar_alpha1, -uar_alpha2, -uar_alpha3, -uar_alpha4]

uar_simulations = np.zeros((num_sims, len_sims))
uar_sys = (uar_num, uar_den, time_unit)

# Use num and den to do simulation
for i in xrange(num_sims):
    # Draw the epsilon shocks
    uar_eps = np.random.randn(len_sims-1, 1)

    # this stores the discrete impluse response
    t_out, y = sig.dlsim(uar_sys, uar_eps, x0=uar_x0)

    # store x0 in the first position and from position 1:T-1 store the computed impulse repsonse
    uar_simulations[i,0],uar_simulations[i,1:]=uar_x0,y.T

plt.plot(y)
plt.plot(uar_eps)
plt.show()

 def SimulateDynamics(self):
     t_out_x, y_out_x, x_out_x = signal.dlsim(self.SysX_d, self.SimulationParameters.u_vec, t=self.SimulationParameters.t_vec);
     t_out_y, y_out_y, x_out_y = signal.dlsim(self.SysY_d, self.SimulationParameters.u_vec, t=self.SimulationParameters.t_vec);
     t_out_z, y_out_z, x_out_z = signal.dlsim(self.SysZ_d, self.SimulationParameters.u_vec, t=self.SimulationParameters.t_vec);
     self.SimOL = SimulationResults(t_out_x,t_out_y,t_out_z,y_out_x,y_out_y,y_out_z,x_out_x,x_out_y,x_out_z);

 def test_linear_invariant(self, sys, sample_time, samples_number):
     time = np.arange(samples_number) * sample_time
     _, yout_cont, _ = lsim2(sys, T=time, U=time, **self.tolerances)
     _, yout_disc, _ = dlsim(c2d(sys, sample_time, method='foh'), u=time)
     assert_allclose(yout_cont.ravel(), yout_disc.ravel())