Development of a Simulation Environment for Signal Processing Algorithms

The development of a simulation environment for the study of signal processing algorithms is discussed. The project was divided into two parts. The first and the major part was the development of MATLAB, a matrix manipulation program, into a suitable simulation environment for studies in signal processing. The second part used the software tools developed by the first in performing a simulation study involving the use of state space models in parameter estimation. Developing software tools in MATLAB included expanding the original inbuilt set of tools to enable the users to write their own functions effectively. This meant creating a set of sophisticated tools and adding these to the source of the MATLAB program. The method for adding functions to the MATLAB library for use by the whole programming community has been detailed. A large set of signal processing functions have been added to the MATLAB library, thus increasing the power of MATLAB as a signal processing simulation environment. The problem of approximate realization which constitutes the second part of the project, has been dealt with at the simulation stage. Starting out with noisy data, the aim was to develop a state space model which approximated this data. Different ways of arriving at a solution have been considered. The results of the simulation studies are detailed.


The problem in brief
Problem solving is a complex issue. To reach a solution, one can approach it in many ways. Our interest was in the solution of signal processing problems. Analysis is generally the first step in approaching a problem.
However, another way to approach a problem is by a simulation study.
Simulation studies can often tell us what needs to be analyzed. Simulation environments for the study of signal processing algorithms are available in plenty. However, most of the environments have their own limitations.
Our intention was to develop an environment which would enhance the ease of doing simulation studies with signal processing algorithms that are described in the language of linear algebra. Such algorithms deal with arrays and matrices. Writing programs in FORTRAN or any other high level language involves accessing each element of the array or matrix, and this could be quite cumbersome. We were thus in search for a language which would enable us to eliminate this problem of accessing each memory cell and thereby make signal processing problem solving an easier task.
MATLAB, a matrix manipulation program was found to have the desired qualities. A description of the MATLAB language is given later on in this chapter. 1 The project was divided into two parts. The first was to develop software tools which wo11 ld automatically perform signal processing functions.
The second was to use these tools in a simulation study of some signal processing algorithms. The software tools included building a set of signal processing functions which were to be incorporated as primitive MATLAB functions. The simulation studies would involve performance comparison between different state variable algorithms for parameter estimation and signal modeling.

Programming language classifications
Computer programming languages can be broadly classified into two kinds [1]. Statem'ent oriented languages and languages oriented towards math- Of the two, imperative languages are more efficient in terms of execution time because they reflect the structure and operations of the machine. However, they require the programmer to pay attention to machine level details.
Programming using imperative languages is based on naming elementary memory cells, assignments to these cells and repetition of elementary operations. The functional programming style does not depend on these three actions, largely due to the simple and uniform data objects (like arrays, lists etc.) which allow design of data structures without concern for elementary memory cells. Functional programming thus appears to be at a higher level than imperative programming.

MATLAB -its superiority over others
MATLAB is a powerful computer language for matrix manipulations and calculations [2,3]. Of the above mentioned programming language categories, MATLAB has the charecteristics of applicative languages. It has only one data type; a rectangular matrix with complex elements. Scalars and vectors are special cases of this rectangular matrix. MATLAB is a program which can be worked in both the interactive and batch modes.
Among the attractions of interactive working are immediate response and ease of correction of errors. The language is very easy to use and learn since a statement in MATLAB is almost identical to standard mathematical notation. A MATLAB program is easy to follow and one can trace out the path of the program in the first run by simply going through it. One of our objectives in chosing a language has been to consider the speed of writing a problem. The development time of writing a MATLAB program is small compared to writing it in any other high level language like FORTRAN.
There are many reasons for this. A single data type in MATLAB eliminates complex data conversion statements and also eliminates equating elements of different data types as in other languages. As the language is based on the concept of arrays and matrices and processes these as entities, it eliminates the need for many of the loops necessary in other languages.
The language also eliminates the need to dimension matrices. It is easy to append a matrix to another or to cut a larger matrix to a shorter one with any desired dimensions. The programs are thus easier to debug.
MATLAB has a large number of built in functions. An important feature of MATLAB is its built in capability for growth. It is easy to add functions written by the user to his or her own library or to the MATLAB 3 library for use by all users. MATLAB however is an interpreted language and thus is not very efficient. To solve a program with many loops takes a Jong time. However, as explained in chapter 4, this problem can be reduced, if not eliminated, with the help of Extended functions.
MATLAB is thus a valuable tool for most computer users even though it appears to be oriented towards the novice user. It is just like a "super desk calculator" although one with a vast memory and a rich and powerful vocabulary. It is intended to provide its users with the means to quickly formulate problems for computer processing and for obtaining answers in a minimum of time. MATLAB can thus be thought of as a "very high level" language. Chapter 5 deals with an introduction to approximate realization and model reduction using the state space approach, and also the simulation example involving the use of state space models in paramter estimation. 4

Contributions of this thesis
• Gives a detailed discussion on writing and implementing extended functions for MATLAB. Extended functions were developed at URI to enhance the power of MATLAB.
• Documents 27 new signal processing functions which have been installed in the MATLAB language along with an on line HELP facility.
• Documents the procedures for installing new functions and HELP files into MATLAB.
• Demonstrates the usefulness of these tools through a simulation study of state space methods for parameter estimation. Simulation results are giv: n which compare the performance of an existing method with two new methods.
• Documents the usage of the new signal processing functions added to the library through a "users manual" with help files and test cases.
This document is available in the department.

Usage
In the MATLAB environment, the USER (or US Rx) function can be used in the following way. In the \1 ATLAB environment: 4, 5, 6 7, 8, 9 > <> EG=USER(D,S,T) The value of EG would be 11.
We are adding the parameters E and F to A(l,1) in this simple example \ and then returning the value of A(l,1) back to the MATLAB environment.
Note that initially the value of both 0 and P was 3, but in the subroutine, we change it to 1 at the end. This means that we are returning the matrix A with new dimensions lxl to the MATLAB environment. Also note that in the MATLAB environment, the dimensions of the matrix D is not specified in the call to the USER function. This is because the dimensions are picked up by the system inherently and there is no need for the user to specify them. This makes the Language very useful and powerful.
MATLAB allows the input matrix to the USER function to be redimensioned, so that the output matrix can be of a different dimension than the input matrix. This holds true with one constraint. The output matrix can only be smaller than the input matrix and not larger. ff we have a matrix A of dimension lOxlO which is the input matrix to a USER function, the output matrix could have a dimension anywhere from lxl to lOxlO, but not lOxll.
In order to remove the above mentioned constraints of the USER functions, the Extended functions were introduced into the MATLAB environment. 8

What are Extended functions ?
To enhance the usability of the USER functions, the Extended functions were incorporated into the MATLAB environment. The Extended functions are depicted by "EXTx" , where "x" is any english letter. The main differences between the USER and EXTx functions are : • We can input as many matrices into the FORTRAN environment as we wish using the EXTx function, as against only one by the USER function.
• We can input as many scalar variables into the FORTRAN environment using the EXTx function, as against only two by the USER function. TAe scalar variables have to be dimensioned as lxl matrices in the EXTx function. So in effect, the scalar variables are special cases of matrices.
• A complex scalar, which is a special case of a complex matrix, can also be inputted. This would have to be declared as a COMPLEX*16 variable and will have to be dimensioned as lxl. A real scalar can be declared as DOUBLE PRECISION.
• The input matrices and scalars can now contain complex numbers, as against only real numbers in the USER function.
• The output matrix from the EXTx function can be dimensioned smaller or larger than the input matrix, as against a smaller or equal dimensioned output matrix with the USER function. However, there is an upper limit to the dimensioning of the output matrix.

Usage
The EXTx function is defined in the following way in the MATLAB environment.
where FUN2, FUN3, FUN4 are input matrices and A is an input scalar.
F NI is the output matrix.
The value of FUNl on return would be We see that the value of FUN2 is unchanged in the MATLAB environment.
Another example to illustrate the power of the extended function is as follows.
We have an input matrix A which is of dimension 2x2 and another matrix B which is also of dimension 2x2. We would like to place both these matrices side by side and multiply the resultant matrix with a scalar and then return this 2x4 matrix. This is how we would do it.
In the FORTRAN environment, we would have to write the following sub- c Note that the matrix D has not been dimensioned and hence c it could be considered as a scalar. c c Test to be sure that the result array will be smaller than c the maximum allowed number of elements. The number of C elements is passed in RROW. This subroutine has to be in a file extd.f and this file has to be compiled and then linked with MATLAB. A sample "make file" to achieve this is given at the end of this section in figure 2.2.1.
In the MATLAB environment, we have the following.
Another example which will clearly illustrate the usage of the extended function is the formation of a Hankel matrix from a data sequence.
The FORTRAN subroutine is as follows. We have created a 7x8 hankel matrix from the given data sequence.
A sample "make file" to compile and link files to MATLAB is given here.  Table 2.3. Note that as the values of Mand N increase, the "function" call time remains almost constant, whereas the "MATLAB" call time increases exponentially.
One goal of this thesis was to write a set of extended functions for many common signal processing tasks. These extended functions were given descriptive command names and incorporated into the MATLAB language itself. Thus they are available to any MATLAB user in either interactive or batch mode. In addition, on-line HELP files for these new functions were provided.
The documentation on extended functions given in this chapter gives all the information required for a user to write his or her own special purpose extended functions. We see that an extended function for a one line MAT-LAB code consists of quite a few lines of FORTRAN code. However, if the speed is desired, and your MATLAB program, which consists of this line is used often, it may be worth while to make use of the extended function.
The only disadvantage in using the Extended function for DO Loops could be while using matrices for computations inside the loops. The main purpose of using MATLAB, where you do not need to specify the dimensions of matrices during computations, is defeated. In the FORTRAN I M N function MATLAB I on. This chapter also deals with problems of filtering. The lowpass and bandpass designs used here are quite standard [8] and we have explained how these designs were used to make up the MATLAB function library.
The bandpass and lowpass filter designs use the filtering routines which filter data using the ARMA parameters. The function for filtering using the ARMA parameters has been included in the MATLAB library, but has not been detailed here. Its definition is given in chapter 4.

FFT
The use of Fast Fourier Transform algorithms in Digital Signal Processing became widespread as a result of the well known paper by Cooley and Tukey [16]. However, most of the applications were oriented towards using the value of N, the data sequence length, as a power of 2, as this decomposition lead to a highly efficient computation when using Cooley and Tukey's algorithm. For a majority of applications, this restricted choice of N was adequate. The need to have the value of N which was not a power of 2 was encountered by Singleton [13] in spectral analysis of speech and economic time series data. A mixed-radix FFT algorithm was developed by a lot of researchers but the one by Singleton [13] surpassed most others in terms of ease of use, fiexibilitj and efficiency. We have thus considered the mixedradix FFT algorithm by Singleton and have included it in the MATLAB library. Singleton FFT has been used very widely in a lot of applications due to its efficiency and flexibility in chasing different transform lengths [14].
Singleton's algorithm for computing the Fast Fourier Transform is based on the method proposed by Cooley and Tukey, but is a mixed radix algorithm. In his paper, Singleton presents an improved method of computing a transform step corresponding to an odd factor p of N. For details on how his algorithm works, the reader is referred to [13].

FFT related topics
is the linear convolution of x1(n) and x2(n). xs(n) can have at most 2N -1 non zero points. The above convolution of 2 data sequences can be done directly using equation 3.1 or by using discrete fourier transforms to arrive at x 3 (n). Using discrete fourier transforms, xs(n) would be obtained by first finding the transforms of X1 ( n) and x 2 ( n) for 2N -1 points and then finding the inverse transform of their product. The following equations describe this method.
Due to the presence of highly efficient algorithms available for computing the discrete fourier transform of finite duration sequences, it has been found that using the transform method instead of the direct method would lead to savings in time as the number of multiplications would decrease by using the transform method. For short sequences, it might still be better to use the direct method of equation 3.1.
In building the convolution function for the MATLAB library, we have used only the transform method. The overhead by using the transform method for short sequences as compared to the direct method is very small and hence the direct method has not been used. For extremely long sequences, the overlap add and the overlap save methods can be used [B,10,12]. However, we are dealing with a set of finite sequences. Real time data acquisition is not being considered, thus eliminating the use of very long sequences in our applications. MATLAB is an interactive program and the data has to be present to be acted upon. We have thus not used the overlap add or the overlap save methods in our programs.
The MATLAB function CONV has been written to perform convolution.
It is documented in chapter 4.

Correhltion
Calculating the autocorrelation estimate is similar to that of calculating the convolution of two sequences. The autocorrelation of the sequence x(n) corresponds to the circular convolution of x( n) with x( -n). We can thus use the FFT to efficiently compute the autocorrelation estimate The MATLAB function CORR has been written to perform correlation.
It is documented in chapter 4.

Filters
In my approach to creating functions for Digital Filters, I have essentially used the same designs for the Butterworth lowpass filter and the band-pass filter as described by Oppenheim and Schafer [8]. However, a few assumptions have been made and they are detailed ahead.

3.3.l Lowpass Design
The squared magnitude function for an analog Butterworth filter is of the form I H ( ·o)J2 - where De is the analog cutoff frequency of the filter. We know that as the value of N, the order of the filter, increases, the charecteristics of the filter become sharper. In the S plane, the equation 3.3 can be written as ' 1 Ha ( Filtering the data has been done by designing 2nd order sections from the poles of the analog butterworth filter and then using these sections in a cascaded manner. We see that if N, the order of the filter is odd, then the filter has just one first order section given by (3.4) 26 and m == (N -1)/2 second order sections given by 02 c • k = 1,2, ... m (3.5) s2 + 20c cos (7rk/N)s + O~' If N is even, then we have m = N /2 second order sections given by The user provides a normalized digital frequency a which is to be prewarped to Oc for use in the equations 3.4, 3.5 and 3.6. The normalized digital frequency a should have a value between 0 and 1 with both limits excluded.
The prewarped cutoff frequency Oc is given by where T is the sampling period. The solution for a from the above equation to arrive at a digital Butterworth lowpass filter design. Each second order section is then used to filter the input data. Filtering is done using a function which uses the ARMA parameters of each second order section generated by the above equations. Note that when equations 3.7 and 3.8 are substituted in equations 3.4, 3.5 and 3.6, T, the sampling period disappears. We thus see that a change in the value of T has no effect on the design and thus we do not need to input the value of T.
A HELP documentation of the subroutine LPF, which designs the lowpass filter as described above, is given in the chapter 4. 27 3 2 Bandpass Design 3. .
The traditional approach to the design of bandpass filters is to first design a frequency normalized prototype lowpass filter and then use an algebraic transformation, usually the spectral transformation, to derive the desired bandpass filter. We have used this approach in our design of the bandpass filter . Our methodology includes designing a lowpass butterworth filter as described in Section 3.3.1 with an normalized digital cutoff frequency of 0.25 and then using the spectral transformation to arrive at the appropriate design.
One thing to note while using the spectral transformations is that an Ntl" order lowpass filter produces a 2Nth order bandpass filter. Thus if the ' user specifies an Nth order bandpass filter to be designed, a N/2th order lowpass filter is designed and then a corresponding Nth order bandpass filter is constructed. In the event that N is odd, an Ni 1 th order lowpass filter is designed and correspondingly, an (N + l)th order bandpass filter is constructed. This higher value of N has been chosen instead of a lower value as the filter charecteristics of this higher order filter would be sharper than the desired filter, with no loss to the user.
The MATLAB function BAND has been written to design a bandpass filter. It is documented in chapter 4.

28
Chapter 4 Building the Matlab library

Summary
To change the help file, do the following in the given order.  If a change is made to both the 'help' and the 'helper.f' files, then follow the respective procedures in the right order.

Adding Functions
A detailed description of adding functions to the MATLAB library is given here. The procedure detailed here requires editing quite a few files in the MATLAB source. Opening files and looking through it in the same order as I have detailed would clarify the procedure and lead to little or no confusion.
A summary at the end of this section describes the procedure in brief.
FORTRAN is a scientific language used mostly for mathematical computations. It is not built to handle character strings very well. The whole of MATLAB is written in FORTRAN. MATLAB handles strings by first converting them to numbers and then following the necessary procedures.
MATLAB has a character set defined within it which has a one to one correspondence between a character or a special symbol and a number. The character set of MATLAB is given in figure 4.2.
Thus the number corresponding to the alphabet 'A' is 10 and 'B' is 11 and On The string 'ABS' is given by 10,11,28. so   • IMPZ(P,Z,G,N) -Finds the impulse response of a system using Pole Zero variables. P and Z are the poles and zeroes of the system. The gain of the system is specified in G. N is the length of the desired output data sequence (impulse response).

Existing MATLAB functions
• IMTF(B,A,N) -Finds the impulse response of a system using Transfer Functions coefficients. The numerator coefficients are specified in "B" and the denominator coefficients are specified in "A". N is the length of the desired output data sequence (impulse response). The leading coefficients of B and A are inputted, even if they are 1.
• NOSR(S,SNR) -Adds real noise to the Data Sequence S with the desired SNR. The data sequence S can be either real or complex.
The SNR is defined as 1 The following functions are based on programs obtained from [9]  The definition of SNR is the given above.
• SOTA( A) -Sorts an array A of numbers in ascending order. The array A can be real or complex. Complex arrays are sorted in magnitude.
• SOTD(A) -Sorts an array A of numbers in descending order. The array A can be real or complex.Complex arrays are sorted in magnitude.
• ATN(A) -Finds the four quadrant arctangent of a vector or scalar A.
• TSER(A,B,N) -Generates a real AR, MA, or ARMA time series.

'
A is an array of AR filter parameters. B is an array of MA filter parameters. The variance of real excitation noise is assumed to be 1.0.
• TSEC(A,B,N) -Generates a complex AR, MA, or ARMA time series.
A is an array of AR filter parameters. B is an array of MA filter parameters. The variance of real and imaginary parts of the excitation noise is assumed to be 0.5.
• LEVN (R) -Implements the Levinsons recursion. R is an array of autocorrelation samples. On return, we get the solution of the Levinsons recursion.
• LERR(R) -Implements the Levinsons recursion. R is an array of autocorrelation samples. The prediction error powers are returned by this function.
• FFT(A,IS) -Finds the fast fourier transform of the data sequence in A. The value of IS indicates if the FFT or the inverse FFT has to be found. If IS is 1.0, the FFT is found and if the value of IS is -1.0, the inverse is found. The definition of the FFT is given in chapter 3.
• FACT(N) -Determines the next number closest to N which has a prime factor of 23 or less. This is to be used in conjunction with the FFT function as it uses the mixed radix transform for sequences whose length has prime factors less than or equal to 23.
• LPF(X,N,W) -Designs a lowpass Butterworth filter of order N and a normalized digital cutoff frequency W. The cutoff frequency has to be in the range 0 to 1.0, with both extremes excluded. X is the input real or complex data sequence. The definition of the filter is given in chapter 3.
• BAND(X,N,Wl,W2) -Designs a bandpass filter of order N and normalized digi~al frequencies of Wl and W2. Both Wl and W2 have to be in the range 0 to 1, with both extremes excluded. X is the input real or complex data sequence which is to be filtered. The definition of this function is given in chapter 3.
• CONV(X,Y) -Convolve a sequence X with a sequence Y.

Introduction
The problems of model reduction and approximate realization for discretetime, deterministic linear systems may be embedded into one another due to their similar nature. An introduction to solving these problems has been made using the state space models. It has been shown [11] that state space representation has certain advantages for approaching the problems of approximate realization and model reduction and we will be dealing with this approach here. Unlike the traditional scalar valued data dealt with by the autoregressive(AR), moving average(MA) or autoregressive moving average (ARMA) models, the state space approach deals with several variables simultaneously as vector valued variables. The AR, MA or ARMA models can deal with vector valued data, but the problems of model reduction and approximate realization are difficult to solve when compared with the state space approach. Theoretically, the Markovian or state space models are equivalent to the traditional ARMA type models because models in one representation can be transformed into another. However, there are some grounds in thich these two models are not equivalent, like numerical stability, statistical properties of parameter estimates, ease of handling for vector valued or nonstationary time series represented in these two al-44 ternative modes. We will not go into detail about the differences or the likes of these two modes here. The interested reader is referred to [11].
This reference deals with state space models for random processes. This chapter deals with the simpler problem of finding state space models for deterministic signals observed in noise.

The problem of approximate realization
The problem of approximate realization may be stated as follows. Given a data sequence, find a system whose impulse response approximates this data sequence in some sense. To explain this, let us construct a state space model of pth ordw. We will assume that the model is stable, i.e. the roots of the system are within the unit circle. The impulse response of the system can be described with the help of the state matrices A, b and c as where hk is the impulse response for k = 1, 2, ... The dynamic behaviour of the system is given by the following equations. (5.2) where A is p x p, b is p x 1 and c is 1 x p.
The Hankel matrix of a system is formed from the impulse response as follows. where em-1 is the first m -1 rows of e, and e~-1 is them -1 rows of 8 starting from the second row. Note that both these matrices are of rank p. We therefore have A -eL er m-1 m-1 (5.11) where e~_ 1 is the left inverse of em-1· We see that e~-l em-1 = Ip,p if 0m-l is of rank p. Thus the realization which corresponds to a factorization of II is simply (A, b, c) given by equations 5. 6 5.7 and 5.11. From equation 5.10, we see that both em and e i span the same p dimensional vector space. We have here considered only single-input, single-output systems. For multivariate systems, the realization forms are similar.
A factorization which is of prime importance is the singular value decomposition (SVD) of H. On taking the SVD of H, its factorization is of the following form.
where E is a matrix of singular values of H.  (5.13) where E 1 is an error matrix.
A least squares solution to this is proposed. If we can find E 1 to be minimum, then we can get a good approximation to A. Thus, choosing A to minimize (5.14) is clearly a reasonable thing to do. This choice of A was first proposed by Kung [15] and has been used by many others [4,5,6,7].
However, another solution could be found by writing the equation 5.13 in the following way. (5.15) where E2 denotes another error matrix. Thus, chosing B to minimize (5.16) would give us another good solution. Here A= B- 1 • We note that the above solutions project columns of e i onto the space spanned bye, or vice versa. This assumes that e or < §)i space is correct.
But this is not correct as the noise affects both e and e i. So we propose total least squares solution. (5 .17) such that is minimized. The solution to this minimization problem is given in the next section.
Retracing our steps a bit, how would we know that the noisy Hankel matrix is to be approximated with a p 0 • order Hankel matrix and not a (p + l)th order matrix. (Note that this pth order matrix is not really Hankel.) On taking the SYD of the noisy Hankel matrix H, we know that the 2: matrix contains the singular values of the system. Thus 0 0 0 0 0 an with ai > a2 > · · · > an > 0, for Hm,n with m > n. Also we note that Hm,n is of rank n.
In theory, the rank of a Hankel matrix is equal to the number of its non zero singular values. If we have a noisy Hankel matrix, all its singular values will be positive. In order to approximate this Hankel matrix to a smaller rank matrix, we have to determine if ai ~ ai+ 1 for all i = 1, 2, ... n.
Methods for determining this are well established [11]. Thus if ap ~ ap+b then the rank p matrix Hm,n will nearly be Hankel. Thus we get a nearly Hankel matrix of rank p which is a good approximation to H, and so the above procedure will yield a pth order model which is nearly optimal with respect to the spectral norm of the Hankel matrix H. In addition, we thus 50 see that a disparity in the magnitude of the singular values (up ~ O"p+i) gives an indication of the model order to which the given system can be reduced without much change in the input -output behaviour.

Simulation study
The problem of approximate realization has been dealt with in the simulation study. The theoretical concepts have been detailed in the previous section. Given a data sequence, the aim is to find a system whose impulse response approximates this data sequence in some sense. The parameters of interest are the angles of the poles of the system. The given data sequence is noisy. The ~ethods given in the last section for solving the A matrix of a state space model were used. The eigenvalues of A are the poles of the system.

METHOD A
This method uses equation 5.14 to get the desired approximate model.
The solution to this equation can be written as "'L "'j A= 8m-18m-1 (5.18) em-1 and e!n-1 are !ormed as mentioned in the previous section. e~-1 is the left inverse of 8m-1· The smallest p singular values are set to zero, and only the first p columns and rows of U 2 p and V 2~ are retained resulting in the optimal rank p approximating matrix UEVT. The solution then proceeds as follows. yT This would give us the following equation to solve.
and we could get the desired model as In the MATLAB notation, equation 5.20 could be written as \ 5.3.1 Simulation procedure (5.20) The procedure to find an approximate model from the given exact model is as follows. The exact model is given by a set of poles and zeroes which charecterizes the desired system.
• Generate the impulse response of the system from the poles and zeroes using the function IMTF.
• Add noise with a desired SNR to the data sequence inside the loop using function NOSC, which adds complex noise to the desired data sequence.
• Form a Hankel matrix from the noisy sequence inside the loop using the function HANK.
• Use methods A,B and C to find out the poles from the generated approximate model.
• Compare the bias, variance and mean square error of the angles of the generated poles.

Simulation results
Noise has been added to the data sequence in two ways. Adding noise to the data sequence and then forming the Hankel matrix has been one of the ways. We call this the 'correlated noise matrix' method since the additive noise matrix ~as Hankel structure. The other way has been to form a matrix of uncorrelated noise and adding this to the Hankel matrix of the data sequence to form the noisy matrix, which is now no longer Hankel. We call this the 'uncorrelated noise matrix' method. This method generates data similar to that in an array processing problem where the noise present at each sensor is uncorrelated with the noise at other sensors, and the noise at a single sensor is uncorrelated at different instants of time.
The system used for this simulation had poles at 0.6108 and 0. 7854 radians, both with radius 0.99. The length of each data sequence generated was 20 points. The number of trials in this study was 100.
'Bias' is the absolute difference between the sample mean and the actual value. 'Variance' is defined as the sum square of the difference between the estimated samples and the actual data over the number of iterations. 'Mean Square Error' (MSE) is the sum of the bias square and the variance.
The tables 5. 1, 5.2, 5.3, 5.4, 5.5, and 5.6 summarize the results obtained from methods A, Band C. The SNR is the Signal to Noise ratio in decibels.
Var is the variance of the observed data and MSE is the Mean Square Error.
We have only accumulated statistics of the angles 0 1 and 0 2 of the poles of the system.
From the study, the preliminary conclusion is that none of the methods has an edge over the others in solving this problem. The reason that the total least squares solution does not give a better result than the other two is because the noise in the matrices Gm-l and 0~_ 1 is correlated. More analysis needs to be done to arrive at the right solution. We have shown that simulations can be done using the tools available to us in the MATLAB environment.

Simulation Efficiency
The above simulation study has been done using some of the new Digital Signal Processing functions which have been added to the MATLAB library.
A comparision between writing pure MATLAB code without any of the new functions for ~oing the above simulation study and code with the new functions is made here. The saving in time is tremendous.

RECOMMENDATIONS
MATLAB was originally designed to be a general purpose matrix manipulation program, for interactive use, to enable the user to arrive at a quick solution to his problem. As has been shown, it can now be used as a strong tool for solving Signal Processing functions, both in the interactive and batch modes. This has been possible due to the introduction of 'extended' fun<!itions into the MATLAB library. The language can still be further enhanced to rriake it a more powerful utility than what it is now.
The present 'extended' functions return only one matrix and this is a drawback quite often. We need to have 'super extended' functions which would return more than one matrix. It may be possible to construct the 'super extended' functions in a similar way to the extended functions. It is thus recommended that the 'super extended' functions should be created.
It was noticed near the end of the project that calls to the newly added functions would 'bomb out' when the number of arguments in the function call were less or more than the desired number of arguments needed. A routine which returns error messages when such function calls are encountered is recommended to be added to the source of the MATLAB language.
This routine may have to be written in DG assembly code.
MATLAB is an interpreted language, unlike many other languages. This effects the efficiency of the program. It is possible to construct a compiler for the MATLAB language. This would boost its usefulness tremendously.
It would then be definitely a 'super high level' language. Building a compiler is another recommendation of this thesis.