% Copyright INRIA \chapter{Description of the Basic Tools} \section{Introduction} The purpose of this document is to illustrate the use of the Scilab software package in a signal processing context. We have gathered a collection of signal processing algorithms which have been implemented as Scilab functions. This manual is in part a pedagogical tool concerning the study of signal processing and in part a practical guide to using the signal processing tools available in Scilab. For those who are already well versed in the study of signal processing the tutorial parts of the manual will be of less interest. For each signal processing tool available in the signal processing toolbox there is a tutorial section in the manual explaining the methodology behind the technique. This section is followed by a section which describes the use of a function designed to accomplish the signal processing described in the preceding sections. At this point the reader is encouraged to launch a Scilab session and to consult the on-line help related to the function in order to get the precise and complete description (syntax, description of its functionality, examples and related functions). This section is in turn followed by an examples section demonstrating the use of the function. In general, the example section illustrates more clearly than the syntax section how to use the different modes of the function. In this manual the {\tt typewriter-face} font is used to indicate either a function name or an example dialogue which occurs in Scilab. Each signal processing subject is illustrated by examples and figures which were demonstrated using Scilab. To further assist the user, there exists for each example and figure an executable file which recreates the example or figure. To execute an example or figure one uses the following Scilab command \begin{verbatim} -->exec('file.name') \end{verbatim} which causes Scilab to execute all the Scilab commands contained in the file called {\tt file.name}. To know what signal processing tools are available in Scilab one would type \begin{verbatim} -->disp(siglib) \end{verbatim} \index{toolbox library} \index{library sigblib@library {\tt siglib}} which produces a list of all the signal processing functions available in the signal processing library. \section{Signals} \index{signals} For signal processing the first point to know is how to load and save signals or only small portions of lengthy signals that are to be used or are to be generated by Scilab. Finally, the generation of synthetic (random) signals is an important tool in the development in implementation of signal processing tools. This section addresses all of these topics. \subsection{Saving, Loading, Reading, and Writing Files} \index{signals!saving and loading} Signals and variables which have been processed or created in the Scilab environment can be saved in files written directly by Scilab. The syntax for the {\tt save} primitive is \begin{verbatim} -->save(file_name[,var_list]) \end{verbatim} \index{function syntax!save@{\tt save}} where {\tt file\_name} is the file to be written to and {\tt var\_list} is the list of variables to be written. The inverse to the operation {\tt save} is accomplished by the primitive {\tt load} which has the syntax \begin{verbatim} -->load(file_name[,var_list]) \end{verbatim} \index{function syntax!load@{\tt load}} where the argument list is identical that used in {\tt save}. Although the commands {\tt save} and {\tt load} are convenient, one has much more control over the transfer of data between files and Scilab by using the commands {\tt read} and {\tt write}. These two commands work similarly to the read and write commands found in Fortran. The syntax of these two commands is as follows. The syntax for {\tt write} is \begin{verbatim} -->write(file,x[,form]) \end{verbatim} \index{function syntax!write@{\tt write}} The second argument, {\tt x}, is a matrix of values which are to be written to the file. The syntax for {\tt read} is \begin{verbatim} -->x=read(file,m,n[,form]) \end{verbatim} \index{function syntax!read@{\tt read}} The arguments {\tt m} and {\tt n} are the row and column dimensions of the resulting data matrix {\tt x}. and {\tt form} is again the format specification statement. In order to illustrate the use of the on-line help for reading this manual we give the result of the Scilab command \begin{verbatim} -->help read \end{verbatim} \begin{verbatim} read(1) Scilab Function read(1) NAME read - matrices read CALLING SEQUENCE [x]=read(file-name,m,n,[format]) [x]=read(file-name,m,n,k,format) PARAMETERS file-name : string or integer (logical unit number) m, n : integers (dimensions of the matrix x). Set m=-1 if you dont know the numbers of rows, so the whole file is read. format : string (fortran format). If format='(a)' then read reads a vec- tor of strings n must be equal to 1. k : integer DESCRIPTION reads row after row the mxn matrix x (n=1 for character chain) in the file file-name (string or integer). Two examples for format are : (1x,e10.3,5x,3(f3.0)),(10x,a20) ( the default value is *). The type of the result will depend on the specified form. If form is numeric (d,e,f,g) the matrix will be a scalar matrix and if form contains the character a the matrix will be a matrix of character strings. A direct access file can be used if using the parameter k which is is the vector of record numbers to be read (one record per row), thus m must be m=prod(size(k)). To read on the keyboard use read(%io(1),...). EXAMPLE A=rand(3,5); write('foo',A); B=read('foo',3,5) B=read('foo',-1,5) read(%io(1),1,1,'(a)') // waits for user's input SEE ALSO file, readb, write, %io, x_dialog \end{verbatim} \subsection{Simulation of Random Signals} \index{signals!simulation of random} \index{simulation of random signals} \index{random number generator} The creation of synthetic signals can be accomplished using the Scilab function {\tt rand} which generates random numbers. The user can generate a sequence of random numbers, a random matrix with the uniform or the gaussian probability laws. A seed is possible to re-create the same pseudo-random sequences. Often it is of interest in signal processing to generate normally distributed random variables with a certain mean and covariance structure. This can be accomplished by using the standard normal random numbers generated by {\tt rand} and subsequently modifying them by performing certain linear numeric operations. For example, to obtain a random vector $y$ which is distributed N($m_y$,$\Lambda_y$) from a random vector $x$ which is distributed standard normal (i.e. N(0,I)) one would perform the following operation \begin{equation} y=\Lambda_y^{1/2}x+m_y \label{e.rand1} \end{equation} where $\Lambda_y^{1/2}$ is the matrix square root of $\Lambda_y$. A matrix square root can be obtained using the {\tt chol} primitive as follows \verbatok{\Diary rand1.dia} \end{verbatim} taking note that it is the transpose of the matrix obtained from {\tt chol} that is used for the square root of the desired covariance matrix. Sequences of random numbers following a specific normally distributed probability law can also be obtained by filtering. That is, a white standard normal sequence of random numbers is passed through a linear filter to obtain a normal sequence with a specific spectrum. For a filter which has a discrete Fourier transform $H(w)$ the resulting filtered sequence will have a spectrum $S(w)=|H(w)|^2$. More on filtering is discussed in Section~\ref{filtering}. \section{Polynomials and System Transfer Functions} \index{polynomials} \index{system transfer functions} Polynomials, matrix polynomials and transfer matrices are also defined and Scilab permits the definition and manipulation of these objects in a natural, symbolic fashion. Polynomials are easily created and manipulated. The {\tt poly}\index{poly@{\tt poly}} primitive in Scilab can be used to specify the coefficients of a polynomial or the roots of a polynomial. A very useful companion to the {\tt poly} primitive is the {\tt roots} primitive. The roots of a polynomial {\tt q} are given by : \begin{verbatim} -->a=roots(q); \end{verbatim} \index{function syntax!roots@{\tt roots}} The following examples should clarify the use of the {\tt poly} and {\tt roots} primitives. \verbatok{\Diary poly1.dia} \end{verbatim} Notice that the first polynomial {\tt q1} uses the {\tt 'roots'} default and, consequently, the polynomial takes the form $(s-1)(s-2)=2-3s+s^2$. The second polynomial {\tt q2} is defined by its coefficients given by the elements of the vector. Finally, the third polynomial {\tt q3} calculates the characteristic polynomial of the matrix {\tt a} which is by definition $\mbox{det}(sI-\mtt{a})$. Here the calculation of the {\tt roots} primitive yields the eigenvalues of the matrix {\tt a}. Scilab can manipulate polynomials in the same manner as other mathematical objects such as scalars, vectors, and matrices. That is, polynomials can be added, subtracted, multiplied, and divided by other polynomials. The following Scilab session illustrates operations between polynomials \verbatok{\Diary poly2.dia} \end{verbatim} Notice that in the above session we started by defining a basic polynomial element {\tt x} (which should not be confused with the character string {\tt 'x'} which represents the formal variable of the polynomial). Another point which is very important in what is to follow is that division of polynomials creates a rational polynomial which is represented by a list (see {\tt help list} and {\tt help type} in Scilab). A rational is represented by a list containing four elements. The first element takes the value {\tt 'r'} indicating that this list represents a rational polynomial. The second and third elements of the list are the numerator and denominator polynomials, respectively, of the rational. The final element of the list is either {\tt []} or a character string (More on this subject is addressed later in this chapter (see Section~\ref{sysrep}). The following dialogue illustrates the elements of a list representing a rational polynomial. \verbatok{\Diary poly3.dia} \end{verbatim} \subsection{Evaluation of Polynomials} \label{poleval} \index{polynomials!evaluation} A very important operation on polynomials is their evaluation at specific points. For example, perhaps it is desired to know the value the polynomial $x^2+3x-5$ takes at the point $x=17.2$. Evaluation of polynomials is accomplished using the primitive {\tt freq}. The syntax of {\tt freq} is as follows \begin{verbatim} -->pv=freq(num,den,v) \end{verbatim} The argument {\tt v} is a vector of values at which the evaluation is needed. For signal processing purposes, the evaluation of frequency response of filters and system transfer functions is a common use of {\tt freq}. For example, a discrete filter can be evaluated on the unit circle in the $z$-plane as follows \verbatok{\Diary poly4.dia} \end{verbatim} Here, {\tt h} is an FIR filter of length 9 with a triangular impulse response. The transfer function of the filter is obtained by forming a polynomial which represents the $z$-transform of the filter. This is followed by evaluating the polynomial at the points $\exp(2\pi in)$ for $n=0,1,\ldots,10$ which amounts to evaluating the $z$-transform on the unit circle at ten equally spaced points in the range of angles $[0,\pi]$. \subsection{Representation of Transfer Functions} \label{tfrep} \index{polynomials!representation of transfer functions} \label{sysrep} Signal processing makes use of rational polynomials to describe signal and system transfer functions. These transfer functions can represent continuous time signals or systems or discrete time signals or systems. Furthermore, discrete signals or systems can be related to continuous signals or systems by sampling. The function which processes a rational polynomial so that it can be represented as a transfer function is called {\tt syslin}: \begin{verbatim} -->sl=syslin(domain,num,den) \end{verbatim} \index{function syntax!syslin@{\tt syslin}} Another use for the function {\tt syslin} for state-space descriptions of linear systems is described in the following section. \section{State Space Representation} \label{ssrep} \index{state space representation} The classical state-space description of a continuous time linear system is : \begin{eqnarray} \dot{x}(t)&=&Ax(t)+Bu(t)\nonumber\\ y(t)&=&Cx(t)+Du(t)\nonumber\\ x(0)&=&x_0\nonumber \end{eqnarray} where $A$, $B$, $C$, and $D$ are matrices and $x_0$ is a vector and for a discrete time system takes the form \begin{eqnarray} x(n+1)&=&Ax(n)+Bu(n)\nonumber\\ y(n)&=&Cx(n)+Du(n)\nonumber\\ x(0)&=&x_0\nonumber \end{eqnarray} State-space descriptions of systems in Scilab use the {\tt syslin} function : \begin{verbatim} -->sl=syslin(domain,a,b,c [,d[,x0]]) \end{verbatim} \index{function syntax!syslin@{\tt syslin}} The returned value of {\tt sl} is a list where {\tt s=list('lss',a,b,c,d,x0,domain)}. The value of having a symbolic object which represents a state-space description of a system is that functions can be created which operate on the system. For example, one can combine two systems in parallel or in cascade, transform them from state-space descriptions into transfer function descriptions and vice versa, and obtain discretized versions of continuous time systems and vice versa. The topics and others are discussed in the ensuing sections. \section{Changing System Representation} \index{changing system representation} Sometimes linear systems are described by their transfer function and sometimes by their state equations. In the event where it is desirable to change the representation of a linear system there exists two Scilab functions which are available for this task. The first function {\tt tf2ss} converts systems described by a transfer function to a system described by state space representation. The second function {\tt ss2tf} works in the opposite sense. The syntax of {\tt tf2ss} is as follows \begin{verbatim} -->sl=tf2ss(h) \end{verbatim} \index{function syntax!tf2ss@{\tt tf2ss}} An important detail is that the transfer function {\tt h} must be of minimum phase. That is, the denominator polynomial must be of equal or higher order than that of the numerator polynomial. \begin{verbatim} -->h=ss2tf(sl) \end{verbatim} \index{function syntax!ss2tf@{\tt ss2tf}} The following example illustrates the use of these two functions. \verbatok{\Diary sysrep1.dia} \end{verbatim} Here the transfer function of a discrete IIR filter is created using the function {\tt iir} (see Section~\ref{desiir}). The fourth element of {\tt h1} is set using the function {\tt syslin} and then using {\tt tf2ss} the state-space representation is obtained in list form. \section{Interconnecting systems} \index{interconnecting systems} Linear systems created in the Scilab environment can be interconnected in cascade or in parallel. There are four possible ways to interconnect systems illustrated in Figure~\ref{f.intcon1}. % \begin{figure}[tb] \center{ \begin{picture}(200,300) \put(180,270){\framebox(0,0){\tt s1*s2}} \put(9,270){\circle{2}} \put(10,270){\vector(1,0){20}} \put(30,260){\framebox(30,20){$s_1$}} \put(60,270){\vector(1,0){30}} \put(90,260){\framebox(30,20){$s_2$}} \put(120,270){\vector(1,0){20}} \put(141,270){\circle{2}} \put(180,220){\framebox(0,0){\tt s1+s2}} \put(29,220){\circle{2}} \put(30,220){\line(1,0){20}} \put(50,220){\circle*{2}} \put(50,200){\line(0,1){40}} \put(50,200){\vector(1,0){20}} \put(50,240){\vector(1,0){20}} \put(70,190){\framebox(30,20){$s_2$}} \put(70,230){\framebox(30,20){$s_1$}} \put(100,200){\line(1,0){20}} \put(100,240){\line(1,0){20}} \put(120,240){\vector(0,-1){15}} \put(120,200){\vector(0,1){15}} \put(120,220){\circle{10}} \put(120,220){\framebox(0,0){$+$}} \put(125,220){\vector(1,0){15}} \put(141,220){\circle{2}} \put(180,140){\framebox(0,0){\tt [s1,s2]}} \put(49,160){\circle{2}} \put(49,120){\circle{2}} \put(50,120){\vector(1,0){20}} \put(50,160){\vector(1,0){20}} \put(70,110){\framebox(30,20){$s_2$}} \put(70,150){\framebox(30,20){$s_1$}} \put(100,120){\line(1,0){20}} \put(100,160){\line(1,0){20}} \put(120,160){\vector(0,-1){15}} \put(120,120){\vector(0,1){15}} \put(120,140){\circle{10}} \put(120,140){\framebox(0,0){$+$}} \put(125,140){\vector(1,0){15}} \put(141,140){\circle{2}} \put(180,50){\framebox(0,0){\tt [s1;s2]}} \put(49,50){\circle{2}} \put(50,50){\line(1,0){20}} \put(70,50){\circle*{2}} \put(70,30){\line(0,1){40}} \put(70,30){\vector(1,0){20}} \put(70,70){\vector(1,0){20}} \put(90,20){\framebox(30,20){$s_2$}} \put(90,60){\framebox(30,20){$s_1$}} \put(120,30){\vector(1,0){20}} \put(120,70){\vector(1,0){20}} \put(141,30){\circle{2}} \put(141,70){\circle{2}} \end{picture} } \caption{Block Diagrams of System Interconnections} \label{f.intcon1} \end{figure} % In the figure the symbols $s_1$ and $s_2$ represent two linear systems which could be represented in Scilab by transfer function or state-space representations. For each of the four block diagrams in Figure~\ref{f.intcon1} the Scilab command which makes the illustrated interconnection is shown to the left of the diagram in typewriter-face font format. \section{Discretizing Continuous Systems} \index{discretization of continuous systems} A continuous-time linear system represented in Scilab by its state-space or transfer function description can be converted into a discrete-time state-space or transfer function representation by using the function {\tt dscr}. Consider for example an input-output mapping which is given in state space form as: % \begin{equation} (C) \left\{ \begin{array}{ccc} \dot{x}(t) & = & A x(t) + B u(t) \\ y(t) & = & C x(t) + D u(t) \end{array} \right. \label{cont} \end{equation} % From the variation of constants formula the value of the state $x(t)$ can be calculated at any time $t$ as % \begin{equation} x(t)=e^{At}x(0)+\int_{0}^{t}e^{A(t-\sigma)}Bu(\sigma)d\sigma \label{voc} \end{equation} % Let $h$ be a time step and consider an input $u$ which is constant in intervals of length $h$. Then associated with (\ref{cont}) is the following discrete time model obtained by using the variation of constants formula in (\ref{voc}), % \begin{equation} (D) \left\{ \begin{array}{ccc} x(n+1) & = & A_{h} x(n) + B_{h} u(n)\nonumber \\ y(n) & = & C_{h} x(n) + D_{h} u(n) \end{array} \right. \label{disc} \end{equation} % where $$A_{h}=\exp(Ah)$$ $$B_{h}=\int_{0}^{h}{e^{A(h-\sigma)}Bd\sigma}$$ $$C_{h}=C$$ $$D_{h}=D$$ Since the computation of a matrix exponent can be calculated using the Scilab primitive {\tt exp}, it is straightforward to implement these formulas, although the numerical calculations needed to compute $\exp(A_h)$ are rather involved (\cite{VanLoan}). If we take $$G=\left\lgroup\matrix{A&B\cr0&0\cr}\right\rgroup$$ where the dimensions of the zero matrices are chosen so that $G$ is square then we obtain $$\exp(Gh)=\left\lgroup\matrix{A_h&B_h \cr0&I\cr}\right\rgroup$$ When $A$ is nonsingular we also have that $$B_h=A^{-1}(A_h-I)B.$$ This is exactly what the function {\tt dscr} does to discretize a continuous-time linear system in state-space form. The function {\tt dscr} can operate on system matrices, linear system descriptions in state-space form, and linear system descriptions in transfer function form. The syntax using system matrices is as follows \begin{verbatim} -->[f,g[,r]]=dscr(syslin('c',a,b,[],[]),dt [,m]) \end{verbatim} \index{function syntax!dscr@{\tt dscr}} where {\tt a} and {\tt b} are the two matrices associated to the continuous-time state-space description % \begin{equation} \dot{x}(t)=Ax(t)+Bu(t)\nonumber \end{equation} % and {\tt f} and {\tt g} are the resulting matrices for a discrete time system % \begin{equation} x(n+1)=Fx(n)+Gu(n)\nonumber \end{equation} % where the sampling period is {\tt dt}. In the case where the fourth argument {\tt m} is given, the continuous time system is assumed to have a stochastic input so that now the continuous-time equation is % \begin{equation} \dot{x}(t)=Ax(t)+Bu(t)+w(t)\nonumber \end{equation} % where $w(t)$ is a white, zero-mean, Gaussian random process of covariance {\tt m} and now the resulting discrete-time equation is % \begin{equation} x(n+1)=Fx(n)+Gu(n)+q(n)\nonumber \end{equation} % where $q(n)$ is a white, zero-mean, Gaussian random sequence of covariance {\tt r}. The {\tt dscr} function syntax when the argument is a linear system in state-space form is \begin{verbatim} -->[sld[,r]]=dscr(sl,dt[,m]) \end{verbatim} where {\tt sl} and {\tt sld} are lists representing continuous and discrete linear systems representations, respectively. Here {\tt m} and {\tt r} are the same as for the first function syntax. In the case where the function argument is a linear system in transfer function form the syntax takes the form \begin{verbatim} -->[hd]=dscr(h,dt) \end{verbatim} where now {\tt h} and {\tt hd} are transfer function descriptions of the continuous and discrete systems, respectively. The transfer function syntax does not allow the representation of a stochastic system. As an example of the use of {\tt dscr} consider the following Scilab session. \verbatok{\Diary dscr1.dia} \end{verbatim} \section{Filtering of Signals} \label{filtering} Filtering of signals by linear systems (or computing the time response of a system) is done by the function {\tt flts} which has two formats . The first format calculates the filter output by recursion and the second format calculates the filter output by transform. The function syntaxes are as follows. The syntax of {\tt flts} is \begin{verbatim} -->[y[,x]]=flts(u,sl[,x0]) \end{verbatim} \index{function syntax!flts@{\tt flts}!for state-space} for the case of a linear system represented by its state-space description (see Section~\ref{ssrep}) and \begin{verbatim} -->y=flts(u,h[,past]) \end{verbatim} \index{function syntax!flts@{\tt flts}!for transfer function} for a linear system represented by its transfer function. In general the second format is much faster than the first format. However, the first format also yields the evolution of the state. An example of the use of {\tt flts} using the second format is illustrated below. \verbatok{\Diary Sflts2.dia} \end{verbatim} Notice that in the above example that a signal consisting of the sum of two sinusoids of different frequencies is filtered by a low-pass filter. The cut-off frequency of the filter is such that after filtering only one of the two sinusoids remains. Figure~\ref{f.flts1} illustrates the original sum of sinusoids and Figure~\ref{f.flts2} illustrates the filtered signal. % \input{\Figdir flts1.tex} \dessin{{\tt exec('flts1.code')} Sum of Two Sinusoids}{f.flts1} % % \input{\Figdir flts2.tex} \dessin{{\tt exec('flts2.code')} Filtered Signal}{f.flts2} \section{Plotting Signals} \index{plotting} Here we describe some of the features of the simplest plotting command. A more complete description of the graphics features are given in the on-line help. Here we present several examples to illustrate how to construct some types of plots. To illustrate how an impulse response of an FIR filter could be plotted we present the following Scilab session. \verbatok{\Diary Splot1.dia} \end{verbatim} \index{plotting!impulse response} Here a band-pass filter with cut-off frequencies of .2 and .25 is constructed using a Hamming window. The filter length is 55. More on how to make FIR filters can be found in Chapter~\ref{fir}. % \input{\Figdir plot1.tex} \dessin{{\tt exec('plot1.code')} Plot of Filter Impulse Response}{f.plot1} % The resulting plot is shown in Figure~\ref{f.plot1}. The frequency response of signals and systems requires evaluating the $s$-transform on the $j\omega$-axis or the $z$-transform on the unit circle. An example of evaluating the magnitude of the frequency response of a continuous-time system is as follows. \verbatok{\Diary Splot2.dia} \end{verbatim} \index{plotting!continuous magnitude} Here we make an analog low-pass filter using the functions {\tt analpf} (see Chapter~\ref{iir} for more details). The filter is a type I Chebyshev of order 4 where the cut-off frequency is 5 Hertz. The primitive {\tt freq} (see Section~\ref{poleval}) evaluates the transfer function {\tt hs} at the values of {\tt fr} on the $j\omega$-axis. The result is shown in Figure~\ref{f.plot2} % \input{\Figdir plot2.tex} \dessin{{\tt exec('plot2.code')} Plot of Continuous Filter Magnitude Response}{f.plot2} % A similar type of procedure can be effected to plot the magnitude response of discrete filters where the evaluation of the transfer function is done on the unit circle in the $z$-plane by using the function {\tt frmag}. \begin{verbatim} -->[xm,fr]=frmag(num[,den],npts) \end{verbatim} \index{function syntax!frmag@{\tt frmag}} The returned arguments are {\tt xm}, the magnitude response at the values in {\tt fr}, which contains the normalized discrete frequency values in the range $[0,0.5]$. \verbatok{\Diary Splot3.dia} \end{verbatim} \index{plotting!discrete magnitude} Here an FIR band-pass filter is created using the function {\tt eqfir} (see Chapter~\ref{fir}). % \input{\Figdir plot3.tex} \dessin{{\tt exec('plot3.code')} Plot of Discrete Filter Magnitude Response} {f.plot3} % Other specific plotting functions are {\tt bode} for the Bode plot of rational system transfer functions (see Section~\ref{bode}), {\tt group} for the group delay (see Section~\ref{group}) and {\tt plzr} for the poles-zeros plot. \verbatok{\Diary Splot4.dia} \end{verbatim} \index{plotting!poles and zeros} Here a fourth order, low-pass, IIR filter is created using the function {\tt iir} (see Section~\ref{desiir}). The resulting pole-zero plot is illustrated in Figure~\ref{f.plot4} % \input{\Figdir plot4.tex} \dessin{{\tt exec('plot4.code')} Plot of Poles and Zeros of IIR Filter}{f.plot4} % \section{Development of Signal Processing Tools} Of course any user can write its own functions like those illustrated in the previous sections. The simplest way is to write a file with a special format . This file is executed with two Scilab primitives {\tt getf} and {\tt exec}. The complete description of such functionalities is given in the reference manual and the on-line help. These functionalities correspond to the button {\tt File Operations}.