Python scipy.signal.dlti() Examples

def butter(order, Wn, N=1, btype='lowpass'):
    """
    build MIMO butterworth filter of order ord and cut-off freq over Nyquist
    freq ratio Wn.
    The filter will have N input and N output and N*ord states.

    Note: the state-space form of the digital filter does not depend on the
    sampling time, but only on the Wn ratio.
    As a result, this function only returns the A,B,C,D matrices of the filter
    state-space form.
    """

    # build DLTI SISO
    num, den = scsig.butter(order, Wn, btype=btype, analog=False, output='ba')
    Af, Bf, Cf, Df = scsig.tf2ss(num, den)
    SSf = scsig.dlti(Af, Bf, Cf, Df, dt=1.0)

    SStot = SSf
    for ii in range(1, N):
        SStot = join2(SStot, SSf)

    return SStot.A, SStot.B, SStot.C, SStot.D


# ----------------------------------------------------------------------- Utils 

def freqresp(self, wv):
        """
        Calculate frequency response over frequencies wv

        Note: this wraps frequency response function.
        """
        dlti = True
        if self.dt == None: dlti = False
        return freqresp(self, wv, dlti=dlti) 

def build_SS_poly(Acf, ds, negative=False):
    """
    Builds a discrete-time state-space representation of a polynomial system
    whose frequency response has from:
        Ypoly[oo,ii](k) = -A2[oo,ii] D2(k) - A1[oo,ii] D1(k) - A0[oo,ii]
    where C1,D2 are discrete-time models of first and second derivatives, ds is
    the time-step and the coefficient matrices are such that:
        A{nn}=Acf[oo,ii,nn]
    """

    Nout, Nin, Ncf = Acf.shape
    assert Ncf == 3, 'Acf input last dimension must be equal to 3!'

    Ader, Bder, Cder, Dder = SSderivative(ds)
    SSder = scsig.dlti(Ader, Bder, Cder, Dder, dt=ds)
    SSder02 = series(SSder, join2(np.array([[1]]), SSder))

    SSder_all = copy.deepcopy(SSder02)
    for ii in range(Nin - 1):
        SSder_all = join2(SSder_all, SSder02)

    # Build polynomial forcing terms
    sign = 1.0
    if negative == True: sign = -1.0

    A0 = Acf[:, :, 0]
    A1 = Acf[:, :, 1]
    A2 = Acf[:, :, 2]
    Kforce = np.zeros((Nout, 3 * Nin))
    for ii in range(Nin):
        Kforce[:, 3 * ii] = sign * (A0[:, ii])
        Kforce[:, 3 * ii + 1] = sign * (A1[:, ii])
        Kforce[:, 3 * ii + 2] = sign * (A2[:, ii])
    SSpoly_neg = addGain(SSder_all, Kforce, where='out')

    return SSpoly_neg 

def _test_phaseshift(self, method, zero_phase):
        rate = 120
        rates_to = [15, 20, 30, 40]  # q = 8, 6, 4, 3

        t_tot = int(100)  # Need to let antialiasing filters settle
        t = np.arange(rate*t_tot+1) / float(rate)

        # Sinusoids at 0.8*nyquist, windowed to avoid edge artifacts
        freqs = np.array(rates_to) * 0.8 / 2
        d = (np.exp(1j * 2 * np.pi * freqs[:, np.newaxis] * t)
             * signal.windows.tukey(t.size, 0.1))

        for rate_to in rates_to:
            q = rate // rate_to
            t_to = np.arange(rate_to*t_tot+1) / float(rate_to)
            d_tos = (np.exp(1j * 2 * np.pi * freqs[:, np.newaxis] * t_to)
                     * signal.windows.tukey(t_to.size, 0.1))

            # Set up downsampling filters, match v0.17 defaults
            if method == 'fir':
                n = 30
                system = signal.dlti(signal.firwin(n + 1, 1. / q,
                                                   window='hamming'), 1.)
            elif method == 'iir':
                n = 8
                wc = 0.8*np.pi/q
                system = signal.dlti(*signal.cheby1(n, 0.05, wc/np.pi))

            # Calculate expected phase response, as unit complex vector
            if zero_phase is False:
                _, h_resps = signal.freqz(system.num, system.den,
                                          freqs/rate*2*np.pi)
                h_resps /= np.abs(h_resps)
            else:
                h_resps = np.ones_like(freqs)

            y_resamps = signal.decimate(d.real, q, n, ftype=system,
                                        zero_phase=zero_phase)

            # Get phase from complex inner product, like CSD
            h_resamps = np.sum(d_tos.conj() * y_resamps, axis=-1)
            h_resamps /= np.abs(h_resamps)
            subnyq = freqs < 0.5*rate_to

            # Complex vectors should be aligned, only compare below nyquist
            assert_allclose(np.angle(h_resps.conj()*h_resamps)[subnyq], 0,
                            atol=1e-3, rtol=1e-3) 

def fast_xambg(refChannel, srvChannel, rangeBins, freqBins, shortFilt=True):
    ''' Fast Cross-Ambiguity Fuction (frequency domain method)
    
    Parameters:
        refChannel: reference channel data
        srvChannel: surveillance channel data
        rangeBins:  number of range bins to compute
        freqBins:   number of doppler bins to compute (should be power of 2)
        shortFilt:  (bool) chooses the type of decimation filter to use.
                    If True, uses an all-ones filter of length 1*(decimation factor)
                    If False, uses a flat-top window of length 10*(decimation factor)+1
    Returns:
        xambg: the cross-ambiguity surface. Dimensions are (nf, nlag+1, 1)
        third dimension added for easy stacking in dask

    '''
    if refChannel.shape != srvChannel.shape:
        raise ValueError('Input vectors must have the same length')

    # calculate decimation factor
    ndecim = int(refChannel.shape[0]/freqBins)

    # pre-allocate space for the result
    xambg = np.zeros((freqBins, rangeBins+1, 1), dtype=np.complex64)

    # complex conjugate of the second input vector
    srvChannelConj = np.conj(srvChannel)    

    if shortFilt:
        # precompute short FIR filter for decimation (all ones filter with length
        # equal to the decimation factor)
        dtaps = np.ones((ndecim + 1,))
    else:
        # precompute long FIR filter for decimation. (flat top filter of length
        # 10*decimation factor).  
        dtaps = signal.firwin(10*ndecim + 1, 1. / ndecim, window='flattop')

    dfilt = signal.dlti(dtaps, 1)

    # loop over range bins 
    for k, lag in enumerate(np.arange(-1*rangeBins, 1)):
        channelProduct = np.roll(srvChannelConj, lag)*refChannel
        #decimate the product of the reference channel and the delayed surveillance channel
        xambg[:,k,0] = signal.decimate(channelProduct, ndecim, ftype=dfilt)[0:freqBins]

    # take the FFT along the first axis (Doppler)
    # xambg = np.fft.fftshift(np.fft.fft(xambg, axis=0), axes=0)
    xambg = np.fft.fftshift(fft(xambg, axis=0), axes=0)
    return xambg 

def freqresp(SS, wv, dlti=True):
    """
    In-house frequency response function supporting dense/sparse types

    Inputs:
    - SS: instance of ss class, or scipy.signal.StateSpace*
    - wv: frequency range
    - dlti: True if discrete-time system is considered.

    Outputs:
    - Yfreq[outputs,inputs,len(wv)]: frequency response over wv

    Warnings:
    -  This function may not be very efficient for dense matrices (as A is not
    reduced to upper Hessenberg form), but can exploit sparsity in the state-space
    matrices.
    """

    assert type(SS) == ss, \
        'Type %s of state-space model not supported. Use libss.ss instead!' % type(SS)
    SS.check_types()

    if hasattr(SS, 'dt') and dlti:
        Ts = SS.dt
        wTs = Ts * wv
        zv = np.cos(wTs) + 1.j * np.sin(wTs)
    else:
        print('Assuming a continuous time system')
        zv = 1.j * wv

    Nx = SS.A.shape[0]
    Ny = SS.D.shape[0]
    try:
        Nu = SS.B.shape[1]
    except IndexError:
        Nu = 1

    Nw = len(wv)

    Yfreq = np.empty((Ny, Nu, Nw,), dtype=np.complex_)
    Eye = libsp.eye_as(SS.A)
    for ii in range(Nw):
        sol_cplx = libsp.solve(zv[ii] * Eye - SS.A, SS.B)
        Yfreq[:, :, ii] = libsp.dot(SS.C, sol_cplx, type_out=np.ndarray) + SS.D

    return Yfreq 

def parallel(SS01, SS02):
    """
    Returns the sum (or parallel connection of two systems). Given two state-space
    models with the same output, but different input:
        u1 --> SS01 --> y
        u2 --> SS02 --> y

    """

    if type(SS01) is not type(SS02):
        raise NameError('The two input systems need to have the same size!')
    if SS01.dt != SS02.dt:
        raise NameError('DLTI systems do not have the same time-step!')
    Nout = SS02.outputs
    if Nout != SS01.outputs:
        raise NameError('DLTI systems need to have the same number of output!')

    # if type(SS01) is control.statesp.StateSpace:
    # 	SStot=control.parallel(SS01,SS02)
    # else:

    # determine size of total system
    Nst01, Nst02 = SS01.states, SS02.states
    Nst = Nst01 + Nst02
    Nin01, Nin02 = SS01.inputs, SS02.inputs
    Nin = Nin01 + Nin02

    # Build A,B matrix
    A = np.zeros((Nst, Nst))
    A[:Nst01, :Nst01] = SS01.A
    A[Nst01:, Nst01:] = SS02.A
    B = np.zeros((Nst, Nin))
    B[:Nst01, :Nin01] = SS01.B
    B[Nst01:, Nin01:] = SS02.B

    # Build the rest
    C = np.block([SS01.C, SS02.C])
    D = np.block([SS01.D, SS02.D])

    SStot = scsig.dlti(A, B, C, D, dt=SS01.dt)

    return SStot 

def ss_to_scipy(ss):
    """
    Converts to a scipy.signal linear time invariant system

    Args:
        ss (libss.ss): SHARPy state space object

    Returns:
        scipy.signal.dlti
    """

    if ss.dt == None:
        sys = scsig.lti(ss.A, ss.B, ss.C, ss.D)
    else:
        sys = scsig.dlti(ss.A, ss.B, ss.C, ss.D, dt=ss.dt)

    return sys

# 1/0

# # check parallel connector
# Nout=2
# Nin01,Nin02=2,3
# Nst01,Nst02=4,2

# # build random systems
# fac=0.1
# A01,A02=fac*np.random.rand(Nst01,Nst01),fac*np.random.rand(Nst02,Nst02)
# B01,B02=np.random.rand(Nst01,Nin01),np.random.rand(Nst02,Nin02)
# C01,C02=np.random.rand(Nout,Nst01),np.random.rand(Nout,Nst02)
# D01,D02=np.random.rand(Nout,Nin01),np.random.rand(Nout,Nin02)

# dt=0.1
# SS01=scsig.StateSpace( A01,B01,C01,D01,dt=dt )
# SS02=scsig.StateSpace( A02,B02,C02,D02,dt=dt )

# # simulate
# NT=11
# U01,U02=np.random.rand(NT,Nin01),np.random.rand(NT,Nin02)

# # reference
# Y01,X01=simulate(SS01,U01)
# Y02,X02=simulate(SS02,U02)
# Yref=Y01+Y02

# # parallel
# SStot=parallel(SS01,SS02)
# Utot=np.block([U01,U02])
# Ytot,Xtot=simulate(SStot,Utot)

# # join method
# SStot=join(SS01,SS02)
# K=np.array([[1,2,3],[4,5,6]])
# SStot=join(K,SS02)
# SStot=join(SS02,K)
# K2=np.array([[10,20,30],[40,50,60]]).T
# Ktot=join(K,K2)