OAA_LIB-tfl_lib Library Help

This page was created by the IDL library routine mk_html_help. For more information on this routine, refer to the IDL Online Help Navigator or type:

     ? mk_html_help

at the IDL command line prompt.

Last modified: Tue May 15 16:15:50 2007.


List of Routines


Routine Descriptions

DELAY_TF

[Next Routine] [List of Routines]

 DELAY_TF

 Complex transfer function corresponding to a pure delay time_delay
 at the frequency values stored in the freq_vec vector

 tf = delay_tf(freq_vec, time_delay)

 MODIFCATION HISTORY:

   Written by: A. Riccardi, Osservatorio Astrofisico di Arcetri, Italy

(See tfl_lib\delay_tf.pro)


GZP2PID

[Previous Routine] [Next Routine] [List of Routines]
 GPZ2PID

 convert a filter parametrized as gain-zero-pole into a PID

        (s+z[0])*...*(s+z[nz-1])               1         A
 gain * ------------------------  ==> kp + ki*--- + kd*-----*s
        (s+p[0])*...*(p+p[np-1])               s       (s+A)

 where z and p are the zeros and poles of the fileter.

 In this frame a valid PID is considered having:
    kp,ki,kd >= 0 and A > 0
 hence the constrains on gain, z, and p are:
    gain >= 0 and (nz = np or (nz = 0 and np = 1))
    if nz eq 0 and np eq 1: the pole must be real and not-negative;
    if nz eq 1: the zero and the pole must be real and not-negative;
    if nz eq 2: one pole is zero and the other real and strictly
                positive. The zeros complex conjugated (real
                part not-negative) or real and strictly positive;
    the gain can be zero only if nz = np = 0.



 err = gpz2pid(gain, z, p, kp, ki, kd, A, N_ZEROS=nz, N_POLES=np)


 gain:    real scalar. Filer Gain (not-negative). If gain eq 0.0 then
          the number of zeros and poles is forced to be 0.
 z:       real or complex vector. Vector of zeros. If N_ZERO is
          defined only the first N_ZERO elements are considered.
          z must have 0 (unefined=no zeros), 1 or 2 valid elements.
 p:       real vector. Vector of poles. If N_POLES is
          defined only the first N_POLES elements are considered.
          z must have the same number of valid elements as z.


 err eq 0 => the filter defined by gain, z, p is a valid PID
 err ne 0 => the filter is not a valid PID

 if N_ZEROS is set, only the first nz zeros in the vector z are used.
            For a valid PID nz can be 0, 1 or 2
            (set N_ZEROS to 0 if the filter has not zeros).
 if N_POLES is set, only the first np poles in the vector z are used.
            For a valid PID nz can be 0, 1 or 2
            (set N_POLES to 0 if the filter has not poles).

 NOTE: the GZP to PID conversion work well if poles are different
       from any zero and vice versa.

 MODIFICATON HISTORY

    Written by: A. Riccardi, Osservatorio Astrofisico di Arcetri, ITALY
                riccardi@arcetri.astro.it

(See tfl_lib\gzp2pid.pro)


HOLD_ON_TF

[Previous Routine] [Next Routine] [List of Routines]
 
 NAME

   HOLD_ON_TF


   tf = hold_on_tf(s, T)

 Return the (Laplace) Tranfer Function associated to the Hold-on
 (DAC) effect for a digital signal with period T.
 s vector of complex frequency [rad/s] (s = sigma + i*omega)
 T period [s]

            1 - exp(-s*T)      sinh(s*T/2)
   TF(s) = ---------------- = ----------- * T*exp(-s*T/2)
                  s               s*T/2

 MODIFICATON HISTORY

    Written by: A. Riccardi, Osservatorio Astrofisico di Arcetri, ITALY
                riccardi@arcetri.astro.it

(See tfl_lib\hold_on_tf.pro)


MDS_DELTA_RESP

[Previous Routine] [Next Routine] [List of Routines]
 MDS_DELTA_RESP

   responce of a mass-damper-spring system to a delta . The system is
   described by the laplace transform:

    (1-s*dT)*w0^2/(s^2+2*gamma*w0*s+w0^2)

   for gamma<1 the returned responce is giben by

    (w0^2/wd)*exp(-gamma*w0*t)*[sin(wd*t)+dT*w0*cos(wd*t+theta)]

   where wd=w0*sqrt(1-gamma^2) and theta=arcsin(gamma).

   For gamma>1 the responce is

    [1/(t1-t2)]*[(dT+t1)/t1*exp(-t/t1)-(dT+t2)/t2*exp(-t/t2)]

   where t1=[gamma+sqrt(gamma^2-1)]/w0 and t2=[gamma-sqrt(gamma^2-1)]/w0

   For gamma==1 the responce is

    exp(-t/t1)*[-dT*t1+t*(T+t1)]/t1^3

   where t1=1/w0

 responce = mds_step_resp(t, gamma, w0)

   considering w0=sqrt(k/m) and gamma=c/[2*sqrt(k*m)] if k->0 then w0=0
   and gamma->infinity. To handle this case set the keyword ZERO_K, and
   gamma=c/m (w0 is not considered)

   in this case the laplace transform is:

   (s*dT+1)/[s*(s+c/m)]

   and the step responce:

   1-exp(-t/t1)*(dT+t1)/t1

   where t1=m/c

   TO BE WRITTEN case ZERO_K

(See tfl_lib\mds_delta_resp.pro)


MDS_SET2PAR

[Previous Routine] [Next Routine] [List of Routines]
 mds_set2par, settiling_time, overshoot, w0, norm_damping, SETTLING_LIMIT=settling_level

 example:
 computes the natural frequency w0 [rad/sec] and normalized damping gamma.

 EXAMPLE:
 for a mass-damper-spring system with 1ms of settling time and 10% overshoot.

 mds_set2par, 0.001, 0.1, w0, gamma

 KEYWORDS:
 SETTLING_LIMIT: Threshold relative to the command for the computation of the settling
                 time. It is 0.1 by default (i.e. +/- 10% wrt command).It must be greater
                 then or equal to the overshoot

 NOTE: the equivalent delay for small w (in transfer function terms) is 2*gamma/w0.
       It is about 0.5*settling_time with 0.1 overshoot and 0.1 setting level.

 HISTORY:
  May 2006: written by A. Riccardi (AR). INAF-OAA

(See tfl_lib\mds_set2par.pro)


MDS_STEP_RESP

[Previous Routine] [Next Routine] [List of Routines]
 MDS_STEP_RESP

  resp = mds_step_resp(t, gamma, w0, ZERO_K=zero_k, DELAY=dT)

   unit step responce of a mass-damper-spring system described by the
   laplace transform:

    (1-s*dT)*w0^2/(s^2+2*gamma*w0*s+w0^2)

   for gamma<1 the returned responce is giben by

    1-(w0/wd)*exp(-gamma*w0*t)*[cos(wd*t-theta)+dT*w0*sin(wd*t)]

   where wd=w0*sqrt(1-gamma^2) and theta=arcsin(gamma).

   For gamma>1 the responce is

    1-[1/(t1-t2)]*[(dT+t1)*exp(-t/t1)-(dT+t2)*exp(-t/t2)]

   where t1=[gamma+sqrt(gamma^2-1)]/w0 and t2=[gamma-sqrt(gamma^2-1)]/w0

   For gamma==1 the responce is

    1-exp(-t/t1)*[t1^2+t*(T+t1)]/t1^2

   where t1=1/w0

 responce = mds_step_resp(t, gamma, w0)

   considering w0=sqrt(k/m) and gamma=c/[2*sqrt(k*m)] if k->0 then w0=0
   and gamma->infinity. To handle this case set the keyword ZERO_K, and
   gamma=c/m (w0 is not considered)

   in this case the laplace transform is:

   (s*dT+1)/[s*(s+c/m)]

   and the step responce:

   [t-dT-t1+exp(-t/t1)*(dT+t1)]/t1

   where t1=m/c

   TO BE WRITTEN case ZERO_K

(See tfl_lib\mds_step_resp.pro)


PID2GZP

[Previous Routine] [Next Routine] [List of Routines]
 PID2GZP

 convert a filter parametrized as a PID into a gain-zero-pole


          1         A                  (s+z[0])*...*(s+z[nz-1])
 kp + ki*--- + kd*-----*s  ==>  gain * ------------------------
          s       (s+A)                (s+p[0])*...*(p+z[np-1])

 pid2gzp, kp, ki, kd, A, gain, z, p, N_ZEROS=nz, N_POLES=np

 For a tipical PID:   A >> kp/kd >> ki/kp

 MODIFICATON HISTORY

    Written by: A. Riccardi, Osservatorio Astrofisico di Arcetri, ITALY
                riccardi@arcetri.astro.it

(See tfl_lib\pid2gzp.pro)


PLOT_AMP

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
    PLOT_AMP

 PURPOSE:

    Plot_amp plots (or overplots) the amplitude of a complex
    transfer function. The plot is log-log and the gridding is
    enabled.

 CATEGORY:

    Plotting Routines, Digital Filtering

 CALLING SEQUENCE:

    plot_amp, f_vec, complex_tf[, /DB|AUNITS=str][, /OVERPLOT|/COMPARISON]

 INPUTS:

    f_vec:       real vector. Vector of frequencyes. Frequencies less
                 or equal to zero are not considered.
    complex_tf:  complex vector. Transfer function. The number of
                 elements of complex_tf must be the same as f_vec.

 OPTIONAL INPUTS:

    None.

 KEYWORD PARAMETERS:

    DB:          If set, plot the amplitude axis in deciBel (dB)
                 units.

    AUNITS:      String containing the units of the amplitude axis.
                 It is not considered if the DB keyword is used.
    OVERPLOT:    If set, plot_amp overplots instead of plotting.

    COMPARISON:  Keyword used by plot_bode.pro. It is not considered if
                 the OVERPLOT keyword is used

    All the keywords allowed in plot (or overplot if OVERPLOT is set)
    can be added to the calling sequence.

 OUTPUTS:

    None.

 OPTIONAL OUTPUTS:

    None.

 COMMON BLOCKS:

    plot_amp_block. Just for internal use.

 SIDE EFFECTS:



 RESTRICTIONS:



 PROCEDURE:



 EXAMPLE:



 MODIFICATION HISTORY:

       Nov 1998, written by A. Riccardi 
                 Osservatorio Astrofisico di Arcetri, ITALY

       G. Brusa, Added COMPARISON and AUNITS keywords

(See tfl_lib\plot_amp.pro)


PLOT_BODE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
    PLOT_BODE

 PURPOSE:

    Plot_bode plots (or overplots) the amplitude and phase of a complex
    transfer function. See amp_plot and phase_plot for more details.

 CATEGORY:

    Plotting Routines, Digital Filtering

 CALLING SEQUENCE:

    plot_bode, f_vec, complex_tf[, AUNITS=astr][, AYRANGE=pyr] $
            [, PYRANGE=pyr][, /NOUNWRAP][, /COMPARISON]

 INPUTS:

    f_vec:       real vector. Vector of frequencyes. Frequencies less
                 or equal to zero are not considered.
    complex_tf:  complex vector. Transfer function. The number of
                 elements of complex_tf must be the same as f_vec.


 KEYWORD PARAMETERS:

    DB:          If set, plot the amplitude axis in deciBel (dB)
                 units.

    AUNITS:      String containing the units of the amplitude axis.
                 It is not considered if the DB keyword is used.

    AYRANGE:     Y-axis range in the amplitude plot.

    PYRANGE:     Y-axis range in the phase plot.

    COMPARISON:  Set this keyword to overplot the bode plots over a
                 plot generated by a previous plot_bode command.
                 The previous keywords are not considered if this
                 keyword is set.

    SMOOTH:      set this keyword to a smoothing window size to plot
                 smoothed data

    All the keywords allowed in plot (or overplot if COMPARISON is set)
    can be added to the calling sequence.

 COMMON BLOCKS:

    plot_amp_block. Just for internal use.

 MODIFICATION HISTORY:

       Nov 1998, written by A. Riccardi 
                 Osservatorio Astrofisico di Arcetri, ITALY

       G. Brusa, Added COMPARISON and AUNITS keywords

       Jul 2006, AR: added keyword SMOOTH

(See tfl_lib\plot_bode.pro)


PLOT_PHASE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
    PLOT_PHASE

 PURPOSE:

    Plot_phase plots (or overplots) the phase of a complex
    transfer function. The plot is log in the freq. axis and the
    gridding is enabled. The procedure unwrap the phase by default.

 CATEGORY:

    Plotting Routines, Digital Filtering

 CALLING SEQUENCE:

    plot_phase, f_vec, complex_tf [, /RAD][, FREQ_UNITS=str][, /OVERPLOT|/COMPARISON]
              [, /NO_UNWRAP]

 INPUTS:

    f_vec:       real vector. Vector of frequencyes. Frequencies less
                 or equal to zero are not considered.
    complex_tf:  complex vector. Transfer function. The number of
                 elements of complex_tf must be the same as f_vec.

 OPTIONAL INPUTS:

    None.

 KEYWORD PARAMETERS:

    RADIANTS:    If set, the phase is plotted in radiants (degree by
                 default).

    FREQ_UNITS:  string. Units of the frequency vector. Default value: "Hz".

    NO_UNWRAP:   If set, the phase unwrapping is disabled.

    OVERPLOT:    If set, plot_phase overplots instead of plotting.

    COMPARISON:  Used by plot_bode.pro. Not considered if OVERPLOT is set.

    All the keywords allowed in plot (or overplot if OVERPLOT is set)
    can be added to the calling sequence.

 OUTPUTS:

    None.

 OPTIONAL OUTPUTS:

    None.

 COMMON BLOCKS:

    plot_phase_block. Just for internal use.

 SIDE EFFECTS:



 RESTRICTIONS:



 PROCEDURE:



 EXAMPLE:



 MODIFICATION HISTORY:

       Nov 1998, written by A. Riccardi (AR)
       Osservarorio Astrofisico di Arcetri (OAA), ITALY
       

       G. Brusa, OAA
       Added overplot and comparison features

       Mar 2002, AR
       Fix of the label of the angle axis when radiants are displayed.
       RADIANTS keyword added, DEG keyword suppressed
       From now degree unit is the default

       Apr 2002, AR
       FREQ_UNITS keyword added

       July 2003, AR
       unwrapping phase code moved to unwrap_phase indipendent procedure

(See tfl_lib\plot_phase.pro)


POLY_MULT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
   POLY_MULT

 p = poly_mult(p1, p2 [, DOUBLE=double])

 p1, p2 vector of coeff of polynimials:
           P1(x)=p1[0]+p1[1]*x+p2[2]*x^2+...+p1[n1]*x^n1
           P2(x)=p2[0]+p2[1]*x+p2[2]*x^2+...+p2[n2]*x^n2

 p      vector of the coeffs of the polynomial P1(x)*P2(x)

 if n1 < n2 poly_mult(p1,p2) is faster then poly_mult(p2,p1)

 MODIFICATON HISTORY

    Written by: A. Riccardi, Osservatorio Astrofisico di Arcetri, ITALY
                riccardi@arcetri.astro.it

(See tfl_lib\poly_mult.pro)


POLY_POW

[Previous Routine] [Next Routine] [List of Routines]
 POLY_POW

 coeff_out = poly_pow(coeff_in, e [, DOUBLE=double])

 returns the coeffs of the polynomial given by the e-th power
 of the polynomial of coeffs coeff_in

 MODIFICATON HISTORY

    Written by: A. Riccardi, Osservatorio Astrofisico di Arcetri, ITALY
                riccardi@arcetri.astro.it

(See tfl_lib\poly_pow.pro)


POLY_SUM

[Previous Routine] [Next Routine] [List of Routines]
 POLY_SUM

 coeff_sum = poly_sum(coeff1, coeff2)

 returns the coeffs of the polynomial given by the sum
 of the polynomial of coeffs coeff1 and coeff2

 MODIFICATON HISTORY

    Written by: A. Riccardi, Osservatorio Astrofisico di Arcetri, ITALY
                riccardi@arcetri.astro.it

(See tfl_lib\poly_sum.pro)


QP_DYN_SOLVE

[Previous Routine] [Next Routine] [List of Routines]
 QP_DYN_SOLVE

 tf = qp_dyn_solve(freq, s_res, r_mode, l_mode, inv_A, FK_mat, FC_mat)

 returns the modal transfer function of the complex modes defined by the equation:

 M_mat##d^2x/dt^2 + C_mat##dx/dt + K_mat##x = FK_mat##c + FC_mat##dc/dt

 s_res, r_mode, l_mode and inv_A are returned by the qp_eigenproblem function:

 s_res = qp_eigenproblem(M_mat, C_mat, K_mat, r_mode, l_mode, inv_A)

 the system can be write as (see qp_eigenproblem for inv_A and R definition)

 dq/dt-R##q = inv_A##(FK##c+FK##dc/dt)
 where:
 q=[[x],[dx/dt]]
 FK=[[FK_mat],[zeros]]
 FC=[[FC_mat],[zeros]]

 q=r_mode##a (a= coeff column vector of decomposition of q over r_modes,
              see qp_eigenproblem for r_mode definition)

 a[i]=transpose(l_mode[i,*])##q/(transpose(l_mode[i,*])##r_mode[i,*]))

 in the Laplace (L{f}(s)=Laplace transform of f(t)) space:

 (s*I-R)##r_mode##L{a}(s) = inv_A##(FK_mat+s*FC_mat)##L{c}(s)

 multipling on the right with transpose(l_mode)

 (s-diag(s_res))##diag(h)##L{a} = transpose(r_mode)##inv_A##(FK_mat+s*FC_mat)##L{c}(s)
 where: h[i] = transpose(l_mode[i,*])##r_mode[i,*]

 N=n_elements(s_res)=2*n_elements(x)
 M=n_elements(c)

 we have NxM transfer functions: tf[s,i,j]=L{a[0,i]}(s)/L{c[0,j]}(s)

 qp_dyn_solve returns tf[f,i,j] where s=2*!PI*complex(0,1)

 Note: c cannot be expanded in the r_mode basis

 HISTORY:
   Jul 2003: written by A. Riccardi (AR), riccardi@arcetri.astro.it
   Jun 2006: AR, help written

(See tfl_lib\qp_dyn_solve.pro)


RECURSIVE_DF

[Previous Routine] [Next Routine] [List of Routines]
 RECURSIVE_DF

 filt_sig = recursive_df(sig, numd, dend)

 apply the digital recursive filter defined by numd, dend on the digital signal sig.
 The filtered signal is returned. See tustin.pro for a definition of numd and dend.

 MODIFICATON HISTORY

    Written by: A. Riccardi (AR), Osservatorio Astrofisico di Arcetri, ITALY
                riccardi@arcetri.astro.it
    15 Jul 2004, AR & M. Xompero
      bug when filter has no poles in z^-1 fixed

(See tfl_lib\recursive_df.pro)


TUSTIN

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
    TUSTIN

 PURPOSE:

 The Tustin (bilinear) transform is a mathematical mapping of variables.
 In digital filtering, it is a standard method of mapping the s
 (or analog) plane into the z (or digital) plane. It transforms
 analog filters, designed using classical filter design trchniques,
 into their discrete equivalents.

 The bilinear transformation maps the s-plane into the z-plane by:

           Hd(z) = Hc[s=2*fs*(z-1)*(z+1)] ,

 where Hc(s) is the analog filter (continous domain), Hd(z) if
 the discrete filter and fs is the frequency of the sampling.

 This transformation maps the j*W axis of the
 s=S+j*W plane into the unit circle in the z-plane by:

           w = 2*fs*atan(W/(2*fs)) ,

 where z=exp(i*w).

 Bilinear can accept an optional keyword PREWARP containing a
 frequency for which the frequency responces before and after mapping
 match exactly. In prewarped mode, the bilinear tranformation maps
 the s-plane into the z-plane with:

           Hd(z) = Hc[s=2*pi*fs/tan(pi*fp/fs)*(z-1)/(z+1)] .

 With the prewarping option, bilinear maps the j*W axis around the
 unit circle by:

           w = 2*atan[W*tan(pi*fp/fs)/(2*pi*fp)] .

 In prewarp mode, bilinear maches the frequency 2*pi*fp [rad/s] in
 the s-plane to the normalized frequency 2*pi*fp/fs [rad/s] in the z-plane.

 CATEGORY:
    Signal processing


 CALLING SEQUENCE:

    tustin, num_c, den_c, fs, num_d, den_d [, PREWARP=fp][, DOUBLE=double]
 
 INPUTS:

    num_c:      real vector. Vector of real coeffs of the numerator
                of the analog filter tranfer function (ascending
                powers of s).
    den_c:      real vector. Vector of real coeffs of the denominator
                of the analog filter tranfer function (ascending
                powers of s).
    fs:         real scalar. Sampling frequency [Hz].
      
 KEYWORD PARAMETERS:

    PREWARP:    real scalar. Optional. If defined contains the
                prewarping frequency [Hz]
    DOUBLE:     if set, force double precision computation.

 OUTPUTS:

    num_d:      real vector. Vector of real coeffs of the numerator
                of the discrete filter (ascending powers of z^-1).
    den_d:      real vector. Vector of real coeffs of the numerator
                of the discrete filter (ascending powers of z^-1).

 EXAMPLE:



 MODIFICATION HISTORY:

       Oct 1998, written by Armando Riccardi OAA 

(See tfl_lib\tustin.pro)


UNWRAP_PHASE

[Previous Routine] [Next Routine] [List of Routines]
 UNWRAP_PHASE

 unwrap_phase, phase

 unwrap phase vector (in radians). The unwrapped phase is overwritten

 HISTORY
   07-2003: written by A. Riccardi, Osservatorio di Arceri, ITALY
            riccardi@arcetri.astro.it

(See tfl_lib\unwrap_phase.pro)


ZERO2COEFF

[Previous Routine] [List of Routines]

 ZERO2COEFF

  c = zero2coeff(z)

  returns the coeffs of the polynomial:

    (x+z[0])*(x+z[1])* ... *(x+z[n-1])

  where n>0. z can be a real or a complex vector.

  The coeffs are ordered from the lowest to the highest power of x

 MODIFICATON HISTORY

    Written by: A. Riccardi, Osservatorio Astrofisico di Arcetri, ITALY
                riccardi@arcetri.astro.it

(See tfl_lib\zero2coeff.pro)