# 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

- DELAY_TF
- GZP2PID
- HOLD_ON_TF
- MDS_DELTA_RESP
- MDS_SET2PAR
- MDS_STEP_RESP
- PID2GZP
- PLOT_AMP
- PLOT_BODE
- PLOT_PHASE
- POLY_MULT
- POLY_POW
- POLY_SUM
- QP_DYN_SOLVE
- RECURSIVE_DF
- TUSTIN
- UNWRAP_PHASE
- ZERO2COEFF

# 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)