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.
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)
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)
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
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, 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
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
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)
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)
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)
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)
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
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
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
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
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)
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
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
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)