Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Estimation and control of paper machine variables using wavelet packet analysis Chun, John Byung-Kyu 1998

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_1998-0123.pdf [ 11.93MB ]
Metadata
JSON: 831-1.0065014.json
JSON-LD: 831-1.0065014-ld.json
RDF/XML (Pretty): 831-1.0065014-rdf.xml
RDF/JSON: 831-1.0065014-rdf.json
Turtle: 831-1.0065014-turtle.txt
N-Triples: 831-1.0065014-rdf-ntriples.txt
Original Record: 831-1.0065014-source.json
Full Text
831-1.0065014-fulltext.txt
Citation
831-1.0065014.ris

Full Text

ESTIMATION AND CONTROL OF PAPER MACHINE VARIABLES USING WAVELET PACKET ANALYSIS by JOHN BYUNG-KYU CHUN B . A . S c . f E . E . ) , University o f British Columbia, 1995 A THESIS S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in T H E F A C U L T Y O F G R A D U A T E STUDIES DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING  We accept this thesis as conforming 1  to the required standard  T H E UNIVERSITY OF BRITISH COLUMBIA March 1998 © John Chun, 1998  In presenting this thesis in partial fulfillment of the requirements of an advanced degree at the University of British Columbia, I agree that the library shall make it free available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial shall not be allowed without my written permission.  Department of Electrical and Computer Engineering  The University of British Columbia 2356 M a i n M a l l Vancouver, B C V 6 T 1Z4 Canada  Date: March 1998  Abstract T h e thesis discusses a new method in paper machine signal analysis using the wavelet packet transformation. It has been shown that this technique can be applied to three main areas in paper machine control such as the paper machine data analysis including estimation and compression, control of basis weight or moisture variations, and control performance monitoring. Compared to the traditional windowed Fourier transform, other filtering techniques, and wavelet transform, this method offers greater  flexibility  in a n a l y z i n g signals with suitable time-frequency resolution. T h e Best B a s i s A l g o r i t h m has been introduced to obtain a best decomposition tree structure in decomposition of a given signal.  T h e paper machine model for use in C D control has been derived from experimental bump test data. It has been clearly shown that the new method along with wavelet estimator produces results superior to standard analysis methods by reducing the error of estimation. In estimating more complex signals, however, the wavelet packet method is superior than the wavelet analysis since the best decomposition tree structure is chosen enabling the estimation to be more precise in the regions where frequency changes. When applied to industrial data, the wavelet packet estimation has provided more visual results by distinguishing the bump locations more evidently. Moreover, it has been shown that wavelet packet compression is more efficient than the wavelet case since it has achieved more compression ratio by y i e l d i n g the combination of more zero coefficients and lesser mean square error.  The thesis also presents the control scheme of the paper machine as well as an efficient technique to monitor the control performance on-line using both wavelets and wavelet packets. It has been shown that the calculation of performance index using the wavelet packet method is more precise since it offers finer frequency resolution.  ii  Table of Contents Abstract  •  List of Table  vi  List of Figures  •  Acknowledgments  1  2  3  ii  vii x  Introduction  1  1.1  Overview of Estimation Methods  3  1.2  Motivation  4  1.3  Outline of Thesis  5  Paper Machine Model for Use in CD Control  6  2.1  M o d e l i n g of a Paper Machine  6  2.2  Identification  9  2.2.1  Estimating the Interaction Matrix G  11  2.2.2  Identification of Dynamic and Time Delay  12  Estimation Theory  13  3.1  Brief History of Wavelet Theory  13  3.2  Wavelet Analysis  14  3.3  3.2.1  Wavelet Bases  16  3.2.2  Properties  18  3.2.3  Discrete Wavelet transform and Multi-Resolution Analysis  19  Wavelet Packets 3.3.1  21  Wavelet Packet Functions  21  III  3.3.2  The Best Basis Algorithm  3.4  Two-Dimensional Analysis  3.5  some Practical Issues  3.6  De-Noising 3.6.1  :  23 27  .1...  27  ."  29  Threshold Rules  30  3.6.1.1 Threshold Selection Methods 3.7  4  31  Data Compression .  •  Data Analysis  34  4.1  Estimation of Simulated Data  35  4.2  Estimation of Two-Dimensional Paper Machine Data  42  4.2.1  Smooth Waves  42  4.2.1.1 Estimation 4.2.2  '  Signal Containing Several Frequency Components 4.2.2.1 Estimation  4.2.3  4.2.4 4.3  44 49 50  Signal with Bumps and a Streak 4.2.3.1 Estimation  5  32  55 56  Industrial Data Analysis  61  Data Compression  66  CD Control of Basis Weight  69  5.1  Dahlin Controller  69  5.2  C D Profile Assessment  71  5.2.1  Performance Index  73  5.2.1.1 Frequency Characteristics of Wavelet and Wavelet Packet  iv  73  5.3  6  C D Basis Weight Control  77  Conclusion  80  References  83  Appendix A. Gram Polynomial Estimation  87  V  List of Tables 3.1  Wavelet families and associated properties  17  3.2  M i n i m a x thresholds for various sample size n (Donoho and Johnstone [17])  32  4.1  Best M S E of four signal using different methods  36  4.2  M S E with the best option sets of wavelets  36  4.3  M S E with the best option sets of wavelet packets  36  4.4  M S E for the wavelet estimation of the smooth waves  45  4.5  M S E for the wavelet packet estimation of the smooth waves  45  4.6  M S E with best option sets for wavelet estimation (signal containing several frequency components)  4.7  ;  52  M S E with best option sets for wavelet packet estimation (signal containing several frequency components)  52  4.8  M S E with best option sets for the wavelet estimation (signal with bumps and a streak) ....58  4.9  M S E with best option sets for the wavelet packet estimation (signal with bumps and a streak)  58  4.10 Data compression results  68  vi  List of Figures 1.1  A schematic of a paper machine  1  1.2  Scanning path on stationary sheet  2  1.3  A n example of Gram-polynomial estimation: the solid line is the basis weight and the dashed line is the estimated  4  2.1  Bump responses of different paper quality  9  2.2  (a) up: Actuator set point changes (b) bottom: bump responses data  10  2.3  (a) Aligned bump responses (b) C D actuator position and averaged bump response  11  2.4  Estimated first order dynamic response  12  3.1  (a) Gaussian window, (b) Time-frequency resolution of the S T F T (c) Time-frequency resolution of the wavelet analysis, (d) Randomly chosen time-frequency resolution of the wavelet packet analysis. (Each rectangle corresponds  to the localization of a basis  function.)  15  3.2  Haar mother wavelet and its translated and scaled version  16  3.3  Discrete wavelet families. The filter coefficients are used to generate decomposition (or reconstruction) highpass (or lowpass)  filters  17  3.4  B l o c k diagram of one-dimensional D W T  20  3.5  Wavelet decomposition tree  22  3.6  Wavelet coefficients indiscrete and continuous transform  22  3.7  The Haar wavelet packets  23  3.8  Wavelet packet decomposition o f sine with jumps  26  3.9  B l o c k diagram of two-dimensional D W T  28  3.10 Decomposition using D W T and D W P T of two-dimensional image  vii  28  3.11 Data analysis using wavelets  :  30  3.12 Hard and soft thresholding with A, = 0.5  31  3.13 B l o c k diagram of compression and decompression schemes  33  4.1  One-dimensional signals to be estimated, (a) Bumps (b) Pure sine (c) Signal containing different frequency components (d) Blocks  35  4.2  Gram-polynomial estimation  37  4.3  Wavelet and wavelet packet estimation  39  4.4  Wavelet tree structure for signal c  39  4.5  First 50 low M S E resulted using wavelets  40  4.6  First 50 low M S E resulted using wavelet packets  41  4.7  Initial actuator settings and resulting data set (Smooth waves)  43  4.8  Generated paper machine basis weight data  4.9  M S E for estimation using different mother wavelets  "4.10 The smooth wave estimation profiles  :  ,  43 46 47  4.11 The smooth wave estimation profiles (Images)  48  4.12 Cumulative error along C D for the smooth wave  49  4.13 Initial actuator settings and resulting data set (Signal containing several different frequency components)  51  4.14 Generated paper machine basis weight data  51  4.15 The estimated profiles of the signal containing several frequency components  53  4.16 The estimated profiles of the signal containing several frequency components (Images) ...54 4.17 Cumulative error along C D for the signal containing several frequency components  55  4.18 Initial actuator settings and resulting data set (Signal with bumps and a streak)  57  viii  4.19 Generated paper machine basis weight data  ;  57  4.20 The estimated profiles of the signal with bumps and a streak  59  4.21 The estimated profiles of the signal with bumps and a streak (Images)  60  4.22 Cumulative error along C D for the signal with bumps and a streak  61  4.23 The estimated profiles of the industrial data (Images)  63  4.24 The estimated profiles o f the industrial data  64  4.25 Decomposition tress resulted from the Best Basis Algorithm  65  4.26 Truncated wavelet coefficients before and after the thresholding  67  4.27 First 2500 coefficients before thresholding  67  4.28 Images before and after the data compression  68  5.1  Dahlin controller  70  5.2  Different cut-off frequencies or wavelengths for a given actuator width  72  5.3  Approximation and details from wavelet decomposition and their spectrums  75  5.4  Approximation and details from wavelet packet decomposition and their spectrum  76  5.5  B l o c k diagram of basis weight control scheme using wavelet analysis  77  5.6  Control of C D profile: control action took place from scan 10  78  5.7  Performance index during C D control  79  ix  Acknowledgments I would like to thank my thesis supervisor Professor M i c a e l S. Davies for his kind support and guidance during my research. I also would like to thank Professor G . A . Dumont and friends in Pulp and Paper Center as w e l l as in Electrical Engineering Department, in particular, M i n g Zhang, Lahoucine Ettaleb, Roger Shirt, Huiping Chen.  I also thank Jahan Ghosfraniha and Zoran Nesic for their sharing of valuable knowledge. Thanks to T a e - K y u n g K i m , Sinclair C h i , Chris K i m , and Peter K i m for their friendship and encouragements.  Special thanks to my parents who endeavored with me and helped me in everything during my research. Finally, I thank my G o d .  CHAPTER 1 Introduction  A paper machine, with several stages, continuously produces a sheet of paper with a width up to 10 m wide and at a speed of up to 90 km/h [37]. The mixture of about 5% fibre and 95% water by weight is ejected from the headbox unto the moving wire as uniformly as possible. It passes over the moving wire, through the press section, and to the dryer section. The mixture becomes paper as water is removed from the fibre. A t the dry end, where a sheet is wound onto the reel, the water to fibre ratio will be reversed to roughly 5% moisture and 95% fibre. A simplified schematic of a paper machine is shown in Figure 1.1.  Figure 1.1 A schematic of a paper machine  1  CHAPTER 1 Introduction  Sensor off-sheet Machine direction  sensor path  Cross direction Figure 1.2 Scanning path on stationary sheet A traversing sensor mounted on an O-frame measures three sheet properties: the basis weight denned as weight per unit area, the moisture content denned as percent water by weight, and the thickness of the manufactured paper or caliper. In this thesis, the measurement always refers to dry basis weight. A s shown in Figure 1.1, the sensor moves to and fro over the rapidly moving sheet and in each scan collects up to 3000 points of data across the sheet in about 20 seconds. A s a result of the sensor motion across, the moving sheet, measurements are taken at points that would be traced in a zig-zag pattern on the stationary  sheet (See Figure 1.2). T h e actual pattern of points would make an angle o f 1 ° with the sheet edge. Burkhard and Wrist classified basis weight variations into three components, machine direction ( M D ) , cross direction ( C D ) , and random residual variations (R) [5]. The M D component, varying with respect to time, is the variation in the machine direction. It depends mainly on the headbox consistency and pressure changes, which influences the entire width of the sheet. T h e C D component is the variation in the cross machine direction, assumed to be time-invariant, and changes slower than the M D variation. T h e final component, R, contains measurement errors. F o r C D control, the C D profile is extracted from the raw measurement and it is used to manipulate the actuator positions accordingly, which ideally make a uniform distribution o f basis weight across the sheet. T h e scan average o f M D samples during each scan is subtracted from the raw data to leave zero mean C D profiles and noise. T h i s method is called  2  scan  CHAPTER 1 Introduction  average estimation. [28].  1.1 Overview of Estimation Methods The computer control of paper machine requires an efficient estimation method and a proper control scheme. In industrial mills, the most common estimation method is an exponential filtering (EXPO) technique that can be illustrated as follows [11].  y"(f) = (1 -a)y (t-I) n  where x (t) n  + ax (t) n  is the zero mean raw measurement and y (t) n  (1.1)  is the filtered C D profile obtained at  C D position n during scan t and a is the weighting factor, 0 < a < 1. The M D profile can be extracted by averaging y"(r) over one scan period. Generally, the separation of the C D and M D using this method is slow due to the low value of a that is used. For example, if a = 0.2 is used, 20% of the most recent measurement will be added to 80% of the historical value to obtain the trended measurement [36]. Consequently, it would take 10 scans to register 90% of a step-wise profile change.  More advanced filtering techniques were also introduced by Chen [6] and Goodwin and Sin [20]. Algorithms that use a modified least-squares parameter identifier for C D estimation and a Kalman filter for M D estimation are presented by Dumont [18], Wang [37], and Natarajan [31].  Modeling profiles with discrete orthogonal polynomials has been used in Sendzimir steel rolling mills [19]. More recently, Kristinsson also presented an estimation method using Gram polynomials, which are defined discrete equidistant orthogonal polynomials with equal weighting on every data points [25]. An example in Figure 1.3 shows the estimation of industrial paper basis weight data using Gram2  polynomials. The data used in the estimation is C D profile and noise after the scan average of 45.19 g/m is subtracted from the raw profile. The details of this method are given in the Appendix A . In this thesis, 3  CHAPTER 1 Introduction  the E X P O method and the Gram polynomial method are adapted in order to compare them with the wavelet and wavelet packet methods.  Estimation of basis wieght using Gram-polynomial  _i  i  i  0  100  .  i  200  i  __i  300 400 Measurements  :  i  i  i  500  600  700  Figure 1.3 A n example of Gram-polynomial estimation: the solid line is the basis weight and the dashed line is the estimate.  1.2  Motivation On a paper machine, the C D control provides benefits such as increased production, improved  paper quality, reduction of raw material usage, and a decreased number of sheet breaks [5] [10] [25]. The development of computer and networking technologies with fast computation have enabled more powerful estimation systems to be set on-line with sophisticated control loops. In addition, as higher resolution sensors appeared in the market place, measurements of paper qualities with shorter sampling intervals have become available. As a result, the need for more precise estimation and control methods has become evident. On the other hand, data compression becomes more necessary in order to store the large amount of industrial data for future analysis. The performance of the control loops has to be monitored on a regular basis so that an operator acquires suitable information and is able to react to a sudden change of manufacturing conditions effectively. The main objectives of this thesis are first to develop a novel method to  4  CHAPTER I Introduction  estimate C D components of the paper basis weight from noisy raw data, second to study the control of a paper machine using the identified paper machine model, and finally to create an efficient technique to monitor the control performance on-line.  1.3 Outline of Thesis In Chapter 2, a paper machine model is developed and then parameters of the model are identified from the actual input-output measurements. Estimation theory using wavelets and wavelet packets are presented in chapter 3. Also in this chapter, the noise reduction procedure and the compression method will be illustrated. Chapter 4 presents the error analysis of simulated data with different methods such as the E X P O analysis, wavelet estimation, and wavelet packet estimation. The estimation and compression of industrial basis weight data will also be displayed. Chapter 5 presents the control scheme of the paper machine as well as the C D profile performance assessment methods using both wavelets and wavelet packets. Finally, conclusions are drawn in chapter 6.  5  CHAPTER 2 Paper Machine Model for Use in CD Control  A C D basis weight control system uses the measured profile read at the sensor to manipulate the slice lip actuator settings.,The computer analyzes the scanned data, calculates proper control actions, and then generates the actuator set point changes accordingly. The control simulations in later chapters are carried out using the paper machine model developed in this section.  2.1 Modeling of a Paper Machine The paper machine model which gives a relationship between the actuator settings and the C D basis weight responses can be derived mathematically using experimental bump test data. The discrete time mathematical model for a paper machine includes actuator dynamics, time-delay, and interactions [25]. The major assumptions made in modeling the actuator response are  (i) The responses to all actuators have the same dynamics (Wilhelm and Field, 1983) (ii) The response shape is symmetric. (iii) The total profile response is the sum of the responses to each other. (iv) A sampling time is known as  T. s  Usually, all actuators are represented by the same first order lag dynamic with the gain k  a  time constant p  a  and the  as follows •  K — r  (2.1)  (i-/vr ) ]  Since the basis weight scans are taken some distance down the machine direction from the actuators, the  6  CHAPTER 2  Paper Machine Model for Use in CD Control  model should include a time-delay q~ corresponding to d • T . It is assumed that the speed of a paper d  s  machine is constant as is the time delay.  Movement of an actuator causes the basis weight to change along the profile. This profile response spreads to the location of neighboring actuators resulting a coupling between actuator spaces. Therefore, assuming the same spatial response for all actuators, a paper machine should be modeled by a banded diagonal matrix called the interaction matrix G as shown in equation (2.2). Except near edges of the matrix, each column consists of zeros and the bump response data streams {P] ---P\, Pq, P\, ... \. l  Pjl  The  height of the matrix is the number of measurements, and the width of the matrix is equal to the number of actuators.  Po P\  '•' Pk  ••  0  0  Pi Po Pi -  G =  P  • 0  Pi k  P  k  o  \ P  x  W  •••  Pi Po Pi  0 ••• 0 p  k  •••  P  l  P  o  As discussed in [25] and [34], some problems exist in building the interaction matrix G . The interaction matrix should be constructed to be nonsingular, otherwise the condition number will grow infinite and it can lead to a poor performance of C D control. The condition number is defined as the measure of the severity of expected control difficulty [26]. More specifically, Kristinn defines the condition number of G as the ratio between its largest and smallest singular value [25]. Singular values are obtained from:  eig{G'  • G)  (2.3)  It can be shown that i f the control number grows to infinite, the machine is uncontrollable because  CHAPTER 2  Paper Machine Model for Use in CD Control  corresponding actuator positions to the response profiles cannot be determined. (If singular, G is not invertible.) The interaction matrix should be built with care avoiding its singularity for better control performance.  The model parameter uncertainties exist due to the changes occurred in machine speed, pulpwood quality, and so on during normal operation. This uncertainty can also lead to poor control performance [26]. The measurements are all aligned with the actuator locations, only if the number of measurements is the integer multiple of the actuator number.  Consequently, the overall dynamic model for the paper machine is given as  y(0 = G  u(t) 1  " /vr  (2.4)  1  The actuator responses vary with fibre quality basis and machine speed. Figure 2.1 displays how a unit change of the actuator affects the basis weight response for different types of paper machines. As the paper weight increases from newsprint to board, the response peak has smaller amplitude, however, the coupling region becomes much greater. The different C D responses downstream from p  Q  equation (2.2) are reported in [26].  to p  k  in  CHAPTER 2  Paper Machine Model for Use in CD Control  Measurements Figure 2.1 Bump responses of different paper quality [10]  2.2 Identification Basis weight profile, measured at the O-frame, is controlled by the slice lip actuators at the headbox, where control actions are influenced by several control loops. Identification of paper machine response model is desirable before building any control loops. The parameters of the model in equation (2.4), which depend on machine speed and grade, can be identified using the input-output measurements of bump test data. As mentioned before, actuator dynamics, time delay and interaction matrix are three components that make up a paper machine model and they are to be identified.  The deposition of coating material on to a base sheet is controlled by screw adjustments to shape the coater blade, and leads to a similar response model to that of basis weight response due to headbox slice lip adjustments. In this thesis, the model is first identified using the bump response data from the coater, and the identified model is used as a paper machine model. However, as explained later, it is required to change a parameter value of time-delay so that the identified model from the coater gives a similar response to that from the headbox. The following bump test data, taken as provided by Measurex Devron Inc. is used as an example to identify the response model that will be used in the simulation.  9  CHAPTER 2  Paper Machine Model for Use in CD Control  The parameters of test are as follows:  (i) (ii)  Number of minislice (high resolution) weight measurement points: 300 Number of date scans gathered: 47  (iii)  Number of actuators: 100  (iv)  Actuator spacing: 3 mm  (v) Coat weight target 7 gsm (vi) (vii)  Positive bump actuator locations: 10, 30, 50, 70, 90 Negative bump actuator locations: 20, 40, 60, 80 In the coater, the set points of ten actuators were manipulated by one step at the specified locations  for 30 scan period. The sensor collected 300 measurement points per scan for 47 scans. Five positive and four negative steps at the actuators, as shown in Figure 2.2 (a), created 10 bumps at the output measurement. As Figure 2.2 (b) reveals, the profiles of all bump responses have similar shapes.  T  Actuator settings  and t h e i r responses  i  i  i  i  i  i  r  100 Actuator p o s i t i o n  0 "  50  100  150 Measurement  >  200  250  >  Figure 2.2 (a) Up: actuator set point changes, (b) Bottom: bump responses data  10  300  CHAPTER 2  Paper Machine Model for Use in CD Control  2.2.1 Estimating the Interaction Matrix G As shown in Equation (2.2), the interaction matrix G describes the cross machine direction steady state actuators coupling in one simple matrix. The coefficients of G can be approximated by averaging the positive bump responses. Each bump in Figure 2.2 (b) is relocated so that the centers of the responses are aligned, and they are drawn in Figure 2.3 (a). The databoxes reveal that the peaks of every bump responses are moved to left by 1 measurement point from the bump locations. A l l bump responses are then averaged. The resulting averaged bump response is plotted in Figure 2.3 (b) by assuming the actuator position at zero. It is seen that the single actuator movement affects about 5 actuator positions. The figure reveals that the coefficients of G i n this case are [p , p, ...,p ] Q  Centered  k  = [1.71, 1.42, 1 . 0 9 , 0 . 5 7 , 0 . 1 , 0 . 0 1 , - 0 . 1 5 ] .  bump responses  Actuators'and  databoxes  B a s i s weight  CD a c t u a t o r p o s i t i o n s  Figure 2.3 (a) The aligned bump responses, (b) C D actuator position and averaged bump response  11  CHAPTER 2 Paper Machine Model for Use in CD Control  2.2.2  Identification of Dynamic and Time Delay In M D , the distance between the location of actuators and the sensor creates the time delay. The  M D step responses are extracted at the C D location corresponding to the peak bump response magnitude. Then, they are averaged and approximated by the first order dynamic. The result gives the value of time delay d = 1 scans, K  a  = 0.0384, and p  a  = 0.732 . The first order approximation is very close to the  real shape of the step response as shown in Figure 2.4. The time-delay from the coater to the sensor is found to be less than one since the coater is located near to the sensor. The identification was performed under the assumption that the coater bump responses are similar to those from headbox actuators. Because of the longer distance from a headbox to a sensor, the time-delay of 1 is not an appropriate value for the paper machine model. To generate more sounding responses, the time-delay d = 2 will be employed in the simulation.  Overall, the paper machine model is identified using the coat weight data as  y(0 = G  0  5  10  2  1 - (0.732)$"  Step response  1.5  0  (0.0384)g-  u(t)  estimated i n f i r s t  15  20  (2.5)  order  ;25  Scans  Figure 2.4 Estimated first order dynamic response  12  u 10  35  CHAPTER 3 Estimation Theory  On a paper machine, a sensor mounted on an O-frame collects measurements as it travels to and fro over the moving paper. The most common measurements are the dry basis weight and the moisture content of sheet. In addition to the measured value, the basis weight data contains high frequency sensor noise as well as random environmental disturbance. Control action and other operations are dependent on accurate measured data, and it is important to obtain data with a minimum of noise. More precise sensors and faster computer systems provide efficient ways of acquiring correct data. With the help of improved data measurement and processing techniques, including de-noising and compression techniques, it is possible to improve our knowledge of process variations. In this chapter, estimation and compression methods using wavelets and wavelet packets are presented. Traditional methods, for example Grampolynomial and Expo-filter approximations, are described in the Appendix A and in the previous chapter respectively.  3.1 Brief History of Wavelet Theory The underlying theory of wavelets can be traced back to the nineteenth century when Joseph Fourier presented the frequency analysis method that is now known as the Fourier analysis [27]. Later the short-time Fourier transform (STFT) was developed for time-varying frequency domain analysis. More recently, wavelet analysis, recognized as a simple tool with a great variety of possible applications, has been applied in many different fields of mathematics, physics, and engineering [14]. The term 'wavelet' was first introduced by Alfred Haar in 1909 [33]. However, it was not until 1982, when Morlet published his work, that wavelets became a method for signal analysis [29]. Grossmann and Morlet [21] in 1984 carried out the detailed mathematical study of the continuous wavelet transform, followed by the work in the discrete case by Daubechies, Grossmann, and Meyer [13]. In 1988, Daubechies [12] constructed a set of wavelet orthonormal base functions that drew much academic attention to the subject of wavelets.  13  CHAPTER 3  Estimation Theory  Moreover, wavelet packets were constructed by Coifman and Meyer in 1991. More references for wavelet packets can be found in Coifman, Meyer, and Wickenhausen [7] and in [32]. Currently, researchers are actively searching for different applications in various areas of signal processing, image processing, data compression, and so on.  3.2 Wavelet Analysis The classical Fourier analysis represents a signal in terms of trigonometric sine and cosine functions in different frequencies, allowing analysis in the frequency domain. Since the trigonometric basis functions extend over all time, this analysis is unsuitable for dealing with signals whose frequency contents are changing in time. To address this problem, Gabor (1946) proposed the Short-time Fourier transform, or windowed transform [33], as  STFT(x,f)  where g(t)  = \x(t)g*(t-x)e- dt 2jnft  (3.1)  is a window that truncates sinusoids to form basis functions with a particular width. For  example, Figure 3.1 (a) shows a Gaussian window g(t) that is given as  8 (t) a  = ^ ~ '  2  /  (  4  C  °  (3-2)  This representation allows changes in signal spectrum to be tracked over time. For the STFT, the choice of g(t) fixes the resolution over the entire time-frequency plane since the same window is used at all frequencies. Figure 3.1 (b) shows the resolution of the STFT in a time-frequency plot. Here, it is clear that the time and the frequency resolutions are fixed in all ranges. To compare resolution with the wavelet and the wavelet packet bases in Figure 3.1 (c) and (d), some definitions are to be made. Given the window function g(t) and its Fourier transform G(f),  A  the time and the frequency resolutions are defined as  ,2  \f Wfdf 2  \\G{ffdf  14  CHAPTER 3  and At  2  \t\g(tfdt  —  = J  Estimation Theory  (3.4)  \\g{t)\dt  where the denominator is the energy of g(t). At and Af are related through an uncertainty principle [33], similar to the Heisenberg inequality which illustrates that the product of time and frequency spreads are lower bounded as A,.A/>;1  (3.5)  As seen from equation (3.5), the area of each rectangle in Figure 3.1 (b), (c) and (d) which is proportional to At • Af is near constant. This implies that a signal can be analyzed either with good time resolution or with good frequency resolution, but not both.  Figure 3.1 (a) Gaussian window, (b) Time-frequency resolution of the S T F T (c) Time-frequency resolution of the wavelet analysis, (d) Randomly chosen time-frequency resolution of the wavelet packet analysis. (Each rectangle corresponds to the localization of a basis function.)  As shown in Figure 3.1 (c), for wavelet analysis the resolutions At and Af vary in the timefrequency plane, resulting in time resolution that is finer at high frequency and frequency resolution that becomes better at low frequency. Therefore, in comparison to the Fourier transform, the wavelet transform  15  CHAPTER  3  Estimation Theory  has the flexibility to estimate a signal with time domain discontinuity and sharp spikes, giving better timefrequency localization. The wavelet packet case will be discussed later in this chapter.  3.2.1 Wavelet Bases Wavelet analysis uses wavelet functions as basis functions, where the STFT employs truncated sinusoids. Wavelets are fixed basis functions, localized in time and last for a few cycles. In wavelet decomposition, a mother wavelet which is defined below is dilated and shifted to represent various signals of different frequency and shape. For example, Haar wavelets in different scales and positions, as shown in Figure 3.2, are constructed from the basic building block, ^(x), known as the "mother wavelet". The discrete operation uses the dyadic scaling index j and integer translation index k, which can be expressed as  "Vj (x) = 2-i M{2-ix  -k),  /lx  k  jeZ,keZ  (3.6)  The good time localization of Haar wavelet is traded with the poor localization property in the frequency domain. Several smoother wavelet bases, that result in better frequency localization, have been introduced. For example, Meyer, Daubechies, Symlets, and Coiflets wavelet families are shown in Figure 3.3. Similar to Equation (3.6), the mother wavelet in continuous analysis is defined in [30]. Wavelet analysis in this thesis refers to only discrete operation, unless it is mentioned otherwise.  4  j=4, k=13  2  j=3, k=2 j=0. k=0  0 -2 -4 0 (a) Haar basis function  0.2  0.4  0.6  (b) Haar wavelet examples  *P(JC)  Figure 3.2 Haar mother wavelet and its translated and scaled version  16  0.8  1.0  CHAPTER 3 Estimation Theory  0.8  1  2  3  4  2  coifl  4  6  8  2  4  coif2  6  8  10  12  5  coif3  10  15  5  10  coif4  15  coif5  Figure 3.3 Discrete wavelet filter families. Those filter coefficients are used to generate decomposition (or reconstruction) highpass (or lowpass) filters.  Meyer Tiller length Compactly supported orthogonal  •  Symmetry  lliiar  Diiuhecliies  Symlets  Coiflets  2  2N  2N  6N  • •  •  •  •  •  • •  • •  N  N  2N  •  •  •  •  Asymmetry Near symmetry  •  Smoothness # o f vanishing moments o f Exact reconstruction  ^(t)  •  •  Table 3.1 Wavelet families and associated properties (N is the order of wavelet e.g. N = 3 for db3.)  17  20  CHAPTER 3  Estimation Theory  3.2.2 Properties The qualities of different types of wavelets vary one from another, and a suitable basis wavelet should be selected to obtain the best results in estimation. The wavelet families and their associated properties are illustrated in Table 3.1. In this section, those wavelet properties are denned as follows;  • Compact support. If the wavelet  is compactly supported, the wavelet filters, that will be  explained in the next section, are finite impulse response (FIR) filters. Consequently, the wavelet transform becomes finite. • Symmetry. Symmetry is essential for the niters to have linear phase; this property is required to avoid phase distortion. • Smoothness: The smoothness of wavelets plays an important role in compression and estimation (de-noising) applications. If the wavelet is not smooth, the error can be easily detected from the filtered signal except for the signal which has sharp edges. M degrees of smoothness implies that the M  th  derivative of a function is continuous at all points. A higher degree of smoothness offers bet-  ter frequency localization of the niters. • Number of vanishing moments. A sequence {\\f } is said to have M vanishing moment if k  •XV*'*' = °>  f o r  0</<M.  (3.7)  k  Large number of vanishing moments give smoother ^ ( f ) . Thus, this property and the smoothness are related each other. This property for each mother wavelet is denned in Table 3.1. • Time-frequency localization. If *?(£) decays rapidly, the wavelet has good time localization. Good frequency localization is achieved if ^ ( c o ) decays fast as CO —> °° and if ^((fl) is sufficiently flat near co = 0. • Orthogonality. The wavelets given in equation (3.6) are orthogonal to each other in all scale:  18  CHAPTER 3 Estimation Theory  X ^ C O ^ C O  = S(j-J)-S(k-K)  (3.8)  2  This property links the L norm of a function directly to the norm of its wavelet coefficients by ll/ll =  /5>7*  (3-9)  If orthogonal wavelets are used, the fast wavelet transform is a unitary transformation. Under this condition as in equation (3.8), numerical calculation is stable since the transformation will not increase error in initial data. This is an important property in numerical calculation.  More details on wavelet properties are presented in [24] [27] [30] and [38].  3.2.3 Discrete wavelet transform and Muti-Resolution Analysis For analyzing discrete-time samples in paper machine data, the discrete wavelet transform (DWT) is used in this thesis with Mallet's multi-resolution analysis (MRA). In M R A , the mother wavelet *¥(x) is determined by the highpass filter G, which produces the details of the wavelet decomposition. Moreover, the scaling function <!>(x), determined by the lowpass filter H , gives the approximations of the wavelet decomposition. Thus, at level j, the signal Cj_ 1  l  l  is decomposed into the approximation c • and the k  detail dj . I is 1 to the length of the signal. k  Cj,k = < / ( * ) . * * * ( * ) > = 5 > / - 2 * - ^ _ ,  d  Jtk  = </(*),¥;,*(*)> = 2 > / - 2 * - c ; -  where the filter coefficients are related as g  jtk  (3-10)  U  (3.11)  = (-1 ) • h _ . The inner product is defined as k  k  (f(x),V (x))  f /  {  k  = jf(x)V (x)dx jk  19  *  (3.12)  CHAPTER 3 Estimation Theory  Figure 3.4 illustrates the multi-level decomposition and reconstruction schemes of DWT. In the decomposition scheme, the sequence f(n) with a length n is first filtered using H and G and then downsampled by a factor of 2. One level decomposition produces the n/2 approximation coefficients C j and k  the n/2 detail coefficients d, . h  Details  H H  •  f(n)-  H  h(n)  Appro.  •f(n)  H*  12  Decomposition Scheme  T2  h(n)  Reconstruction Scheme  Figure 3.4 Block diagram of one-dimensional DWT With details and approximations, perfect reconstruction is achievable using the quadrature mirror reconstruction filters H * and G * , which is summarized in the following equation,  hk ~ X  C  k-2l j+\,l  h  C  +  Sk-2l j+\,l d  (3.13)  Here, a system formed by the low and high pass decomposition filters H and G , together with their associated reconstruction filters H * and G * is called quadrature mirror filter. The down-sampling in decomposition introduces a distortion called aliasing. However, if the filters are chosen carefully to relate the decomposition and the reconstruction,.it is possible to cancel out the effects of aliasing. Thus, perfect reconstruction is possible. Those chosen filters are quadratic mirror filters. This breakthrough was made possible by the work of Ingrid Daubechies and others. The shape of the wavelets in Figure 3.3 were determined by various kinds of quadrature mirror decomposition filters. The technical details of how to design these filters are discussed in [35].  20  CHAPTER 3  Estimation Theory  As an example, a sine wave with two jumps in Figure 3.6 (a) is decomposed using the M R A . The decomposition tree, in this ease up to level 5, is shown in Figure 3.5. The detail coefficients at each level are plotted in Figure 3.6 (b) after the D W T is performed to level 5. At level 1, the wavelet coefficients contain the high frequency part of the signal. The analysis clearly located the jumps at this level. At higher levels, lower frequency components are extracted as shown in the figure. In Figure 3.6 (c), the wavelet coefficients of the continuous wavelet transform are displayed just to,compare it with the discrete case.  3.3 Wavelet Packets As shown in Figure 3.1 (b) and (c), the time-frequency resolution for different basis functions such as the STFT basis and the wavelet basis are predetermined. For efficient estimation of different types of signals, the basis of the signal, and so the time-frequency resolution, should be selected in a signaldependent manner. Discrete wavelet packet transform (DWPT) to be discussed in this section, the generalization of wavelets bases, offers greater flexibility because it allows the time-frequency resolution to change in a signal decomposition. Figure 3.1 (d) shows one possible wavelet packet basis among many possible options. A n adaptive procedure for choosing the most suitable set of basis functions is described later in this section.  3.3.1  Wavelet packet functions This section introduces wavelet packets following the appearance of [27]. A wavelet packet  function is indexed by three parameters: W -  m  (x).  As in wavelets,./ can be interpreted as a scale  k  parameter and k as a position parameter:  *W(*)  = 2- W (2- x-k) J/2  J  m  (3.14)  The extra index m = 0, 1,... is called the oscillation parameter. The wavelet packet functions W  (x) is defined as  21  CHAPTER  3  Estimation Theory  (0,0)  (1,0)  (1,1)  (2,0)  (2,1)  (3,0)  (4,0)  ,0)  (3,1) (4,1)  (5,1) Figure 3.5 Wavelet decomposition tree  (a)  0  S i n e wave w i t h jumps t o be a n a l y z e d .  100  200  (b) D i s c r e t e T r a n s f o r m ,  100 (c)  200  300  Continuous Transform,  100  300 wavelet  200  wavelet  300  400  500  coefficients.  400  500  coefficients.  400  Figure 3.6 Wavelet coefficients in discrete and continuous transform  22  CHAPTER 3 Estimation Theory  2N-1  W Jx) 2  = 2 S  h W (2x-k) k  m  (3.15)  k=0  2N-  W  2m  +  (x)  l  1  = 2 X  g W (2x-k) k  (3.16)  m  /t = o  where W (x) Q  =  O ( x ) is the scaling function and  length of the filters, h  k  WJ(JC)  =  *¥(x) is the wavelet function. 2N is the  and g . The Haar wavelet packets for m= 0 to 7 are computed and plotted in k  Figure 3.7. There are many possible combinations of grouping wavelet packet bases as a result of the added degree of freedom. The criteria for the selection of the best basis and wavelet coefficients will now be discussed.  0.5  1  0  0.5  Figure 3.7 The Haar wavelet packets  3.3.2  The Best Basis Algorithm In wavelet analysis, a signal is split into an approximation and a detail component. Then the  approximation itself is split into the next level approximation and the detail. This procedure is repeated  23  CHAPTER 3 Estimation Theory  until it reaches the desired level as illustrated by the tree structure in Figure 3.5. In wavelet packet analysis, however, both the detail and the approximation can be decomposed. The completed decomposition tree structure is shown in Figure 3.8 (a). Coifman and Wickerhauser introduced the Best Basis Algorithm to suggest the best decomposition for a particular signal [8]. This algorithm finds the most efficient tree structure for a given signal in the sense that the signal will be most economically represented. Before splitting at each node, the algorithm calculates an entropy value (to be defined in this section) and decides whether further splits are necessary in order to yield a minimum-entropy decomposition tree. Entropy is a measure of concentration, describing information related property for an accurate representation of a given signal [27]. If entropy value is small, more coefficients are negligible. A n example of this procedure is given at the end of the section. The information cost functions that calculate entropy are now defined; other entropy definitions are presented in [2] and [41]. Four alternative entropy-like measure of a sequence are: (i) Shannon entropy. The non-normalized entropy of a sequence s = {i -} depends upon the loga1  rithm of the squared value of each signal sample and defined by H  l I \ i\  = -5>>g(x, )  (3.17)  2  s  2  S  where X; =  and i f x = 0 then set xlogx  = 0 . It has been shown that e x p / / „ is propor-  , INI tional to the number of coefficients needed to represent the signal within a fixed mean square error [9] [41]. (ii)  The number of samples above a threshold level. The number of samples in a sequence s = {S;} whose absolute value exceeds chosen threshold X is counted, giving the number of coefficients needed to transmit the signal to a precision X.  (iii)  Concentration in l  P  s = {5.} , the f  norm. The f  norm with 1 < p < 2 is calculated. For a sequence  norm with chosen value of p is ^ K i ^ -  T  n  e s m a l l e r  m  e  ^  n  o  r  m  o  f a m  n  c  "  i  tion for a given energy level, the more concentrated that energy is into a few coefficients [41].  24  CHAPTER 3  (iv)  Logarithm of energy.  Estimation Theory  The log energy is defined as H  = 2log(s,: ), 2  s  (3.18)  setting log 0 = 0 whenever necessary. This is an additivity-type of function.  For example, using the complete tree representation of the wavelet packet decomposition shown in Figure 3.8 (a), the sine wave with jumps is decomposed up to level 3. Using 8 terminal nodes from (3,0) to (3,7), the wavelet coefficients are displayed in Figure 3.8 (b) with the high frequency coefficients at the top and the low frequency towards the bottom. Using Shannon entropy, the calculated entropy values are shown at each node in Figure 3.8 (a). The Best Basis Algorithm first splits all approximation nodes to the desired level. Here, the approximation nodes are (1,0), (2,0) and (3,0). Then at each detail node where a split is possible, (1,1) and (2,1) in this case, the entropy value is compared with the sum of the next level entropy values. At (1,1), for example, the entropy 0.3007 is compared with the value 0.3385 that is the sum of entropies at (2,2) and (2,3). Since the entropy at (1,1) is smaller than the compared entropy, the algorithm will not split this node. Similarly all detail nodes are examined. The final result of the Best Basis Algorithm is given in Figure 3.8 (c) with the corresponding wavelet coefficient plot. In this figure, it is clear that fewer coefficients are used in the representation of the signal using wavelet packets. For an n N sample signal, there exist 2  different possible tree structures or bases in the library. In Chapter 4, it will  be shown that the Best Basis decomposition result is a more compressed representation following thresholding.  25  CHAPTER 3 Estimation Theory  (0,0) -636.8562 -637.1569  0.3007  (1,0)  (1,1) 0.5493  -637.7061  (2,0)  (2,1)  -638.5740 10.8679  (3,0)  0.1003  (3,2)  (3,1)  0.2495  0.0890  (2,2) 10.4490  (3,3)  j 0.0746  (3,4)  (2,3) 0.0144  0.0880  (3,5)  (3,6)  0.2115  (3,7)  (a) Complete wavelet packet decomposition tree with entropy S i n e w i t h jumps.  i 0.5 0 -0.5 -1 -1.5 0  100  200  300  400  500  I  100  200  300  400  500  (b) Wavelet packet coefficients (0,0)  (1,0)  (1,1)  (2,0)  (3,0)  (2,1)  (3,1)  100  (3,2)  (3,3)  200  300  400  500  (c) The best tree and its decomposed coefficients Fi gure 3.8 Wavelet packet decomposition o f sine with jumps  26  CHAPTER 3 Estimation Theory  3.4 Two-dimensional Analysis The two-dimensional D W T and the D W P T algorithm can be constructed by two separate onedimensional transforms. A s shown in Figure 3.9, the rows of the two-dimensional matrix c • are decomposed first and then the columns of the resulted sets are decomposed, resulting 4 submatrices. In Figure 3.10, the wood image is decomposed to level 3 using D W T and DWPT. At the first level of decomposition, the transformation gives 4 sets of coefficients. As mentioned before, at the next stage the D W T operates on the approximation only, where the DWPT is performed using both the approximation and the details. Thus, the fixed pattern of decomposition in the left hand image of Figure 3.10 is replaced by the right hand decomposition in which component of the mosaic may be further transformed.  3.5 Some Practical Issues In practice, the computational complexity of the fast wavelet decomposition or reconstruction algorithm is O(n) for an n sample signal, while the fast Fourier transform requires In the wavelet packet case, the computational cost increases to  0(nlog2  n  )  ar  0(nlog^n  )  operations.  >d the search for the best  basis uses an additional O(n) operations [30].  Since the length of the filters H and G is larger than 2 and the signal sequence is finite, a boundary problem occurs. The filter lengths for different wavelets are defined in Table 3.1. In the convolution of the signal and the filters, some sequence values have to be padded at the ends of signal. There are three usual ways to handle the boundaries:  symmetric, periodic,  and zero padding. The symmetric boundary approach  involves reflecting the function of interest about the boundaries. In the periodic case, the signal is repeated at the edge. Finally, the simplest way is to pad zeros at the edge. Periodic and zero padding solutions may introduce a discontinuity at the edge, which leads to poor estimation results. In this thesis, the first approach will be used. Other solutions to the boundary problem, such as Meyer's boundary wavelets, dyadic boundary wavelet, and orthogonal wavelets on the interval, are discussed in [14] [24] and [32].  There is no formal criteria exist for selecting the best mother wavelet. Nevertheless, a  priori  knowledge of the nature of the signal may be used to select the mother wavelet based on its smoothness  27  CHAPTER  3  Estimation Theory  [2].  columns  H  rows  columns rows  G  columns  Figure 3.9 B l o c k diagram o f two-dimensional D W T  2D Data  200  250 50  100  150  200  2D Wavelet Decomposion  250  2D W. Packet Decomposition  200  250 50  100  150  200  250  50  100  150  200  Figure 3.10 Decomposition using D W T and D W P T o f two-dimensional wood image  28  250  CHAPTER 3 Estimation Theory  3.6 De-Noising Traditional noise removing techniques remove high frequency noise by low-pass filtering, which smooths the function and also removes high frequency signal components. As discussed in Section 3.2 and 3.3, wavelet and wavelet packet analyses are superior to these de-noising methods, especially in estimating the signals with jumps, spikes and other non-smooth features. The wavelet-based approach to de-noising is non-linear since it is based on thresholding of coefficients. In the wavelet packet framework, the de-noising approach is identical to that of the wavelet framework, and is illustrated as follows. The model for the noisy signal is assumed to be  s = / . + a e , . , i = 0, . . . , n - l t  |  .  '  (3.19)  where £• ~iid N(0,1) is a Gaussian white noise, and a is a noise level. In extracting the signal / , from noise, the following two goals are suggested in the approximation, v .  /: • To represent  in terms of only a relatively small number of wavelet coefficients (i.e. achieve  compression) • To minimize the expected mean-square error n-1  R(f,f) = ^Eil-ff.  (3.20)  i= 0  Donoho and Johnstone in 1992 have proposed a de-noising procedure in three steps [17]. First, decompose the discrete signal / = { / , , . . . , / „ } into wavelet coefficients co at level N . Second, apply the thresholding technique to the detail coefficients and produce the set CO*. Third, reconstruct / using co* which consists of the original approximation coefficients of level N and the modified detail coefficients of level from 1 to N . The procedure is also displayed in Figure 3.11.  29  CHAPTER  3  Estimation Theory  They explain their motivation to propose this principle "selective wavelet reconstruction" as follows [17]. First, the decomposition of the real signal / is concentrated into small number leaving large value of nonzero coefficients. Second, the wavelet transform of the white noise given in Equation 3.19 is also a white noise which is evenly spread over all coefficients.  CO w  1: Decomposition  CO * w  2: Thresholding  w  3: Reconstruction  w  Figure 3.11 Data analysis using wavelets  3.6.1 Threshold Rules Wavelet decomposition produces two sets of wavelet coefficients: approximations and details (See Figure 3.5). Under the assumption given in Equation (3.19), deleting small coefficients in details does not affect the original signal significantly. Thresholding refers to a way of cleaning out unimportant details considered to be noise (i.e. white noise). Figure 3.12 illustrates hard and soft thresholding of the straight line s for a fixed threshold X = 0.5 . The hard and soft thresholding strategies are:  .hard  [0,  d <X  d,  d >X  jk  o f t  _ fsign(d )(\d \ jk  jk  1  jk  - X), 0,  30  (3.21)  jk  if\d \ > X jk  if\d \<X jk  (3.22)  CHAPTER 3 Estimation Theory  0.5 " h  -0.5  HARD  SOFT  Figure 3.12 Hard and soft thresholding with A = 0.5  3.6.1.1 Threshold Selection Methods Most of the work on thresholding is due to Donoho and Johnstone, who have introduced several threshold selection methods that can be grouped into two categories: dependent thresholding  global  thresholding  and level-  [27] [32]. The former chooses a single value A to be applied globally to all detail  coefficients. The latter adapts a different threshold value A • for each wavelet level j. This data adaptive wavelet thresholding is used in handling of a nonwhite noise sequence. The thresholding value A is defined by  A  y  =  djl  (3.23)  where the median of absolute deviation (MAD) is given as  (3.24)  and j is fixed to 1 for global thresholding. The following threshold selection methods, implemented in Matlab™ toolbox, are used to find A . The methods are given as:  31  CHAPTER 3 Estimation Theory  (i)  Minimax thresholding depends on the sample size n. X is chosen to minimize the upper bound of mean square error involved in estimating a signal [32]. No closed form is exist for this method. Its numerical approximation for various sample sizes are given in Table 3.2.  n  X  n  X  64  1.474  2048  2.414  128  1.669  4096  2.594  256  1.860  8192  2.773  512  2.047  16384  2.952  1024  2.232  32768  3.131  Table 3.2 Minimax thresholds for various sample size n (Donoho and Johnstone [17])  (ii)  Universal thresholding The thresholding value is known as X = j2\og(n)  where n is the sam-  ple size. The universal threshold value X is substantially larger than that in the minimax case, so the reconstruction is smoother. This method is also called VisuShrink. (iii)  SURE thresholding is based on Stein's Unbiased Estimation of Risk. The threshold value X is derived by minimizing and estimating the expected error value [32]. (See Equation (3.20)).  (iv)  Hybrid SURE thresholding If the signal to noise ratio is small, the S U R E estimate may be very noisy. If such a situation is present, universal thresholding is used; otherwise the S U R E criterion is used [32]. As will be seen in Chapter 4, the differences in results using those methods when applied to paper  machine data will not be significant and the threshold method is possibly fixed to a particular case in simulations.  3.7 Data Compression In industry, it is desirable that the estimated paper machine data be compressed and stored for a future analysis. The wavelet-based data compression schemes have been reported by researchers, ranging  32  CHAPTER 3 Estimation Theory  from simple entropy coding to more complex techniques such as vector quantization [2] [23]. In the wavelet transform, the compression feature appears naturally due to the large number of zero coefficient that are present after thresholding. The fact that a signal can be approximated with the decreased number of wavelet coefficients leads to the wavelet data compression. The wavelet compression of a sequence s = {S[j} is shown in Figure 3.13. The compression is accomplished by applying a wavelet transform or a wavelet packet transform with a thresholding method as described in the previous section. Coefficient amplitude quantization is optional, but leads to much higher compression ratio. The resulting coefficients are encoded to be stored in a data file. The compressed signal is recovered to form a sequence J = by the procedure illustrated in bottom of Figure 3.13.  COMPRESSOR  Transform*  w  Quantizer  Lossy  Lossy  w  Encoder  Lossless  Storage Device  Inverse Transform* U-  DeQuantizer 4  Decoder  DECOMPRESSOR  Figure 3.13 Block diagram of compression and decompression schemes In the simulations described here, the quantization and dequantization and the encoding and decoding will not be implemented. However, in this thesis, to measure the compression performance, nonzero coefficients are counted after thresholding the data, and then compression ratio is calculated. In the reconstruction process, the mean square error, R , defined in equation (3.20) is also used to assess the performance. The compression procedure using wavelet packets is identical to that of wavelet analysis. The simulation results are displayed in the next chapter.  33  CHAPTER 4 Data Analysis  Paper machine data includes bumps, spikes, and smooth variations at different frequencies. In this section, various methods for estimation of scanned paper machine data are examined and compared. Data compression using wavelet and wavelet packet analyses will also be discussed. In one-dimensional analysis, three different methods are considered: Gram-polynomial basis function, wavelet, and wavelet packet. The E X P O method is compared with wavelet analysis in the two-dimensional estimation. For simulation, these methods are tested using Matlab™. Several software packages are available for Matlab™ wavelet use including: WaveLab (Stanford University) [3], Wavelet Toolbox by the Math Work Inc. [27], Rice Wavelet Toolbox (Rice university), and the U B C PPC wavelet toolbox [38]. The Matlab Wavelet Toolbox is used in this thesis and has been adapted to process two-dimensional paper machine data for estimation and compression.  Estimation using the E X P O and Gram-polynomial methods is simple since only one parameter need be manipulated; a weighting factor in the E X P O algorithm and the number of polynomial order for gram decomposition. However, wavelet and wavelet packet analyses offer many different options to examine different signals. The parameters used in simulation for both wavelet and wavelet packet analyses are given briefly as follows:  • Boundary handling: symmetrical padding at the edge • Level of decomposition: 1 to 5 • Threshold type: Hard, Soft • Threshold selection: Minimax, Universal, SURE, SURE hybrid  34  CHAPTER 4 Data Analysis  • Thresholding level: global threshold, level dependent threshold • Wavelet basis function selection: 'Haar', 'Daubechies', 'Symlets', 'Coiflets'  For the wavelet packet analysis, there is one more option, that is: • Best basis algorithm: Shannon, threshold, norm, log energy, SURE, no best basis algorithm  4.1 Estimation of One-Dimensional Data In this simulation, four signals with different characteristics are considered as shown in Figure 4.1. The first signal contains bumps and high frequency disturbances of the type often found in paper machine data. The second signal, a sine wave, is estimated and compared with the third signal which contains three different frequencies. The last signal is piece wise constant. To measure the error occurred in each estimation, the mean square error (MSE) defined in Equation (3.20) is employed. White noise is added to each signal before estimation carried out.  D i f f e r e n t ID s i g n a l s  Noisy s i g n a l s  l  0 0  50  100  150  200  250  0  50  100  150  200  250  0  100  200  300  400  500  0  100  200  300  400  500  0  100  200  300  400  500  0  100  200  300  400  500  l •a  f  0.5 1:  0 -0.5  1  U: 0  100  U 200  0.5  •  300  •  400  ¥1  0  :  -0.5  •  500  0  100  200  300  400  Figure 4.1 One-dimensional signals to be estimated, (a) Bumps (b) Pure sine (c) Signal containing different frequency components (d) Blocks  35  500  CHAPTER 4 Data Analysis  Table 4.1 displays the mean square errors between the estimated signals and the original signals. The first row shows the mean square errors calculated using the raw signals before estimation. Note that the noise variation is higher for the two middle cases. The numbers shown inside the parenthesis in the Gram-polynomial case are the polynomial orders used in the estimations. As mentioned above, wavelet and wavelet packet have many option choices. For the wavelet and wavelet packet analyses, an exhaustive search of 5280 estimates was used to test nearly all possibilities of combined options. The best results are listed in Table 4.2 and Table 4.3.  Estimation Methods  Signal a  Signal b  Signal c  Signal d  No Estimation  10.3  41.3  40.2  10.3  Cirum-poKiiomkil  7.7 (118)  1.6(16)  92.4 (148)  7.5 (118)  Wave lei  5.1  2.0  12.0  2.5  Wave Id packd  6.4  2.0  9.7  2.6  Table 4.1 Best M S E of four signals using different methods {unit • 10 ) 3  Signal T\ pes  Levels  Noise Methods  Threshold Types  Threshold level  MSE  \Va\tlet Types  (io" 3  Signal a  3 •  SURE  Hard  level dependent  sym7  5.1  Signal h  5  SURE  Hard  global  db9  2.0  Signal l"  4  Universal  Soft  level dependent  sym5  12.0  Signal d  5  S U R E hybrid  Soft  level dependent  haar  2.5  MSE  Table 4.2 M S E with the best choice of wavelets  Signal Types  Levels  Signal a  Noise Methods  Best Tree Algorithm  Threshold level  \Va\elet Types  i !()-'.  2  SURE  Shannon  level dependent  sym7  6.4  Signal h  5  SURE  Threshold  global  db7  2.0  Signal c  4  SURE  Threshold  global  coif2  9.7  Signal d  5  SURE  Shannon  level dependent  haar  2.6  Table 4.3 M S E with the best choice of wavelet packets  36  CHAPTER  4 Data Analysis  Now, a comparison of the estimations using the three methods is given more in detail. In Table 4.1, it is evident that after estimation, the mean square errors for all signals have been decreased significantly except in the Gram polynomial case where the' high order polynomials are tracking the high frequency noise. The right side of Figure 4.2 shows how the M S E changes as the order of polynomial increases. For the signals a, c, and d, M S E decreases abruptly at low orders and decreases gradually at high orders. Using this method, high order Gram-polynomials, which require more computation, must be used to reduce the error closer to the level achieved by wavelet analysis. In spite of using high order Gram-polynomials, the estimation cannot obtain the M S E as low as the value resulting from the wavelet method. The Grampolynomial method works particularly poorly in estimating the signal c, consisting of sine waves of different frequencies. As expected, however, the Gram-polynomial method gives the lowest M S E for the simple sine wave because the basis function has a similar shape. The left column in Figure 4.2 shows the estimated signals for the order marked by ' X ' in the M S E vs. order plot. The dashed line indicates the original signal, and the bold line is the estimate.  Gram-polynomial E s t i m a t i o n  MSE v s . polynomial  0.08  1  0.06  II  0.04  50  100  150  200  250  0  100  200  300  400  500  1Illllllflfllgtggnflimyo.m, L.J II n nn.no n .  0.02 0  order  1  0  '•  nninnnf-finnnnnnnnn 150  I XI  o -l 150  2 0.6 u  0  0.4 0.2  -2 0  100  200  300  400 ..  0  500  0  50  100  150  0.1 0.05' 0  100  200 300 400 measurements >  0  500  nnnnnnnnnninninnn^nnnnn  0  50 100 polynomial orders  Figure 4.2 Gram-polynomial estimation.  37  150 >  CHAPTER 4 Data Analysis  Generally, the wavelet and wavelet packet methods give reasonable mean square errors for all of the signals, as shown in Table 4.1. For signals a, b, and d, both wavelet and wavelet packet work similarly and give similar results. Theoretically wavelet packet works better than wavelet analysis. However, the results from wavelet analysis for signals a and d are just slightly better in Table 4.1 since some options are fixed for the simulations using wavelet packet method. Table 4.2 and Table 4.3 shows that these methods use the same mother wavelets and similar thresholding. In estimating the signal c, the wavelet packet method is superior to the other methods since the best decomposition tree structure is chosen enabling the estimation to be more precise in regions where the frequency changes. In the wavelet estimation, however, the tree structure is fixed as shown in Figure 3.5. Figure 4.4 displays the tree structure for the wavelet packet estimation of the signal c. In general, the wavelet packet method works better in more complex situations because more detailed analysis at a particular frequency region is possible. In Donoho's paper [16], the mathematical analysis comparing efficiencies of those two methods is given.  As explained, in wavelet and wavelet packet analyses, the combination of different options results in different MSE's. To examine which option sets yield the best estimation for the wavelet method, the best 50 -70 combinations that result in the lowest MSE's are selected from among 5280 estimation trials. The number of each option usage in the 50 ~ 70 best estimations is counted and plotted in Figure 4.5. The left upper columns of 5 figures display the examples for signal a. Here, the plots show that the best choice in estimating the first signal for decomposition level, noise method, threshold type, threshold level type, and wavelet type are level 3, S U R E , Hard, Level Dependent Thresholding, and Symmlet 7 respectively. The details of different options are explained at the beginning of this section. The cases for signals b, c, and d are also shown in the bottom and left columns of Figure 4.5. The best level of decomposition varies from 2 to 5 in all signals. It is also clearly shown that the first noise method (SURE), hard thresholding, and level dependent thresholding are best for the most signals generally, allowing these options to be set at fixed value. The choice of mother wavelet is also an important factor. For blocks, the Haar mother, wavelet performs best due to its similarity in shape.  38  CHAPTER 4 Data Analysis  For the wavelet packet estimation, the results are similar and shown in Figure 4.6. The additional option, the best basis criterion, has clear preferred setting for signals a and d. In the analysis of twodimensional data, since more computational effort is needed, some options are fixed and only the significant variables are manipulated in the next section.  Wavelet  Estimation  Wavelet packet E s t i m a t i o n  0  50  100  150  200  250  0  50  100  150  200  250  0  100  200  300  400  500  0  100  200  300  400  500  0  100  200  300  400  500  0  100  200  300  400  500  1 •a  0.5 0 -0.5 .0  1  lm 100  0.5 0 200  300  400'  -0.5  500  -4 0  "in  100  200  300  Figure 4.3 Wavelet and wavelet packet estimation  (0,0)  0.0)  (1,1)  (2,0)  (3,0)  (4,0)  (2,1)  (3,1)  (4,1)  (4,2)  (3,2)  (4,3)  (3,3)  (4,6)  Figure 4.4 Wavelet packet tree structure for signal c  39  (4,7)  400  500  CHAPTER 4 Data Analysis  (a) S i g n a l a: T o t a l count = 50  Decomposition 4J  C 20  o U  (b) S i g n a l b: T o t a l count = 76  0  level  SURE  0  10  Minimax 1 1  Univ.  n  a  Minimax  Noise Method  40  •  20  •  Hard  Soft  Threshold 20  level dependent  0  <—Db4  n  <- —Db9  type  20  n .  10 15 K i n d of Wavelet '  (c) S i g n a l c: T o t a l  level dependent  THreshold l e v . |<—Sym7  n  n  5  Global  0  Threshold lev.  10  count  k — D b 9 & DblO  0  20  54  5  10 15 Kind of Wavelet  20  (d) S i g n a l d: T o t a l count = 54  20 \  4J  a 10  10  O  u Decomposition SURE  10  Univ.  0  level  SURE hybrid  0 o  Decomposition  20  SURE  u 10 • 0  Noise Method  40  0 30  nt  0 20  SURE hybri  40  O O  0  Umv.  level  0  4J  5  SURE  0  T h r e s h o l d type  50  Decomposition  20  level  . ' SURE Univ hybrid  Minimax  Noise Method 4-1  20  Hard  a 20  Soft  0 u  0 T h r e s h o l d type  0  Global 1  1  Threshold  3 0  u  level dependent^  0  4->  S 20  10 15 Kind of Wavelet  type  a  level dependent  Threshold l e v .  5  Soft  0 50  40 20  Hard "  0 0  20  Figure 4.5 Count of options used by first 50  k—Haar D a  • _  5  Threshold l e v . 1: Haar 2-10: Daubechies 11-17: Symlets 18-22: C o i f l e t s i  •—•! I —  10 15 Kind o f Wavelet  —  20  70 best results that give lowest M S E using wavelets  40  CHAPTER 4 Data Analysis  (a) S i g n a l a: T o t a l count  = 56  (b) S i g n a l b: T o t a l count  =' 60  20 10 0 Decomposition  15 10  SURE Univ.  5  level  SURE hybri  Decomposition 10  Minimax  0  0  Univ.  SURE  Noise Method 20  Shannon  10  0  Best B a s i s A l g o r i t h m  50  a  a  Global  o  level | dependent) THreshold l e v .  CJ  0 20 10  •  Threshold l e v . k'—Db4  4-1  5  10 15 Kind of Wavelet  (c) S i g n a l c: T o t a l count  -  -Db9 <—DblO  1 10  <—Db7  •  20  •  < — Sym.8 .  u  <—Sym7 0  SURE  log  Ithrtt  4-)  3 o  CJ  Minimax  no-bestj  Shannon  0  Best B a s i s A l g o r i t h m  50  SURE hybri  Noise Method  10  J->  level  n  20  5  = 57  10 15 K i n d of Wavelet  20  (d) S i g n a l d: T o t a l count §.  = 54  2 0  o 10 Decomposition  level  Decomposition 4J  % 10  SURE  o o  0  Noise Method thre.  20  log  g  SURH no-bestj  0 40  4->  20  Global  thre. 1  level dependent!  o  CJ  0 Threshold lev.  a 20 3  o  Threshold l e v . : ' 1: Haar 2-10: Daubechies 11-17: Symlets k—Haar 18-22: C o i f l e t s  CJ  5  10 15 Kind of Wavelet  20  SURH  1  Best B a s i s A l g o r i t h m  50  a  level dependent  Minimax  Noise Method  1 0  0  Best B a s i s A l g o r i t h m  SURE hybri  Shar non  a  10 SharW  Univ.  level  o  5  m  10 15 K i n d o f Wavelet  -Q. 20  Fi: ;ure 4.6 Count of options used by first 50 ~ 70best results that give lowest M S E using wavelet packets  41  CHAPTER 4 Data Analysis  4.2 Estimation of Two-Dimensional Paper Machine Data After being sampled by a sensor, the basis weight measurements are stored into a matrix. Each row of the matrix contains the data collected during one scan and each column stores the profile measurement at an actuator position. By adjusting actuator settings in the paper machine model given in Chapter 2, various types of two-dimensional data sets can be created for error analysis of the different estimation methods. In this section, three different paper machine data sets are created from the model and then they are estimated using the three different filtering techniques: exponential filter, two dimensional wavelet method, and two dimensional wavelet packet method. Then, the estimated or de-noised signal is compared with the desired signal which is the original signal without measurement noise. Similar to the onedimensional analysis, mean square error between these two signals is used to assess the efficiency of each method.  4.2.1  Smooth waves The first data set that has been generated from the paper machine model is made up with a smooth  sine wave across the sheet, a fluctuation in M D , and Gaussian noise. More specifically, during the scanning period 1 to 40, actuators are adjusted to a sine wave as shown in Figure 4.7 (a). From scan 41 onward, all actuators are set to zero. A two dimensional view of the actuator setting is depicted in Figure 4.7 (b). The basis weight profiles, as shown in Figure 4.7 (c), result from the time delay of 2-scans as well as the first order dynamic in the machine direction. For the cross direction, there are interactions between actuators, approximated by the interaction matrix G as given in equation 2.2. The M D signal in Figure 4.8 (a) is added to the previously obtained pure C D profile as shown in Figure 4.7 (c). As mentioned before, the machine direction fluctuation of basis weight results mainly from variations of consistency in the headbox, and from the pump speed changes. The M D signal is approximately periodic waves. Adding white Gaussian noise z  f  • ~ N(0,G  2  ) where a  2  = (0.2)  2 '  to the combination of M D and C D yields the  final paper machine data set that will be used in the simulation (See figure 4.8 (c)).  42  CHAPTER 4 Data Analysis  (a) Actuator set p o i n t changes  0  10  20  30  40 50 60 Actuators l o c a t i o n s  70  80  (b) Actuator adjustments i n a d i f f e r e n t view  (c) R e s u l t i n g data s e t  Actuator l o c a t i o n s  Measurements  Figure 4.7 Initial actuator settings and resulting data set (Smooth waves)  (a) MD  Figure 4.8 Generated paper machine basis weight data  43  90  100  CHAPTER 4 Data Analysis  4.2.1.1. Estimation The raw data given in Figure 4.8 (c) is processed using different wavelet and wavelet packet estimators, in addition to the E X P O method. Before estimation, the M D components are extracted from the raw data by subtracting scan averages of M D samples. Thus, the C D components combined with noise are used in the estimations. The simulations using various options in wavelet and wavelet packet estimations have yielded the results, in terms of M S E , as shown in Table 4.4 and Table 4.5. The lowest M S E is given at the top of the tables, and the result from use of the E X P O is given at the bottom. In comparison to the E X P O , it is shown that both wavelet and wavelet packet methods decrease the errors to less than 5% of that for the E X P O . However, there are no significant differences in mean square error between the wavelet and the wavelet packet methods.  For both methods, the choice of wavelets is an important factor in decreasing the estimation errors, as shown in Figure 4.9. It is not surprising that the smoother wavelets perform better for the estimation of the smooth profile. In effect, the best results using the Haar wavelet are about 5 times larger than the estimation errors using the smooth wavelets. The threshold values are automatically chosen according to the noise level of the signal. Figure 4.10 displays the clean C D profile and the estimated profiles in 3D view. The images of the estimated data, as shown in Figure 4.11, also illustrate that the wavelet and wavelet packet methods are superior to the E X P O method.  Figure 4.12 presents the cumulative estimated errors along the cross direction for each method. Here, errors along M D are added together to show these error profiles. The signal amplitude is higher in the error profile from the E X P O method. These figures also reveal that the error signals for both wavelet and wavelet packet have the saw-tooth form. This comes about because the original signal that is used for estimation has been created using the interaction matrix G that introduces high frequency (saw-tooth like) actuator responses. Error calculation that is obtained by subtracting the estimated signal from the original signal does not remove this high frequency component since the estimated signal using wavelets and  44  CHAPTER 4 Data Analysis  wavelet packets has a smooth surface. Consequently, the error signals in Figure 4.12 show this saw-tooth effect. It can also be seen that the edge effects are present in the wavelet and wavelet packet estimations:  Levels  Noise Methods  Threshold Types  Threshold level  Wavelet Types  MSK  4  SURE  Soft  Level dependent  sym6  4.3  3  SURE  Soft  Level dependent  sym6  4.6  2  SURE  Soft  Global  sym6  5.6  4  S U R E hybrid  Hard  Level dependent  coif4  6.5  3  SURE  Soft  Global  db5  5.1  4  Universal  Hard  Global  db8  9.7  4  SURE  Soft  Level dependent  haar  26.6  Exponential filtering estimation  124.5  Table 4.4 M S E for the wavelet estimation of the smooth waves  MSK  Levels  Noise Methods  Best Tree Algorithm  Threshold level  3  SURE  Shannon  Global  sym5 .  4.2  3  SURE  Threshold  Global  sym5  4.6  3  SURE  Norm  Global  sym5  6.4  3  SURE  Log energy  Global  db4  6.5  4  SURE  Shannon  Global  sym7  4.6  3  SURE  Shannon  Global  haar  19.9  Wavelet Types  Exponential filtering estimation Table 4.5 M S E for the wavelet packet estimation of the smooth waves  45  (in  \  124.5  CHAPTER 4 Data Analysis  Wavelet Estimation using Different Wavelets  Wavelet Packet Estimation using Different Wavelets  Figure 4.9 M S E for estimations using different wavelets.  46  CHAPTER 4 Data Analysis  (>  (a) Clean CD Profile  b  measurements  scans  M  S  E  0  Figure 4.10 Estimated profiles of the smooth waves  47  e  x  p  0  o = °-  1  2  4  5  measurements  CHAPTER 4 Data Analysis  Figure 4.11 Estimated profiles of the smooth waves (Images)  48  CHAPTER 4 Data Analysis  Error along the sensor measurements  150  200  250  300  150  200  250  300  150 measurements  200  250  300  100  Figure 4.12 Cumulative error along C D for the smooth wave  4.2.2 Signal Containing Several Frequency Components The paper machine data set, depicted in Figure 4.14 (c), has been created to test the performance of the different methods for estimating signal containing several frequency components. The C D and M D profiles given in Equation 4.1 and Equation 4.2 are corrupted by the white Gaussian noise z- • ~  N(0,C  2  2 . 2  ) where G = (0.2) . This M D profile is shown in Figure 4.14 (a).  CD = s i n 0  c m  + sin0  49  C £ ) 2  + sin6  C D 3  (4.1)  CHAPTER 4 Data Analysis  where 0 CDl G [0, 271 ] for the actuator positions [1, 100], 0 CDl G [0, 30TC] for the actuator positions [70, 9 0 ] , 0  C Z ) 3  G [0, 507U] for the actuator positions [20, 40]  MD = sinO MD  where 0 MD = [0,2K]  (4.2)  The corresponding actuator positions for generating the C D profile (Equation (4.1)) are shown in Figure 4.13 (a) and (b). Note that there are three different frequency components existing in the actuator setting. Two different high frequency components are added during the actuator positions from 20 to 40 and from 70 to 90. The resulting data set without noise is also given in Figure 4.13 (c).  4.2.2.1. Estimation The estimations have been processed by using the wavelet, wavelet packet, and E X P O methods and their results are presented in Table 4.6 and Table 4.7. Compared to M S E from the E X P O estimation, M S E for the wavelet estimation has been decreased to 15% and that for the wavelet packet to 11%. The wavelet packet estimation decompose a signal according to its frequency characteristics using the Best Basis Algorithm as presented in Section 3.3.2. Thus, this method analyzes the signal with different frequencies in a more efficient way. As a result, the overall MSE's for the wavelet packet estimation are lower than those found from other methods. The resulting data sets estimated from those three methods are displayed in Figure 4.15 (3D view) and Figure 4.16 (images). The cumulative error along C D , as shown in Figure 4.17, reveals that the E X P O method is not able to detect the frequency changes in the signal. Consequently, the errors at the locations where the signal has a different frequency content, are high and seem to be affected by the amplitude of the original signal. In Figure 4.17 (b), the wavelet analysis has decreased the errors significantly. Nevertheless, the wavelet packet analysis is superior to the wavelet method in detecting the frequency differences in the signal. The wavelet packet almost flattened the errors  50  CHAPTER 4 Data Analysis  1  1 i n  -TTnTtirTTTTlT^  •  (a) Actuator set point cihanges !  ir in in U u  i  i 10  0  i  20  i  1 rTTTrinrv-^ —^ULUJLUJLULIJJULI i 40  i  fl  :U  i 30  i  50  60  70  Actuators locations  L  J | I I 1 80  [JJXILILM-^  90  (b) Actuator adjustments in a different view  (c) Resulting data set  Actuator locations  Measurements  100  Figure 4.13 Initial actuator settings and resulting data set (Signal containing several frequency components)  (a) MD  0  10  20  30  40  50  Scans (b) Clean MD + CD  (c) Noisy signal to be estimated  Figure 4.14 Generated paper machine basis weight data  51  CHAPTER 4 Data Analysis  throughout the C D range, as illustrated in Figure 4.17 (c).  MSE  Noise Methods  Threshold Types  Threshold level  Wavelet Tj pes  (10 ')  2  SURE  Hard  Level dependent  coif5  9.1  3  SURE  Hard  Level dependent  coif5  11.4  4  SURE  Hard  Level dependent  coif5  11.9  2  S U R E hybrid  Soft  Global  db6  10.9  Hard  Level dependent  sym8  9.4  Soft  Level dependent  haar  18.6  Le\els  2 2  SURE • S U R E hybrid  Exponential filtering estimation  57.1  Table 4.6 M S E with the best option sets for wavelet estimation (signal containing several frequency components)  MSE  Noise Methods  Best Tree Algorithm  Threshold level  4  SURE  SURE  Global  coif4  6.7  4  S U R E hybrid  SURE  Global  coif4  6.7  3  SURE  Log energy  Global  coif4  7.2  4  SURE  Log energy  Global  db6  7.5  4  S U R E hybrid  SURE  Global  sym7  7.6  3  S U R E hybrid  Shannon  Global  haar  18.3  Levels  Exponential filtering estimation  Wa\clet Types  no  S  57.1  Table 4.7 M S E with the best option sets for wavelet packet estimation (signal containing several frequency components)  52  CHAPTER 4 Data Analysis  (a) Clean CD Profile  <) b  measurements  s c a n s  M S E  0  e  xpo =  0  0  0  5  7  1  measurements  Figure 4.15 The estimated profiles of the signal containing several frequency components  53  CHAPTER 4 Data Analysis  (a) Noisy CD  (b) EXPO estimation  measurements  measurements  (c) Wavelet estimation  (d) Wavelet Packet estimation  measurements  measurements  Figure 4.16 The estimated profiles of the signal containing several frequency components (Images).  5 4  CHAPTER 4 Data Analysis  Error along the sensor measurements  300  250  300  250  300  0.2 h  a Q_  g -0.1  100  150 measurements  200  Figure 4.17 Cumulative error along C D for the signal containing several frequency components  4.2.3  Signal with Bumps and a Streak On a paper machine, the scanned data often contains response to actuator step changes and streaks  resulting when a few adjacent actuators are manipulated. To test the estimation of data containing such bump responses and streaks, the actuator settings shown in Figure 4.18 (a) and (b) are used. The resulting C D profile (Figure 4.18 (c)) from these actuator settings is added with the M D profile (Figure 4.19 (a)) 2  together with the white Gaussian noise z^  2  2  ~ N(0, G ) where CJ = (0.2) . Using a paper machine  model, two positive bumps, two negative bumps and one streak are created at the following locations.  55  CHAPTER 4 Data Analysis  • Positive bumps: (amplitude 1) actuator number 23-26 during scan 30-40 actuator number 74-77 during scan 10-20 • Negative bumps: (amplitude-1) actuator number 23-26 during scan 10-20 actuator number 74-77 during scan 30-40 • Streak: (amplitude 0.5)  actuator number 49-51 during scan 1-60  Figure 4.19 (c) displays the resulting paper machine data used in the simulation.  4.2.3.1. Estimation Table 4.8 and Table 4.9 show results obtained from the simulation. Again, the best results are shown at the top of the tables. Since the data set in Figure 4.19 (c) has null values in most locations except those where bumps and a streak occur, the errors from the E X P O estimation are relatively low when compared with the previous results from this method. In estimating bumps and a streak, the wavelet packet estimation is slightly better than the wavelet case. Moreover, both smooth and non-smooth (Haar) wavelets accomplish estimations that give reasonable mean square errors as a result of the shapes of bumps and a streak, which appear sharp in the cross direction and smooth in the machine direction. The estimated profiles resulting from simulations are displayed in Figure 4.20 (3D view) and Figure 4.21 (images). The cumulative error along C D , as shown in Figure 4.22, shows that the E X P O method gives high estimation errors along the measurement locations, and the error is especially large at the location of the streak; that is between measurement points 147 and 153. Nevertheless, the errors from the wavelet and wavelet packet estimations are lower, i  56  CHAPTER 4 Data Analysis  (a) Actuator set point changes  as ¥  0  ->—•  o <  -1t= 10  20  30  40 50 60 Actuators locations  (b) Actuator adjustments in a different view  70  8G  90  100  (c) Resulting data set  20  I  40  60  1  20  40 60 80 Actuator locations  100  100 200 Measurements  300  Figure 4.18 Initial actuator settings and resulting data set (Signal with bumps and a streak)  Figure 4.19 Generated paper machine basis weight data  57  CHAPTER 4 Data Analysis  MSE  Noise Methods  Threshold Types  Threshold level  Wavelet Types  3  SURE  Soft  Global  db4  5.4  2  SURE  Soft  Global  db4  6.3  4  SURE  Soft  Global  db4  8.7 .  3  Minimax  Hard  Level dependent  sym5  6.1  3  Minimax  Hard  Level dependent  coif2  6.2  5  SURE  Hard  Level dependent  haar  8.7  Levels  Exponential filtering estimation  (10  \  38.7  Table 4.8 M S E with the best option sets for wavelet estimation (signal with bumps and a streak)  Noise  MSE  Methods  Best Tree Algorithm  Threshold level  SURE  Threshold  Global  db4  5.2  3  S U R E hybrid  Threshold  Global  db4  5.2  3  SURE  Shannon  Global  db4  7.0  3  SURE  Norm  Global  sym7  6.1  3  S U R E hybrid  Log energy  Global  Coif3  6.2  3  SURE  Shannon  haar  8.5  Levels  3'  Global .  Exponential filtering estimation  Wavelet Types  (10  !  38.7  Table 4.9 M S E with the best option sets for wavelet packet estimation (signal with bumps and a streak)  58  CHAPTER 4 Data Analysis  (a) Clean CD Profile  () b  measurements  s c a n s  0  MSE eX  0  po =  0  0  3  7  measurements  Figure 4.20 The estimated profiles of the signal with bumps and a streak  5 9  8  CHAPTER 4 Data Analysis  (a) Noisy CD  (b) EXPO estimation  measurements  measurements  (c) Wavelet estimation  (d) Wavelet Packet estimation  measurements  measurements  Figure 4.21 The estimated profiles of the signal with bumps and a streak (Images)  60  CHAPTER 4 Data Analysis  Error along the sensor measurements  O  Q_ X W  -0.05  150  200  250  300  200  250  300  200  250  300  0.15F  150 measurements  Figure 4.22 Cumulative error along C D for the signal with bumps and a streak  4.2.4 Industrial Data Analysis The industrial data set, presented in Figure 4.23 (b), has been obtained from a coater with several actuator manipulations during 15 scanning periods. The data represents changes in sheet basis weight as a result of the addition of coating material. The actuator set point changes are shown in Figure 4.23 (a). The specifications of this data are given as follows:  • # of high resolution weight measurements: 360 # of scans: 37  61  CHAPTER 4 Data Analysis  • # of actuators: 78 • Coat weight target: 7gsm • Four bumps during scan 10-26 resulted from following actuator settings: Positive bumps: (amplitude shown inside parenthesis) actuator number 37 (0.9), 38 (2.2), and 39 (0.9) actuator number 67 (1.9), 68 (2.2), and 69 (0.9) Negative bumps: (amplitude shown inside parenthesis) actuator number 22 (-0.9), 23 (-2.2), and 24 (-0.9) actuator number 52 (-0.9), 53 (-2.2), and 54 (-0.9)  The estimation of the noisy data has been conducted using three different methods, the E X P O , wavelet, and wavelet packet analyses. The images in Figure 4.23 illustrate that it is hard to distinguish the locations of bumps in the E X P O estimation, while the wavelet and wavelet packet methods clearly identify two positive bumps and two negative bumps. The simulation results are also shown in Figure 4.24 in 3D plots. Note that the basis weight value is normalized in the range of 0 to 1 in this figure.  As mentioned before, the wavelet packet method gives results that are similar to the wavelet method. However, the wavelet packet method is more flexible since it selects the best decomposition tree structure among many possible options. For example, the Best Basis Algorithm using five different entropy types generates the five different best tree structures as shown in Figure 4.25 (b) to (f). If the options 'threshold', 'norm', and ' S U R E ' are selected in entropy calculation, the decomposition tree structures will differ depending upon what parameter values are chosen. In Figure 4.25 (c), the threshold value 1.5 was chosen and the tree structure becomes the same as that for the wavelet method as displayed in Figure 4.25 (a). For 'norm' and ' S U R E ' , the values 1.8 and 15 are used respectively. The different tree structures result in different estimations. In this examples using the above tree structures, the difference in accuracy of the  •  62  CHAPTER 4 Data Analysis  (a) Actuator Adjustments  30  40 Actuator positions (b) R a w Data  5 10 15 20 25 30 35 50  150  200  350  measurements — > (o) E X P O Filter Estimation 5 10 15  .3 CO  20 25 30 35 150  200 measurements — > (d) Wavelet Estimation  150  200 measurements  250  >  (e) Wavelet Packet Estimation 5 10 15 20 25 30 35 150  200 measurements -  Figure 4.23 The estimated profiles of the industrial data (Images)  63  CHAPTER 4 Data Analysis  estimated data is not comparable in this case. It is clear from Figure 4.25 that the wavelet packet tree structure analyzes the high frequency region more precisely than the wavelet.  (a) Raw Data  (b) Expo Filter  400  scans  0  o  400  scans  measurements  (c) Wavelet  0  0  (d) Wavelet Packet  400  scans  0  0  measurements  400  scans  measurements  0  0  Figure 4.24 The estimated profiles of the industrial data  64  measurements  CHAPTER 4 Data Analysis  (a) Wavelet method  (b) Shannon method  (c) Threshold method: 1.5  (d) Norm method: 1.8  (e) Log energy method  (f) Sure method: 15  Figure 4.25 Decomposition trees resulted from the Best Basis Algorithm  65  CHAPTER 4 Data Analysis  4.3 Data Compression In investigating data compression, the estimated data sets obtained in the previous section are tested with the compression techniques described in Section 3.7. The data compression procedures are summarized into three steps: decomposition, thresholding of wavelet coefficients, and reconstruction. In the decomposition process, the wavelet and wavelet packet analysis transforms the data into fewer wavelet coefficients. For example, the estimated smooth wave with the size 60 by 300 is decomposed and its wavelet coefficients are displayed in Figure 4.26. Figure 4.27, shown with a different axis scale, illustrates that first 1000 coefficients out of total coefficients flunctuate in amplitude within the range ± 2 0 , and these concentrated coefficients have larger values compared to the remaining coefficients. The thresholding reduces to zero all coefficients between the dashed lines in Figure 4.26 (a). A l l coefficients except those that have relatively large values become null values, providing the compression. In this thesis, the compression ratio and the mean square error in Equation (3.20) are used to determine how well data is compressed. To achieve a high compression ratio, it is desired to have more zero coefficients and small mean square error after recovering signal from data compression.  Table 4.10 illustrates the efficiency of the data compression in terms of compression ratio and mean square error between the signals before and after the compression, using different options for the given data sets. For all cases, the wavelet packet compressions achieve higher compression ratios with lower mean square errors.  The original data sets before data compression and the recovered signals after compression are shown in Figure 4.28 in an image view. The data compression is lossy since the recovered data has lost some values in the thresholding process. Nevertheless, the recovered signals have retained most of their energy while keeping the compression ratio high so that the recovered signals in Figure 4.28 show very similar images to the original signals.  66  CHAPTER 4 Data Analysis  (a) Wavelet c o e f f i c i e n t s  f o r t h e smooth wave b e f o r e  3  4  (b) A f t e r  1.5  thresholding  thresholding  7 x 10  4  1 0.5 0 0.5 -1 _l_  1.5  7 x 10  Figure 4.26 Truncated wavelet coefficients before and after thresholding  Magnified view of Figure 4.26(a)  -10  -20 0  500 1000 1500 2000 2500  Figure 4.27 First 2500 coefficients before thresholding  67  4  CHAPTER  Estimation  Signal used  Compression  Wavelet types  methods wavelet  db3  wavelet packet  4 Data Analysis  M S E (10  ratio 38.81 : 1  9.4  coif3  47.14: 1  2.5  wavelet  sym4  19.58 : 1  15.3  wavelet packet  coif5  37.39 : 1  5.1  wavelet  sym4  37.03 : 1  5.4  wavelet packet  db4  59.23 : 1  1.5  wavelet  sym5  31.73 : 1  0.01  wavelet packet  sym5  64.00: 1  0.01  )  3  Smooth wave  Signal with different frequencies  Signal with bumps and streak  Industrial data  Table 4.10 Data compression results using wavelet and wavelet packet methods  Before Compression  100  200  After Compression: Wavelet  300  100  After Compression: Wavelet Packet  100  200  300  200  300  100  200  300  200  300  100  200  300  20 40 60 100  100  200  200  300  300  Measurements  100  100  If m  200  300  Measurements  100  Figure 4.28 Images before and after the data compression  6 8  200  300  Measurements  Chapter 5 CD Control of Basis Weight  After the estimation of the true profile using data that is measured by a sensor, actuators in a headbox are manipulated in order to achieve a uniform distribution of basis weight in the cross direction of a paper sheet. It has been shown in chapter 2 that the process model has a first-order response with a timedelay, and that the couplings between actuators may be handled by the interaction matrix G . Using these characteristics, the Dahlin controller, which is a dead time compensator, is implemented because of its simple design. Later in this chapter, a method for profile control performance assessment is also discussed.  5.1 Dahlin Controller The Dahlin algorithm can be applied to an arbitrary process, forcing the closed loop response of the feedback system to be a unit gain first-order process with dead-time [25]. The closed loop dead-time is equal to the open loop value. This algorithm has the only one tuning parameter, p, that is the closed loop pole location and should be positive and less than one. It is assumed that the process dynamic of a discrete plant, H (q  1  QL  ,-')=. "° ^;' - <v (i+  Hoi(  A{q  )  l+a,q  +  +fc  +---+a q N  . (5..)  "  M~ )  -l A controller, C(q  ) , can be described by  B  ) = \(f  l  — can be designed using feedback as shown in Figure 5.1 so that the )  closed loop system becomes a unit gain first-order process with a time delay. The controller first cancels out all the stable poles and zeros of H  QL  which means that the process has to be stable. Since the model  69  Chapter 5  CD Control of Basis Weight  used for basis weight control has a stable first-order dynamic with a delay, this controller can be used. The desired closed loop system is given as  HCL^' )  = ^ l-pq  1  x  < l  (5-2)  d  where p is a tuning parameter. From Figure 5.1, the closed loop system can be derived to be H,  _ i  HM  •C  n  ) =  C  ]  +  l  +  ° t  BB  L r  t  O L  =  —r q  c  d  AA  C  1 —n  j  = - ^d  (5-3)  [  + BB q  c  -H  L  l-pq  c  Using this equation, the controller is obtained:  c(,-y=^X—— .*kh t  ydesired '  5.  - o  ( 5  .4)  -dB A  C  H OL  Figure 5.1 Dahlin controller Assume that the plant model is given as a first-order process with delay of d, as  (5.5) 1 - aq Then, the controller becomes  C(q~l  1 - aq  =  i  -pq~  l  -(i  70  -p)q~  (5.6)  Chapter 5 CD Control of Basis Weight  For a time delay equal to one, the controller in recursive form can be written as u (t) = ±j£[e(t)-ae(t-l)]  + u (t-l)  c  (5.7)  c  and if a delay is two, the controller becomes u (t) = i^P.[e(t)-ae(t-l)) c  + pu (t-l) c  + (\-p)u {t-2) c  (5.8)  5.2 C D Profile Assessment In two-dimensional wavelet analysis, the number of scans required for the estimation should be greater than 2 to 3 filter lengths [38]. In particular, the filter length for Daubechies 4 as shown in Figure 3.3 is eight samples. Hence, the required number of scans is in the range of 16 to 24 scans. The previous data boxes that include 16 to 24 scans should be estimated and then the most recent estimates are taken in the calculation of control action.  On the basis of the current profile, the next actuator settings are determined automatically in the computer control setting. However, it is necessary for an operator to monitor those processes continuously in order to respond to the changing conditions adequately. Therefore, a method with which an operator can assess the control performance is developed in this section. The industrial settings used in this simulation are as follows: The basis weight profile data were collected on a newsprint machine for 130 scans. The paper machine is equipped with N  m  - 350 measurements, N' = 130 actuators, and 3.25 inches of the  actuator width.  The profile measurements contain various frequency components that can be divided into three categories: controllable components, uncontrollable components, and noise. If there exist controllable frequency components in the profile, control actions are required to remove those components. The uncontrollable frequency components and noise are those variations that can not be removed by manipulat-  71  Chapter 5  CD Control of Basis Weight  ing the actuators. The spatial frequency that divides the controllable and the uncontrollable regions is defined as a cut-off frequency. This frequency is determined by the actuator width and the single actuator response for a particular machine. It was found that the spatial cut-off frequency should be less than 0.5 or that the wavelength be greater than 2 [40]. Here, the unit wavelength is denned to be an actuator width. The frequency is inversely proportional to the wavelength. In Figure 5.2, the dotted lines display actuator widths and the solid lines display the sine wave for the cut-off frequencies which are identified in the spectrum plot. The sine wave is plotted only up to 40 measurement points. In plotting these figures, the industrial settings given at the beginning were employed. The top two plots show the cut-off frequency 0.38 Hz or the wavelength 2.7 when the single actuator response is broader than the given actuator width. Also, two plots in the bottom illustrate the sine wave when the spatial cut-off frequency is 0.5 or the wavelength 2. This is the case when the cut-off frequency is equal to one actuator width. In the simulation, the wavelength value in the later case is used to calculate a cut-off frequency.  Input Signal: Shown in the range of first few actuator spaces  Spectrum of the left signal  Wavelength = 2.7 Frequency = 0.38  0  5  10  15  20  25  30  35  40  0  0.2  0.4  0.6  0.8  1  1.2  Wavelength = 2.0 Frequency = 0.50  0  5  10  15 20 25 30 Measurements (up to 350)  35  40  0  0.2  0.4 0.6 0.8 spatial frequency  1 >  Figure 5.2 Different cut-off frequencies or wavelength for a given actuator width.  72  1.2  Chapter 5  5.2.1  CD Control of Basis Weight  Performance Index In the wavelet and the wavelet packet analyses, some features such as their abilities to decompose  a signal into different frequencies and orthonormal characteristics of the decomposed coefficients, can be 2  used to build the performance index that measures the control performance. Let (5  be the noise-free  pr  f  2  overall profile variance and O  o n t  be the variance of the controllable component. Then, the control perfor-  mance index, C, is denned [40] as  2 C  =  V 2  'pr  G  2  < - > 5  9  ®cont  If C = 1, it denotes that the control action has removed all controllable components. C > 1 indicates that the controllable variations remain in the profile. C < 1 is not possible.  5.2.1.1 Frequency Characteristics of Wavelet and Wavelet Packet To implement the performance index in Equation (5.9), it is crucial to divide the profile data into the controllable and uncontrollable components using the given cut-off frequency. For an example, the profile from the data specified at the beginning of Section 5.2 is decomposed down to level 3 using wavelet analysis. The profile is shown at top of Figure 5.3. Also the approximation and details with their spectrums in the frequency domain are displayed in this figure. The spectra at the right hand side show that the wavelet decomposition divides the data into different frequency regions. The dashed line shows the expected ranges of frequency that is resulted from the decomposition. Moreover, the wavelength ranges of the signals are displayed between the signal and the spectrum plots. Similarly, the wavelet packet decomposition using the full tree structure, as shown in Figure 3.8 (a), results in Figure 5.4. It is clear that the wavelet packet analysis divides the signal into smaller frequency components. Hence, the wavelet  73  Chapter 5  CD Control of Basis Weight  packet provides the more accurate value in calculating the performance index. For example, the performance index calculation using the wavelet analysis gives C = 4.94 while that for the wavelet packet analysis yields C = 7.81. The lower value of C, 4.94, in the wavelet case implies better control performance than the value 7.81. However, this lower value resulted from less accurate frequency division. More specifically, using the spatial cut-off frequency 0.5, the controllable components in the wavelet decomposition are located at (3,0) and (3,1) as shown in Figure 5.3. Here, some controllable components in (2,1) are not included in the calculation. However, the controllable variations from the wavelet packet analysis includes the components at (3,2) in the addition to those at (3,0) and (3,1). It is clear that the wavelet packet yields a more precise calculation.  74  Chapter 5  CD Control of Basis Weight  Chapter 5  CD Control of Basis Weight  Chapter 5  5.3  CD Control of Basis Weight  C D Basis Weight Control The control strategy for C D basis weight control using wavelets or wavelet packets is shown in  block diagram in Figure 5.5. The inverse of the interaction matrix, G, is used to estimate present actuator settings before the estimated output proceeds to the controller, C. The parameter values for G as found in Chapter 2 are used. The collected profile data by a sensor is de-noised using wavelet and wavelet packet methods that were presented in Chapter 3 and 4. For assessment of control performance, the performance index defined in Section 5.2 is calculated with the help of wavelet and wavelet packet decompositions.  e ydes  —>o  -d  c  U  G  W  c  Dahlin Controller  W  G  q  B A  y w  Wavelet Analysis  ->  y  Paper Machine  Controllability Assessment Figure 5.5 Block diagram of basis weight control scheme using wavelet analysis  In control simulations, it is assumed that the control action takes place at every scan even though industrial practice is to take action every 2 to 3 scans. In [25], it is explained that the results from those two cases have small differences. To simplify the estimation process in this simulation, the estimator analyzes one-dimensional data in every scan instead of estimating the two-dimensional data as illustrated in Chapter 4.  In the simulation of C D control, the paper machine data is first generated. The actuators are fixed as the shape shown in Figure 4.7 (a) for 9 scan periods to generate the basis weight data. From the scan 9, the controller is turned on to calculate appropriate actuator adjustments. The simulation result is displayed in Figure 5.6. The basis weight reaches its peak at scan 9 and then gradually decreases toward zero after  77  Chapter 5  CD Control of Basis Weight  the controller is turned on. Figure 5.7 (a) and (b) show the calculation of the performance index based on wavelet analysis and on wavelet packet analysis respectively. The final performance index value reaches around to 1.9 if wavelet analysis is used and around to 2.4 for wavelet packet. It seems that the wavelet analysis yields more desirable control performance. However, as mentioned before, the performance index calculation based on the wavelet packet analysis yields more accurate results due to its ability to decompose the signal in finer frequency resolution.  As a result, the control actions of the designed controller have reduced the basis weight variation significantly and the value of performance index coefficient has been decreased dramatically after a controller has been turned on at scan 9.  Figure 5.6 Control o f C D Profile: control action took place from scan 9  7<S'  Chapter 5  Performance Index using (a) Wavelet (b) Wavelet Packet  Figure 5.7 Performance index during C D control  79  CD Control of Basis Weight  CHAPTER 6 Conclusion 6.1 Summary The thesis discusses a new method for paper machine signal analysis using the wavelet packet transformation. It has been shown that this technique can be applied to three main areas in paper machine control: paper machine data analysis including estimation and compression, control of basis weight or moisture variations, and control performance monitoring.  The paper machine model for use in C D control has been derived from experimental bump test data. This model gives a relationship between the actuator settings and the C D basis weight response. It has been shown that the model consists of actuator dynamic, time-delay, and interaction. Coupling between actuator spaces has been approximated using an interaction matrix.  To improve our knowledge of process variation from noisy sensor measurements, there exist several traditional methods; Expo-filter and Gram-polynomial approximation, and more recently wavelet estimation. A new approach, discrete wavelet packet transform (DWPT), introduced in this thesis, offers greater flexibility because it allows time-frequency resolution to change in a signal decomposition where other methods use fixed time-frequency resolution. The Best Basis Algorithm has been introduced to obtain a best decomposition tree structure in decomposition of a given signal.  It has been shown clearly that the new method, along with the wavelet estimator, produces results superior to standard analysis methods by reducing the error of estimation. In estimating a more complex signal, however, the wavelet packet method is superior to the wavelet analysis since the best decomposition tree structure is chosen enabling the estimation to be more precise in the regions where frequency changes. From the exhaustive search in the estimation of scanned paper machine data, it has been found that options  80  CHAPTER 6  Conclusion  in wavelet and wavelet packet analyses can be fixed depending on the characteristics of a given signal. When applied to industrial data, the wavelet packet estimation has provided more visual results by distinguishing the bump locations more evidently. Moreover, it has been shown that wavelet packet compression is more efficient than the wavelet case since it has achieved a higher compression ratio by yielding the combination of more zero coefficients and lower mean square error.  The control scheme has been developed using Dahlin controller. Also the control performance assessment tool developed by using wavelet and wavelet packet has been introduced. Wavelet packet provide more precise values in calculating the performance index because this method can divide the frequency range of a signal with more finer resolution. It has been shown that the control actions of the designed controller have reduced the basis weight variation significantly and the value of performance index coefficient has been decreased dramatically after a controller has been turned on.  6.2 Future Work There are some possible areas for further research. Some topics are given as follows:  •  In this thesis, only white noise was used as a noise model, but the paper machine measurements can be corrupted by various types of noise (i.e. colored noise). It is necessary to analyze the effect of estimation and compression carried out using different noise models.  •  On the assumption that the basis weight variations include time-varying components, the timevarying wavelet packet estimator can be adapted, which would achieve better estimation with reduced error. This method can also be applied to data compression.  •  More detailed analysis can be carried out in boundary handling problems using wavelet and wavelet packet methods.  81  CHAPTER 6  Conclusion  Data compression based on the wavelet technique can be further developed in paper machine data compression. For example, many researchers reported wavelet-based data compression using vector quantization, adaptive transforms, tree encodings, and edge-based coding.  82  References [I]  M . Antonni, M . Barlaud, P. Mathieu, and I. Daubechies. Image coding using wavelet transfrom. I E E E Trans Image Processing, l(2):205-220, 1992.  [2]  B . R. Bakshi and G . Stephanopoulos. Compression of Chemical Process Data by Functional Approximation and Feature Extraction. AIChE Journal, volumn 42, N o . 2, Feb. 1996.  [3]  I. Buckeit and D . Donoho. WaveLab and Reproducible Reaserch. To Appear, Wavelets and Statistics, Anestis Antoniadis, ed. Springer-Verlag Lecture Notes, 1995.  [4]  A . Bruce, D . Donoho, H . Y. Gao. Wavelet Analysis. I E E E Spectrum. Oct. 1996.  [5]  G . Burkhard and P. E . Wrist. The evaluation of paper machine stock systems by basis weight analysis. Tappi, 37(12), Dec. 1954  [6]  S. C . Chen. Kalman filtering applied to sheet measurement. In American Conference, volume 1 of 3, pages 643-647, Atlanta, Georgia, June 15-17 1988.  [7]  R. R. Coifman, Y. Meyer, S. Quake and M . V. Wickerhauser. Signal processing and compression with wavelet packets. Process in Wavelet Analysis and Applications, Meyer and Roques eds., Editions Frontiers, Toulouse, France, pages 77-93, 1992.  [8]  R. R. Coifman and M . V. Wickerhauser. Best-adapted wave packet bases with exponential decay, preprint, Yale University. 1990.  [9]  R. R. Coifman and M . V. Wickerhauser. Entropy-Based Algorithms for Best Basis Selection. IEEE Trans, on Inform. Theory, 38(2), 713, 1992.  [10]  K . Cutshall. Cross-direction control. In B . A . Thorp and M . J. Kocurek, editors, Paper Machine Operations, volume 7 of Pulp and Paper Manufacture, chapter X V I I I , pages 472-506. Tappi and C P P A , 1991.  [II]  E . Dahlin. Computational methods of a dedicated computer system for measurement and control on paper machines. In 24th Engineering Conference, TAPPI, pages 62.1-62.42, San Francisco, California, Sep 1969.  83  Control  References  [12]  I. Daubechies. Orthogonal Basis of Compactly Supported Wavelets. Comm. Pure & Applied Math. 41, pages 909-996, 1988.  [13]  I. Daubechies, A . Grossmann and Y. Meyer. Painless non-orthogonal expansions. J. Math. Phys., 27, pages 293-309, 1986.  [14]  I. Daubechies, S. Mallat and A . S. Willsky eds. Special Issue on Wavelet Transforms and Multiresolution Signal Analysis, IEEE Trans. Information Theory, 38, 2, 1992.  [15] "  D . L . Donoho. Nonlinear Wavelet Methods for Recovery of Signals, Densities, and Spectra from Indirect and Noisy Data. Proceeding of Symposia in Applied Mathematics. volumn 00, 1993. -  [16]  D . L . Donoho, and I. M . Johnstone. Ideal Denoising in an Orthonormal Basis Chosen from a Library of Bases. Technical Report 461 ,Standford University, Department of Statistics, Sept. 1994. .  [17]  D . L . Donoho and I. M . Johnstone. M i n i m a x estimation via wavelet shrinkage. Technical Report 402, Standford University, Department of Statistics, July 1992.  [18]  G . A . Dumont, M . S. Davies, K . Natarajan, C . Lindeborg, and E . M . Heaven. Estimation of Moisture Variation on Paper Machines. IEEE Transaction on Control System Technology, 1(2): 101-113, June 1993.  [19]  K . Dutton and J. F. Barrett, and M . J. Grimble. The shape control problem for a Sendzimir cold steel rolling m i l l . In Symposium on application of multivariable systems theory, pages 217-227, Oct. 1982.  [20]  G . Goodwin and K . Sin. Adaptive Filtering, Prediction and Control. Englewood Cliffs, N J : Prentice-Hall, 1984.  [21]  A . Grossmann and J. Morlet. Decomposition of hardy functions into square integrable wavelets of constant shape. SIAM J. Math. Anal, 15 pages 723-736. 1984.  [22]  C . Herley, J. Kovacevic, K . Ramchandran and M . Vetterli. Tilings of the Time-Frequency Plane: Construction of Arbitrary Orthogonal Bases and Fast Tiling Algorithm. IEEE Trans, on Signal Processing, volumn 41, N o . 12, Dec. 1993.  [23]  M . Hilton, B . Jawerth, and A . Sengupta. Compressing Still and M o v i n g Images with Wavelets. Multimedia Systems, 3(2), 1994.  84  References  [24]  B . Jawerth and W . Sweldens. A n overview of wavelet based multiresolution analyses. Technical Report 1993:1, Industrial Mathematics Initiative, department of Mathematics, University of South Carolina, 1993. Submitted to S I A M Review.  [25]  K . Kristinsson. Cross Directional Control of Basis Weight On Paper Machines using Gram Polynomials. Ph. D Thesis. U B C 1994  [26]  D . L . Laughlin, M . Morari, and Richard D . Braatz. Robust Performance o f Crossdirectional Basis-weight Control i n Paper Machines. In American Control Conference, volume 3, pages 2122-2128 June 21-23 1989.  [27]  M . M i s i t i , Y. M i s i t i , G . Oppenheim, and J. M . Poggi. Wavelet Toolbox User's Guide. The Math Works Inc. March 1996.  [28]  S. T. Morgan. Estimation and Identification for Machine Direction Control of Basis Weight and Moisture. Master's Thesis. U B C 1994.  [29]  J. G . Morlet, G . Arens, I. Fourgeau and D . Giard. Wave Propagation and sampling theory. Geophysics, 47, pages 203-236. 1982.  [30]  R. L . Mortard and B . Joseph. Wavelet Applications in Chemical Engineering. Kluwer Academic Publishers, 1994.  [31]  K.Natarajan, G . A . Dumont, and Michael S. Davies. A n algorithm for estimating cross and machine direction moisture profiles for paper machinesin Proc. 8th IFAC/IFORS Symp. Identification Syst. Parameter Estimation, Beijing, PRC, A u g . 1988, pp 1091-1096.  [32]  R. T. Ogden, Essential Wavelets for Statistical Applications and Data Analysis, Birkhauser, Boston, 1997, 206p.  [33]  O. R i o u l and M . Vetterli. Wavelet and Signal Processing. I E E E Signal Processing Magazine, 8(4): 14-38, Oct. 1991.  [34]  S. Skogstad, M . Morari, and J. C . Doyle. Roubust control o f ill-conditioned plants: Highpurity distillation. IEEE Transaction on automatic control, 33(12): 1092-1105, Dec. 1988.  [35]  G . Strang and T. Nguyen. Wavelets and Filter Banks. Wellesley-Cambridge Press. 1996  [36]  B . F. Taylor. Optimal separation of machine-direction and cross-direction variation. Tappi Journal, Feb. 1991.  85  product  References  [37]  X . G . Wang, Guy A . Dumont, and Michael S. Davies. Estimation in Paper Machine Control. IEEE Control System, pages 34-43, A u g . 1993.  [38]  N . Zoran. Paper Machine Data Analysis using wavelets. Master's Thesis. U B C 1996.  [39]  N . Zoran, M . S. Davies, G . Dumont, Paper Machine Data Analysis Using Wavelets. TAPPI Journal (to appear)  [40]  N . Zoran, M . S. Davies, G . Dumont, and D . Brewster. CD Control Diagnostics Wavelet Toolbox. Pulp and Paper Center, U B C . 1997.  [41]  M . V. Wickenhauser. Lectures on Wavelet Packet Algorithm. Washington University, Nov. 1991. •  86  Using a  Department of Math. • „  Appendix A. Gram Polynomial Estimation By definition, Gram polynomials are the discrete equidistance orthogonal polynomials with equal weighting on every data points. This technique is used to approximate the noisy measurements. The general form of the polynomial is described as  P  ( i  m\ >  r  l  \ _  V ( - 1 ) (m + j)  r  ~ mN ZJ  2  c  y = 0  i  (EQ 1)  F7i  (j\)\N-l)  ij)  The polynomial order is given as m and i is defined over the interval [0, N - 1 ] . When c  N  = 1 , the Gram polynomials in a recursive form are written as  p  m m K  =  '  (N-\){2m-\)f m(N-m)  {  2i \ N-lJ  ~  m  1  (m-l)(N-l + m) m(N-m)  (EQ 2)  a n d s e t P ( i ) = 1 and/?_,(*') = 0 . 0  The m orthogonal polynomial evaluated at N points can be constructed and it is defined as  P (0)  P,(0)  PJO)  2  PJO)  ^i(l)  P.iN-l)  P (N-l) 2  -  (EQ 3)  PJN-l)  Now, let the profile y(t) be a sequence measured at / = 0, ..., N - 1. The profile then can be T  expressed in.terms of Gram-polynomial as y = X y Q  c  + e where y  c  is the polynomial coefficients and e  the estimation residue that is consists of terns of order m + 1 and higher. The coefficients y is obtained  87  Appendix A. Gram Polynomial Estimation  from the least square estimation, that is  (EQ 4)  88  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065014/manifest

Comment

Related Items