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
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Estimation and control of paper machine variables using...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Estimation and control of paper machine variables using wavelet packet analysis Chun, John Byung-Kyu 1998
pdf
Page Metadata
Item Metadata
Title | Estimation and control of paper machine variables using wavelet packet analysis |
Creator |
Chun, John Byung-Kyu |
Date Issued | 1998 |
Description | The 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 analyzing signals with suitable time-frequency resolution. The Best Basis Algorithm has been introduced to obtain a best decomposition tree structure in decomposition of a given signal. The 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 yielding 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. |
Extent | 12514057 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-04-30 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0065014 |
URI | http://hdl.handle.net/2429/7769 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1998-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
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
Cite
Citation Scheme:
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>
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