← Back to branch summary

~ubuntu-branches/ubuntu/raring/python-scipy/raring-proposed

« back to all changes in this revision

Viewing changes to Lib/signal/ltisys.py

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2007-01-07 14:12:12 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070107141212-mm0ebkh5b37hcpzn
http://bugs.debian.org/402783
* Remove build dependency on python-numpy-dev.
* python-scipy: Depend on python-numpy instead of python-numpy-dev.
* Package builds on other archs than i386. Closes: #402783.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/signal/ltisys.py
3
3
#
4
4
 
5
5
from filter_design import tf2zpk, zpk2tf, normalize
6
 
from Numeric import product, zeros, concatenate, \
7
 
     array, dot, transpose, arange, ones, Float
8
 
import Numeric
 
6
import numpy
 
7
from numpy import product, zeros, array, dot, transpose, arange, ones, \
 
8
    nan_to_num
9
9
import scipy.interpolate as interpolate
10
10
import scipy.integrate as integrate
11
11
import scipy.linalg as linalg
12
 
from scipy_base import r_, c_, eye, real, atleast_1d, atleast_2d, poly, \
 
12
from numpy import r_, eye, real, atleast_1d, atleast_2d, poly, \
13
13
     squeeze, diag, asarray
14
 
from Matrix import Matrix as Mat
15
14
 
16
15
def tf2ss(num, den):
17
16
    """Transfer function to state-space representation.
25
24
      A, B, C, D -- state space representation of the system.
26
25
    """
27
26
    # Controller canonical state-space representation.
28
 
    #  if M+1 = len(num) and K+1 = len(den) then we must have M <= K 
 
27
    #  if M+1 = len(num) and K+1 = len(den) then we must have M <= K
29
28
    #  states are found by asserting that X(s) = U(s) / D(s)
30
 
    #  then Y(s) = N(s) * X(s) 
 
29
    #  then Y(s) = N(s) * X(s)
31
30
    #
32
31
    #   A, B, C, and D follow quite naturally.
33
32
    #
34
33
    num, den = normalize(num, den)   # Strips zeros, checks arrays
35
34
    nn = len(num.shape)
36
35
    if nn == 1:
37
 
        num = asarray([num],num.typecode())
38
 
    M = num.shape[1] 
39
 
    K = len(den) 
 
36
        num = asarray([num], num.dtype)
 
37
    M = num.shape[1]
 
38
    K = len(den)
40
39
    if (M > K):
41
40
        raise ValueError, "Improper transfer function."
42
41
    if (M == 0 or K == 0):  # Null system
43
 
        return array([],Float), array([], Float), array([], Float), \
44
 
               array([], Float)
 
42
        return array([],float), array([], float), array([], float), \
 
43
               array([], float)
45
44
 
46
45
    # pad numerator to have same number of columns has denominator
47
 
    num = c_[zeros((num.shape[0],K-M),num.typecode()), num]
 
46
    num = r_['-1',zeros((num.shape[0],K-M), num.dtype), num]
48
47
 
49
48
    if num.shape[-1] > 0:
50
49
        D = num[:,0]
51
50
    else:
52
 
        D = array([],Float)
 
51
        D = array([],float)
53
52
 
54
53
    if K == 1:
55
 
        return array([], Float), array([], Float), array([], Float), D
56
 
    
 
54
        return array([], float), array([], float), array([], float), D
 
55
 
57
56
    frow = -array([den[1:]])
58
57
    A = r_[frow, eye(K-2, K-1)]
59
58
    B = eye(K-1, 1)
103
102
    if MD != MC:
104
103
        raise ValueError, "C and D must have the same number of rows."
105
104
    if ND != NB:
106
 
        raise ValueErrro, "B and D must have the same number of columns."
 
105
        raise ValueError, "B and D must have the same number of columns."
107
106
 
108
107
    return A, B, C, D
109
 
    
 
108
 
110
109
def ss2tf(A, B, C, D, input=0):
111
110
    """State-space to transfer function.
112
111
 
122
121
    """
123
122
    # transfer function is C (sI - A)**(-1) B + D
124
123
    A, B, C, D = map(asarray, (A, B, C, D))
125
 
    # Check consistency and 
 
124
    # Check consistency and
126
125
    #     make them all rank-2 arrays
127
126
    A, B, C, D = abcd_normalize(A, B, C, D)
128
127
 
139
138
 
140
139
    den = poly(A)
141
140
 
142
 
    if (product(B.shape) == 0) and (product(C.shape) == 0):
143
 
        num = Numeric.ravel(D)
144
 
        if (product(D.shape) == 0) and (product(A.shape) == 0):
 
141
    if (product(B.shape,axis=0) == 0) and (product(C.shape,axis=0) == 0):
 
142
        num = numpy.ravel(D)
 
143
        if (product(D.shape,axis=0) == 0) and (product(A.shape,axis=0) == 0):
145
144
            den = []
146
145
        end
147
146
        return num, den
148
147
 
149
 
    num_states = A.shape[0] 
 
148
    num_states = A.shape[0]
150
149
    type_test = A[:,0] + B[:,0] + C[0,:] + D
151
 
    num = Numeric.zeros((nout, num_states+1),type_test.typecode())
 
150
    num = numpy.zeros((nout, num_states+1), type_test.dtype)
152
151
    for k in range(nout):
153
152
        Ck = atleast_2d(C[k,:])
154
153
        num[k] = poly(A - dot(B,Ck)) + (D[k]-1)*den
179
178
    Outputs:
180
179
 
181
180
      z, p, k -- zeros and poles in sequences and gain constant.
182
 
    """    
 
181
    """
183
182
    return tf2zpk(*ss2tf(A,B,C,D,input=input))
184
183
 
185
 
class lti:
 
184
class lti(object):
186
185
    """Linear Time Invariant class which simplifies representation.
187
186
    """
188
187
    def __init__(self,*args,**kwords):
223
222
            self.__dict__['zeros'], self.__dict__['poles'], \
224
223
                                    self.__dict__['gain'] = ss2zpk(*args)
225
224
            self.__dict__['num'], self.__dict__['den'] = ss2tf(*args)
226
 
            self.inputs = B.shape[-1]
227
 
            self.outputs = C.shape[0]
 
225
            self.inputs = self.B.shape[-1]
 
226
            self.outputs = self.C.shape[0]
228
227
        else:
229
228
            raise ValueError, "Needs 2, 3, or 4 arguments."
230
229
 
256
255
                                           self.C, self.D)
257
256
            self.__dict__['num'], self.__dict__['den'] = \
258
257
                                  ss2tf(self.A, self.B,
259
 
                                        self.C, self.D)            
 
258
                                        self.C, self.D)
260
259
        else:
261
260
            self.__dict__[attr] = val
262
261
 
283
282
                  4 (A, B, C, D)
284
283
      U -- an input array describing the input at each time T
285
284
           (linear interpolation is assumed between given times).
286
 
           If there are multiple inputs, then each column of the 
 
285
           If there are multiple inputs, then each column of the
287
286
           rank-2 array represents an input.
288
287
      T -- the time steps at which the input is defined and at which
289
288
           the output is desired.
306
305
    #   being the number of inputs
307
306
 
308
307
    # rather than use lsim, use direct integration and matrix-exponential.
309
 
    if isinstance(system, lti):        
 
308
    if isinstance(system, lti):
310
309
        sys = system
311
310
    else:
312
311
        sys = lti(*system)
313
312
    U = atleast_1d(U)
314
313
    T = atleast_1d(T)
315
314
    if len(U.shape) == 1:
316
 
        U.shape = (U.shape[0],1)
317
 
    sU = U.shape        
 
315
        U = U.reshape((U.shape[0],1))
 
316
    sU = U.shape
318
317
    if len(T.shape) != 1:
319
318
        raise ValueError, "T must be a rank-1 array."
320
319
    if sU[0] != len(T):
323
322
        raise ValueError, "System does not define that many inputs."
324
323
 
325
324
    if X0 is None:
326
 
        X0 = zeros(sys.B.shape[0],sys.A.typecode())
 
325
        X0 = zeros(sys.B.shape[0],sys.A.dtype)
327
326
 
328
327
    # for each output point directly integrate assume zero-order hold
329
328
    #   or linear interpolation.
330
329
 
331
 
    ufunc = interpolate.linear_1d(T, U, axis=0, bounds_error=0, fill_value=0)
 
330
    ufunc = interpolate.interp1d(T, U, kind='linear', axis=0, bounds_error=False)
332
331
 
333
332
    def fprime(x, t, sys, ufunc):
334
 
        return dot(sys.A,x) + squeeze(dot(sys.B,ufunc([t])))
 
333
        return dot(sys.A,x) + squeeze(dot(sys.B,nan_to_num(ufunc([t]))))
335
334
 
336
335
    xout = integrate.odeint(fprime, X0, T, args=(sys, ufunc))
337
336
    yout = dot(sys.C,transpose(xout)) + dot(sys.D,transpose(U))
351
350
                  4 (A, B, C, D)
352
351
      U -- an input array describing the input at each time T
353
352
           (interpolation is assumed between given times).
354
 
           If there are multiple inputs, then each column of the 
 
353
           If there are multiple inputs, then each column of the
355
354
           rank-2 array represents an input.
356
355
      T -- the time steps at which the input is defined and at which
357
356
           the output is desired.
380
379
    U = atleast_1d(U)
381
380
    T = atleast_1d(T)
382
381
    if len(U.shape) == 1:
383
 
        U.shape = (U.shape[0],1)
384
 
    sU = U.shape        
 
382
        U = U.reshape((U.shape[0],1))
 
383
    sU = U.shape
385
384
    if len(T.shape) != 1:
386
385
        raise ValueError, "T must be a rank-1 array."
387
386
    if sU[0] != len(T):
390
389
        raise ValueError, "System does not define that many inputs."
391
390
 
392
391
    if X0 is None:
393
 
        X0 = zeros(sys.B.shape[0],sys.A.typecode())
 
392
        X0 = zeros(sys.B.shape[0], sys.A.dtype)
394
393
 
395
 
    xout = zeros((len(T),sys.B.shape[0]),sys.A.typecode())
 
394
    xout = zeros((len(T),sys.B.shape[0]), sys.A.dtype)
396
395
    xout[0] = X0
397
396
    A = sys.A
398
397
    AT, BT = transpose(sys.A), transpose(sys.B)
400
399
    lam, v = linalg.eig(A)
401
400
    vt = transpose(v)
402
401
    vti = linalg.inv(vt)
403
 
    GT = dot(dot(vti,diag(Numeric.exp(dt*lam))),vt).astype(xout.typecode())
 
402
    GT = dot(dot(vti,diag(numpy.exp(dt*lam))),vt).astype(xout.dtype)
404
403
    ATm1 = linalg.inv(AT)
405
404
    ATm2 = dot(ATm1,ATm1)
406
 
    I = eye(A.shape[0],typecode=A.typecode())
 
405
    I = eye(A.shape[0],dtype=A.dtype)
407
406
    GTmI = GT-I
408
407
    F1T = dot(dot(BT,GTmI),ATm1)
409
408
    if interp:
413
412
        dt1 = T[k] - T[k-1]
414
413
        if dt1 != dt:
415
414
            dt = dt1
416
 
            GT = dot(dot(vti,diag(Numeric.exp(dt*lam))),vt).astype(xout.typecode())
 
415
            GT = dot(dot(vti,diag(numpy.exp(dt*lam))),vt).astype(xout.dtype)
417
416
            GTmI = GT-I
418
417
            F1T = dot(dot(BT,GTmI),ATm1)
419
418
            if interp:
458
457
        vals = linalg.eigvals(sys.A)
459
458
        tc = 1.0/min(abs(real(vals)))
460
459
        T = arange(0,7*tc,7*tc / float(N))
461
 
    h = zeros(T.shape, sys.A.typecode())
 
460
    h = zeros(T.shape, sys.A.dtype)
462
461
    s,v = linalg.eig(sys.A)
463
462
    vi = linalg.inv(v)
464
463
    C = sys.C
465
464
    for k in range(len(h)):
466
 
        es = diag(Numeric.exp(s*T[k]))
467
 
        eA = (dot(dot(v,es),vi)).astype(h.typecode())
 
465
        es = diag(numpy.exp(s*T[k]))
 
466
        eA = (dot(dot(v,es),vi)).astype(h.dtype)
468
467
        h[k] = squeeze(dot(dot(C,eA),B))
469
468
    return T, h
470
469
 
495
494
        vals = linalg.eigvals(sys.A)
496
495
        tc = 1.0/min(abs(real(vals)))
497
496
        T = arange(0,7*tc,7*tc / float(N))
498
 
    U = ones(T.shape, sys.A.typecode())
 
497
    U = ones(T.shape, sys.A.dtype)
499
498
    vals = lsim(sys, U, T, X0=X0)
500
499
    return vals[0], vals[1]
501