SB10YD

Fitting frequency response data with a stable, minimum phase SISO system

Purpose

  To fit a supplied frequency response data with a stable, minimum
  phase SISO (single-input single-output) system represented by its
  matrices A, B, C, D. It handles both discrete- and continuous-time
  cases.

      SUBROUTINE SB10YD( DISCFL, FLAG, LENDAT, RFRDAT, IFRDAT, OMEGA, N,
     $                   A, LDA, B, C, D, TOL, IWORK, DWORK, LDWORK,
     $                   ZWORK, LZWORK, INFO )
C     .. Scalar Arguments ..
      INTEGER            DISCFL, FLAG, INFO, LDA, LDWORK, LENDAT,
     $                   LZWORK, N
      DOUBLE PRECISION   TOL
C     .. Array Arguments ..
      INTEGER            IWORK(*)
      DOUBLE PRECISION   A(LDA, *), B(*), C(*), D(*), DWORK(*),
     $                   IFRDAT(*), OMEGA(*), RFRDAT(*)
      COMPLEX*16         ZWORK(*)

Input/Output Parameters

  DISCFL  (input) INTEGER
          Indicates the type of the system, as follows:
          = 0: continuous-time system;
          = 1: discrete-time system.

  FLAG    (input) INTEGER
          If FLAG = 0, then the system zeros and poles are not
          constrained.
          If FLAG = 1, then the system zeros and poles will have
          negative real parts in the continuous-time case, or moduli
          less than 1 in the discrete-time case. Consequently, FLAG
          must be equal to 1 in mu-synthesis routines.

  LENDAT  (input) INTEGER
          The length of the vectors RFRDAT, IFRDAT and OMEGA.
          LENDAT >= 2.

  RFRDAT  (input) DOUBLE PRECISION array, dimension (LENDAT)
          The real part of the frequency data to be fitted.

  IFRDAT  (input) DOUBLE PRECISION array, dimension (LENDAT)
          The imaginary part of the frequency data to be fitted.

  OMEGA   (input) DOUBLE PRECISION array, dimension (LENDAT)
          The frequencies corresponding to RFRDAT and IFRDAT.
          These values must be nonnegative and monotonically
          increasing. Additionally, for discrete-time systems
          they must be between 0 and PI.

  N       (input/output) INTEGER
          On entry, the desired order of the system to be fitted.
          N <= LENDAT-1.
          On exit, the order of the obtained system. The value of N
          could only be modified if N > 0 and FLAG = 1.

  A       (output) DOUBLE PRECISION array, dimension (LDA,N)
          The leading N-by-N part of this array contains the
          matrix A. If FLAG = 1, then A is in an upper Hessenberg
          form, and corresponds to a minimal realization.

  LDA     INTEGER
          The leading dimension of the array A.  LDA >= MAX(1,N).

  B       (output) DOUBLE PRECISION array, dimension (N)
          The computed vector B.

  C       (output) DOUBLE PRECISION array, dimension (N)
          The computed vector C. If FLAG = 1, the first N-1 elements
          are zero (for the exit value of N).

  D       (output) DOUBLE PRECISION array, dimension (1)
          The computed scalar D.

  TOL     DOUBLE PRECISION
          The tolerance to be used for determining the effective
          rank of matrices. If the user sets TOL > 0, then the given
          value of TOL is used as a lower bound for the reciprocal
          condition number;  a (sub)matrix whose estimated condition
          number is less than 1/TOL is considered to be of full
          rank.  If the user sets TOL <= 0, then an implicitly
          computed, default tolerance, defined by TOLDEF = SIZE*EPS,
          is used instead, where SIZE is the product of the matrix
          dimensions, and EPS is the machine precision (see LAPACK
          Library routine DLAMCH).

  IWORK   INTEGER array, dimension (max(2,2*N+1))

  DWORK   DOUBLE PRECISION array, dimension (LDWORK)
          On exit, if INFO = 0, DWORK(1) returns the optimal value
          of LDWORK and DWORK(2) contains the optimal value of
          LZWORK.

  LDWORK  INTEGER
          The length of the array DWORK.
          LDWORK = max( 2, LW1, LW2, LW3, LW4 ), where
          LW1 = 2*LENDAT + 4*HNPTS;  HNPTS = 2048;
          LW2 =   LENDAT + 6*HNPTS;
          MN  = min( 2*LENDAT, 2*N+1 )
          LW3 = 2*LENDAT*(2*N+1) + max( 2*LENDAT, 2*N+1 ) +
                max( MN + 6*N + 4, 2*MN + 1 ), if N > 0;
          LW3 = 4*LENDAT + 5                 , if N = 0;
          LW4 = max( N*N + 5*N, 6*N + 1 + min( 1,N ) ), if FLAG = 1;
          LW4 = 0,                                      if FLAG = 0.
          For optimum performance LDWORK should be larger.

  ZWORK   COMPLEX*16 array, dimension (LZWORK)

  LZWORK  INTEGER
          The length of the array ZWORK.
          LZWORK = LENDAT*(2*N+3), if N > 0;
          LZWORK = LENDAT,         if N = 0.

  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          = 1:  if the discrete --> continuous transformation cannot
                be made;
          = 2:  if the system poles cannot be found;
          = 3:  if the inverse system cannot be found, i.e., D is
                (close to) zero;
          = 4:  if the system zeros cannot be found;
          = 5:  if the state-space representation of the new
                transfer function T(s) cannot be found;
          = 6:  if the continuous --> discrete transformation cannot
                be made.

  First, if the given frequency data are corresponding to a
  continuous-time system, they are changed to a discrete-time
  system using a bilinear transformation with a scaled alpha.
  Then, the magnitude is obtained from the supplied data.
  Then, the frequency data are linearly interpolated around
  the unit-disc.
  Then, Oppenheim and Schafer complex cepstrum method is applied
  to get frequency data corresponding to a stable, minimum-
  phase system. This is done in the following steps:
  - Obtain LOG (magnitude)
  - Obtain IFFT of the result (DG01MD SLICOT subroutine);
  - halve the data at 0;
  - Obtain FFT of the halved data (DG01MD SLICOT subroutine);
  - Obtain EXP of the result.
  Then, the new frequency data are interpolated back to the
  original frequency.
  Then, based on these newly obtained data, the system matrices
  A, B, C, D are constructed; the very identification is
  performed by Least Squares Method using DGELSY LAPACK subroutine.
  If needed, a discrete-to-continuous time transformation is
  applied on the system matrices by AB04MD SLICOT subroutine.
  Finally, if requested, the poles and zeros of the system are
  checked. If some of them have positive real parts in the
  continuous-time case (or are not inside the unit disk in the
  complex plane in the discrete-time case), they are exchanged with
  their negatives (or reciprocals, respectively), to preserve the
  frequency response, while getting a minimum phase and stable
  system. This is done by SB10ZP SLICOT subroutine.

  [1] Oppenheim, A.V. and Schafer, R.W.
      Discrete-Time Signal Processing.
      Prentice-Hall Signal Processing Series, 1989.

  [2] Balas, G., Doyle, J., Glover, K., Packard, A., and Smith, R.
      Mu-analysis and Synthesis toolbox - User's Guide,
      The Mathworks Inc., Natick, MA, USA, 1998.

  None

Program Text

  None

  None

  None