Fourier Transforms and Wavelet Analysis: Concepts and Examples

FROM FOURIER TRANSFORMS TO WAVELET ANALYSIS: MATHEMATICAL CONCEPTS AND EXAMPLES MAY 15, 2006 Abstract. This paper studies two data analytic methods: Fourier transforms and wavelets. Fourier transforms approximate a function by decomposing it into sums of sinusoidal functions, while wavelet analysis makes use of mother wavelets. Both methods are capable of detecting dominant frequencies in the signals; however, wavelets are more efficient in dealing with time-frequency analysis. Due to the limited scope of this paper, only Fast Fourier Transform (FFT) and three families of wavelets are examined: Haar wavelet, DaubJ, and CoifI wavelets. Examples for both methods work on one dimensional data sets such as sound signals. Some simple wavelet applications include compression of audio signals, denoising problems, and multiresolution analysis. These applications are given for comparison purposes between Fourier and wavelet analysis, as well as among wavelet families. Although wavelets are a recent discovery, their efficacy has been acknowledged in a host of fields, both theoretical and practical. Their applications can be expanded to two or higher dimensional data sets. Although they are omitted in this paper, more information is available in [7] or many other books on wavelet applications. 1. Introduction Wavelets are a recent discovery in mathematics; however, their rapid development and wide range of applications make them more powerful than many other long-existing analytical tools. Conventionally, wavelets are often compared to the Fourier transform to promote their advantages. This paper will take a similar approach in attempt to illustrate wavelet transform in various applications. The Fourier transform makes use of Fourier series, named in honor of Joseph Fourier (1768-1830), who proposed to represent functions as an infinite sum of sinusoidal functions [1]. Joseph Fourier was the first to use such series to study heat equations. After him, many mathematicians such as Euler, d’Alembert, and Daniel Bernoulli continued to investigate and develop Fourier analysis [1]. From the original series, various Fourier transforms were derived: the continuous Fourier transform, discrete Fourier transform, fast Fourier transform, short-time Fourier transform, etc... Fourier analysis is adopted in many scientific applications, especially in dealing with signal processing. As the applications grew more complex over time, the Fourier transform started to reveal its inefficiencies when working with time series or data with certain characteristics. Despite the attempt to tailor the method to different groups of data, the Fourier transform remained inadequate. Consequently, wavelets received more attention as they proved able to overcome the difficulties. Date: May 15, 2006. 1 2 The first block of wavelet theory was started by Alfred Haar in the early 20th century [2]. Other important contributors include Goupillaud, Grossman, Morlet, Daubechies, Mallat and Delprat. Their different focuses helped to enrich the wavelet families and widen the range of wavelet applications. Because of the similarities, wavelet analysis is applicable in all the fields where Fourier transform was initially adopted. It is especially useful in image processing, data compression, heart-rate analysis, climatology, speech recognition, and computer graphics. This paper focuses on only a few aspects of each analysis. The first section discusses Fourier series in different representations: sinusoidal functions and complex exponential. A discretization method is also introduced so as to provide support for the discussion of Fast Fourier Transform (FFT). After illustrating Fourier analysis with concrete examples, the paper will turn to the Fourier transform’s shortcomings, which give rise to wavelets . The second section discusses three families of wavelets: the Haar wavelets, Daubechies wavelets, and Coiflets. Concepts and general mechanisms will be provided in detail for Haar wavelets and omitted for the others. Finally, we will look at the advantages of wavelets over Fourier transform through a number of examples. The paper uses three main references: Course notes in Modeling II, A Primer on Wavelets and their Scientific Applications by James Walker, and A First Course in Wavelets with Fourier Analysis by Boggess and Narcowich. As the paper is aimed at readers at undergraduate level, mathematical background of linear algebra and basic calculus is assumed. 2. Fourier Analysis Fourier analysis, which is useful in many scientific applications, makes use of Fourier series in dealing with data sets. In this section, a few representations of Fourier series and related concepts will be introduced. Consequently, several examples will implement these defined concepts to illustrate the idea of Fourier analysis. 2.1. Fourier Series. 2.1.1. Sine and Cosine Representation. A Fourier series is the expression of any function f (x) : R → R as an infinite sum of sine and cosine functions: (1) f (x) = a0 + ∞ X k=1 ak sin(kx) + ∞ X bm cos(mx). m=1 In this section, we work with continuous functions on [0, 2π], and thus, it is necessary to familiarize ourselves with the vector space of such functions. It is an infinite dimensional vector space, denoted C[0, 2π], where each point on the continuous interval [0, 2π] represents a dimension. Then, an orthogonal basis of C[0, 2π] is: {1, sin(kx), cos(mx)|k, m = 1, 2, 3, ...}. Fourier series can also be employed to write any continuous function f (x) : C[0, 2π] → C[0, 2π]. Concepts such as inner product, norm, and distance between two functions in an infinite dimensional vector space are defined in a similar manner to that of a finite dimensional vector space. However, the infinite sum gives rise to the use of an integral in the definition: WAVELETS 3 - Inner Product: The inner product of two continuous functions f (x), g(x) ∈ C[0, 2π] is Z 2π < f, g >= f (x)g(x)dx. 0 - Norm: Recall that the norm of a vector in finite dimensional vector space is calculated as k ~x k2 =< ~x, ~x >. Norm is defined similarly for a function in infinite dimensional space C[0, 2π]: Z 2π 2 k f k =< f, f >= f 2 (x)dx. 0 - Angle between two functions: The inner product of f and g can also be written as < f, g >=k f kk g k cos(θ). The angle between two functions, θ, can be achieved from the above equation. - Distance Distance between two continuous functions f and g is defined as k f − g k. The inner product formula can be applied to prove that the basis mentioned above is indeed orthogonal. Let usR consider six possible R 2π R 2π inner products of the basis 2π vector functions: 0 sin(kx)dx, 0 cos(mx)dx, 0 sin(kx) sin(mx)dx (k 6= m), R 2π R 2π 0 sin(kx) cos(mx)dx (k 6= m), 0 cos(kx) cos(mx)dx (k 6= m). The following trignometric identities [3] are helpful in proving that all six inner products are equal to 0, implying an orthogonal basis: 1 sin(x) cos(y) = [sin(x + y) + sin(x − y)], 2 1 cos(x) cos(y) = [cos(x + y) + cos(x − y)], 2 1 sin(x) sin(y) = [cos(x − y) − cos(x + y)]. 2 2.1.2. Function Coefficients. Coefficients of f can be obtained by projecting f on the corresponding basis vectors (similar to that of finite dimensional space): Z 2π < f, 1 > 1 f (x)dx, a0 = = < 1, 1 > 2π 0 Z 1 2π < f, sin(kx) > = ak = sin(kx)f (x)dx, < sin(kx), sin(kx) > π 0 Z 2π 1 < f, cos(mx) > = bm = cos(mx)f (x)dx. < cos(mx), cos(mx) > π 0 2.1.3. Complex Form. Instead of taking the integral of individual sine and cosine functions, the complex form representation of Fourier series enables us to compute the various inner product integrals simultaneously. The transformation makes use of Euler’s formula [3] eiθ = cos(θ) + i sin(θ). Simple algebra yields bk − iak ikx bk + iak −ikx e + e . ak sin(kx) + bk cos(kx) = 2 2 4 k k Thus, if we let c0 = a0 /2, ck = bk −ia , and c−k = bk +ia , then the complex form 2 2 representation of the Fourier series is ∞ X f (x) = ck eikx . k=−∞ Then differentiation and integration can be calculated simultaneously [4]: d iθ e = ieiθ = − sin(θ) + i cos(θ), dθ Z 1 eiθ dθ = eiθ = sin(θ) − i cos(θ). i To illustrate the convenience of Euler’s formula, let us compute < x, sin(3x) > and < x, cos(3x) > using the complex exponential. Note that Z 2π Z 2π Z 2π < x, e3ix > = xe3ix dx = x cos(3x)dx + i x sin(3x)dx 0 0 0 = < x, cos(3x) > +i < x, sin(3x) > . On the other hand, integration by parts yields Z 2π Z x  1 xe3ix dx = e3ix − e3ix dx 3i 3i 0 x  2π 1 = e3ix + e3ix 3i 9 0 2π 6πi 1 6πi 1 = e + e − . 3i 9 9 Since e6πi = cos(6π) + i sin(6π) = 1, it follows that Z 2π 2π 2π xe3ix dx = = −i . 3i 3 0 Therefore, 2π 0 < x, cos(3x) > = 0, 2π . 3 Integrating two functions separately would have doubled the work. < x, sin(3x) > = − 2.1.4. Fourier Transformation. Fourier transformation is a method frequently used in signal processing. As the name suggests, it transforms a set of data into a Fourier series. Due to limitation, this paper will only introduce Fast Fourier Transform (FFT). Discretization. Most data sets are available as a collection of data points. Therefore, the assumption that the function is continuous is no longer required. However, a continuous function can always be represented by N points for some positive integer N . A process in which a continuous function is translated into N representative points is called discretization. Discretization translates an infinite dimensional vector space into a finite one [4]. To discretize a function, we choose a finite number of points that represent the function from which the function itself can be sketched. Consider functions y1 = cos( π4 t) and y1 = cos( 7π 4 t) with 0 ≤ t ≤ 64. In Figure 1, these functions WAVELETS 5 a)∆t = 0.1 b) ∆t = 1 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 −0.2 −0.2 −0.4 −0.4 −0.6 −0.6 −0.8 −1 −0.8 y1 = cos((pi/4)*k) y2 = cos(7pi/4)*k 0 1 2 3 4 5 6 7 8 9 −1 10 y1 = cos((pi/4)*k) y2 = cos(7pi/4)*k 0 1 2 3 4 5 6 7 8 9 10 Figure 1. Functions y1 = cos( π4 t) and y2 = cos( 7π 4 t) with different ∆t are plotted with different increments of t (∆t) of 1 and 0.1. When ∆t = 0.1, two functions are different, but when ∆t = 1, the values plotted are essentially the same. This suggests that ∆t = 1 is too large to differentiate these two functions. The same pattern is observed for any pair of sine or cosine functions that are conjugate of each other (for example, cos(αt) and cos(βt), where α + β = 2π). Formulas When a continuous function is discretized into N points for FFT, N is usually a power of 2 for computational reasons. Suppose that N points {f0 , f1 , ..., fN −1 } are spaced out evenly on the interval [0, 2π], then the k-th point in Fourier series form is (2) Fk = N X fn e− 2πi N ·(k−1)(n−1) n=1 for 1 ≤ k ≤ N . Recall that for two vectors x, y in complex plane, their inner product is defined as < x, y >= xT ȳ. Thus, equation (2) is in fact an inner product between F = [f1 , f2 , ..., fN ]T and 2πi 2πi 2πi 2πi [e− N ·0·(k−1) , e− N ·(k−1) , e− N ·2(k−1) , ..., e− N ·(k−1)(N −1) ] for any integer k in interval [1, N ]. Hence, we can represent equation (2) by a product of F and an N ×N matrix C, whose entries are determined as follows: Ckj = e 2πi N ·(k−1)(j−1) . This process is called Fourier transformation. C is the transformation matrix. Instead of dealing with integrals, we now face easier matrix operations:      where    A=  0 0 .. . 0 1 F1 F2 .. . FN     2πi   =e N A×   0 2 ... ... f1 f2 .. . fN    ,  0 N −1 .. . 0 N − 1 2(N − 1) . . . (N − 1)(N − 1)      6 The above cross product suggests that the inverse transformation process is also feasible: N 1 X 2πi (3) fn = Fk e N ·(k−1)(n−1) N k=1 for 1 ≤ n ≤ N . Expansion to function f on the interval [α, β] The process of Fourier transformation can be extended to any function on the interval [α, β] with the appropriate discretization process. Discretize the function f into N equal subintervals, so that each subinterval is β−α (4) ∆x = or β − α = N · ∆x. N and we have N points {x1 , x2 , ...xN } on the interval [α, β]. To perform FFT, we can find an one-to-one correspondence from the function f to [0, 2π]: 2π xj ∈ [α, β] → tj = · (xj − α) ∈ [0, 2π]. β−α The inverse process is: β−α tj ∈ [0, 2π] → xj = α + · tj ∈ [α, β]. 2π Substitute equation (4) in, we get: 2π N · ∆x tj = · (xj − α) or xj = α + · tj . N · ∆x 2π The above equation suggests that discretizing f (x) over [α, β] is the same as f ( N2π ·∆x (xj − α)) over [0, 2π]. Therefore, we can use the function f (x) = a0 + ∞ X ak cos(k k=1 ∞ X 2π 2π (x − α)) + bm sin(m (x − α)) N · ∆x N · ∆x m=1 for Fourier transformation. To summarize, if we have N points with spacing ∆x, labeled as {f0 , f1 , . . . , fN } on the interval [α, β], then the relationship between function f and its FFT is expressed in the following equations: N/2 (5) f n = a0 + X ak cos(k k=1 where a0 ak bm N/2 X 2π 2π xn ) + bm sin(m xn ), N · ∆x N · ∆x m=1 F1 , N 2 = · Real(Fk+1 ), N 2 = − · Imag(Fm+1 ). N = Let ∆x is the distance between two consecutive points. The periodic nature of sines and cosines entails the period length to be at least 2∆x to yield good enough approximations (See Figure 2.1.4. Fix the period to be 2∆x. If the distance between two points is less than half the period, the intersection with the x-axis cannot be WAVELETS 7 Figure 2. ∆x and the period length detected, and hence imprecise estimation is possible. A similar situation occurs when the distance between two points is larger than half the period.) Thus, it is best to have ∆x as the distance between two points, or 2∆x as the period length. For this reason, k and m range from 1 to N/2 in equation (5). From what we know about trignometric functions sin(At) and cos(At), we have: - The functions’ periods are Period = N ∆x 2π = . A k - The functions’ frequencies are Frequency = A k = . 2π N ∆x As k, m = 1, 2, . . . , N2 , the periods measured in time unit per cycle are N · ∆x, N · ∆x, . . . , 2 · ∆x. 2 and the frequencies measured in cycles per time unit are 1 2 1 1 1 1 · , · ,..., · . N ∆x N ∆x 2 ∆x As k changes from 1 to N2 , some functions sin(ktn ) and cos(ktn ) might contribute more to the original function fn than others. To measure each function’s “importance”, we define a new concept, frequency content of k, which is given by q freq(k) = a2k + b2k . A plot of k against the frequency content yields the power spectrum of a signal. 2.2. Applications. In this section, we will examine how FFT is applied in simple examples. 8 4500 4000 3500 Frequency 3000 2500 2000 1500 1000 500 0 0 50 k 100 150 Figure 3. Power spectrum of the sunspot data. 2.2.1. Sunspot Example. Data from sunspots have interested scientists for hundreds of years, as they are the indicator for radiation from the sun and therefore have an effect on many scientific fields [5]. MatLab has a collection of data recording the number of sunspots counted on the sun every year between 1700 and 1987. Using Fourier transformation, we will analyze the data to find periodicity in the number of sunspots. Let the domain be 288 years from 1700 to 1987 (thus N = 288, ∆x = 1). So k will go from 1 to 143 (= N2 ). After the data are loaded, the frequency content of each k is plotted against k (see Figure 3). The plot peaks at the 26th position, which suggests that the component functions with k = 26 contribute the most to the FFT. Hence, the period of the sunspot cycle can be estimated using k = 26, which yields N · ∆x 288 · 1 = ≈ 11.08 years. k 26 Figure 4 confirms this finding as the frequency content peaks approximately every 11 years. To see how well the FFT estimates the real function, we will produce a graph of the original function and the approximate function by using the inverse FFT of the ten most significant component functions (with the largest frequency contents) (see Figure 5). Although only ten component functions were used, the approximation shows close estimation to the original function. Note that the imaginary part resulting from the inverse FFT is negligible since we started with a real function, and the imaginary part therefore should be insignificant. 2.2.2. Interpolation. Fourier transformation can be useful for interpolation of data. Given a set of data points, we can fit them with a continuous function using FFT coefficients. The last part of the sunspot example illustrates this method. First, WAVELETS 9 4000 3500 Frequency content 3000 2500 2000 1500 1000 500 0 2 4 6 8 10 12 Time period 14 16 18 20 Figure 4. Time cycle and frequency content of sunspot data. 200 Number of sunspots 150 100 50 0 Original function Approximate function −50 0 50 100 150 Time 200 250 300 Figure 5. Original function and the approximate function using the ten most significant component functions. Fourier transform is employed to find FFT coefficients. Then, with an appropriate interval length, ∆x, the new continuous function can be constructed. Determining ∆x is the key to obtaining an accurate approximation of the original function, given that there are sufficient number of data points. 10 1 original function appr. function, m=1 appr. function, m=2 0.8 0.6 0.4 f(x) 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 0 2 4 6 8 10 x 12 14 16 18 20 Figure 6. Interpolation of 1000 data points from function f = cos( 7π 8 k). m represents ∆x. 200 Original Data Estimated function 180 160 Number of sunspots 140 120 100 80 60 40 20 0 0 50 100 150 Time 200 250 300 Figure 7. Filtering based on 5% percentage method. Figure 6 is a graph of the original function f = cos( 7π 8 k) and its approximations with different interval lengths ∆x. When ∆x = 1, the approximation matches the original function. 2.2.3. Filtering. In the third part of the sunspot example, only the 10 most significant component functions were taken, and number 10 was indeed arbitrarily determined. The process of choosing certain components over others is called filtering. Filtering can be extremely useful if data points are suspected to include noise. In this implementation, we take a look at two different methods of filtering. - Method 1: The idea of this method is to choose functions whose coefficients contribute more (or less) than a certain percentage of the total sum. For WAVELETS 200 11 Original data Approximated Number of sunspots 150 100 50 0 −50 0 50 100 150 Time 200 250 300 Figure 8. The original data and its approximation, resulted from band-pass filtering method. example, given a function F , we sum up all frequency contents: S= N/2 X k=1 |F (k)|. Then for each k, we determine how much its frequency content contributes to the sum in percentage terms: c(k) = |F (k)|/S × 100%. Figure 7 displays the original function f and its approximation after filtering all component functions whose contribution is less than 5%. f is the function of sunspot data. - Method 2 (band-pass method): The second method is called band-pass filtering, based on the distribution of frequencies in the signal. The idea is to choose (or to not choose) only the frequencies that lie within one standard deviation from the mean. Using the same set of data, Figure 8 depicts the original function f and its approximation function. Figure 9 displays the original power spectrum, the filtered power spectrum, and the filtered part. Note that the latter two add up to the original power spectrum. 2.3. The Drawbacks of Fourier Analysis. Despite its convenience in dealing with different groups of data, Fourier analysis still poses problems in many applications. This section will examine two common issues of FFT: leakage and Gibb’s phenomenon. 2.3.1. Leakage. To best illustrate the problem of leakage, we will look at the function 7π y = cos( t). 8 12 4500 Original Within 1STD Outside 1 STD 4000 Frequency content 3500 3000 2500 2000 1500 1000 500 0 0 20 40 60 k 80 100 120 140 Figure 9. Band-pass filtering method. Power spectrum of the original data and the filtered power spectrum. 1 Function y Sample points 0.8 0.6 0.4 y 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 0 10 20 30 t 40 50 60 70 Figure 10. Function y = cos( 7π 8 t) and 64 sample points. N = 64 sample points from 0 to 63 are taken (∆x = 1). Figure 10 provides us with a graph of y and sample points. A problem data collectors usually encounter is that there are fewer data points available than they wish. For other data mining methods, a simple method of WAVELETS 13 a) 64 sample points b) 52 sample points 35 16 30 14 12 Frequency content Frequency content 25 20 15 10 8 6 4 5 0 10 2 0 10 20 30 t 40 50 60 70 0 0 10 20 30 t 40 50 60 70 Figure 11. Power Spectra of the FFT of function y with a) 64 sample points and b) 52 sample points (the rest are zeroed). zeroing the rest of the sample points can be used. However, it is not the same for Fourier transform: Figure 11(a) plots power spectrum of the FFT of y with 64 sample points, and Figure 11(b) plots the same thing, however, the last 10 sample points are zeroed. The frequency contents in Figure 11(a) reflect function y better; while the frequency contents in Figure 11(b) appear to include some noise. From these graphs, we can conclude that adding zeros does not help in case of sample points shortage. This effect is known as “leakage”: while the frequency content of the signal did not change, the power spectrum did [4]. 2.3.2. Gibb’s Phenomenon. In many cases, an attempt to fit a function with discontinuities or steep slopes using Fourier transform fails. The reason is that FFT only yields smooth functions. This creates a problem known as Gibb’s phenomenon. For example, take a look at a square wave function defined as follow: f (x) =  1, if −1, if 2kπ ≤ x ≤ (2k + 1)π (2k + 1)π ≤ x ≤ 2(k + 1)π for any nonnegative integer k. Then f can be written as sum of odd harmonics [4]: f (x) = sin(x) + = ∞ X k=1 sin(3x) sin(5x) + +... 3 5 1 sin((2k − 1)x). 2k − 1 The larger k is, the better the function f is approximated. Figure 12 illustrates this fact: even with large k, the values at the points of discontinuities cannot be approximated. Thus, a perfect square wave can never be obtained. 14 a) Small k b) Larger k 1 1 k=1 k=3 k=5 k=7 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 −0.2 −0.2 −0.4 −0.4 −0.6 −0.6 −0.8 −1 −0.8 0 500 1000 1500 2000 2500 3000 3500 4000 4500 −1 5000 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Figure 12. Gibb’s phenomenon. Approximation of function f with a) small k and b) larger k 3. Wavelets 3.1. The Transitions from Fourier Analysis. The traditional Fourier transform is implemented only on the interval [0, 2π]. This suggests that Fourier transform will encounter difficulty in coping with large data files with complicated patterns. Fourier analysis has been modified in many cases to comply with the available data. For example, windowing - the method of multiple Fourier transform - is often adopted to solve the problem of abundant data points [4]. 3.1.1. Lag Matrix. Multiple transform is the method where we split a large data file into equal time intervals and apply Fourier transform on each of them. Dividing a data file into smaller intervals requires two specifications: (i) the time interval, and (ii) the number of data points between two consecutive time intervals that will be skipped (overlaps are encouraged for the purpose of precision). To make it easier for computation, data can be put into a matrix, whose column vectors act as an individual data set for a Fourier transform. Such a matrix is called a lag matrix. Another way to visualize the lag matrix is to slide along the data file a window whose length is equal to the time interval. Notice that the window is not slid smoothly but discretely in accordance with the number of skipped data points. The sliding process is called windowing process. For example, consider a time series consisting of 2400 data points. They can be split into 600 individual data sets, each of which contains 16 data points. The skipping step is therefore 4. Hence the data matrix   x1 x5 x9 . . . x2385  x2 x6 x10 . . . x2386    X= . . ..  ..  . x16 x20 x24 . . . x2400 3.1.2. Sunspot Example Revisited. We will apply the windowing process to the sunspot data set: the data was cut into lengths of 32 time units. Recall that power spectrum is a plot of k against its frequency contents. Figure 13 is the spectra of multiple Fourier transforms applied on the sunspot data set, with skipping WAVELETS 15 1000 5 800 10 600 15 20 400 25 200 30 50 100 150 200 250 0 Figure 13. Power spectra of the windowed sunspot data, data length = 32, skipping step = 1. 1100 1000 5 900 800 10 700 600 15 500 20 400 300 25 200 100 30 5 10 15 20 25 30 35 40 45 50 Figure 14. Power spectra of the windowed sunspot data, data length = 32, skipping step = 5. step equal to 1. Figure 14 is the spectra of the same data set; however, the skipping step is now 5. The color bar next to each graph is the guide to decipher power spectrum in each smaller data set. The symmetric property is retained in both graphs. Figure 13 has better resolution, which results from smaller skipping step. However, both graphs represent the same trend of frequency contents. 3.2. Haar Wavelets. Although multiple transforms can be adopted to solve problems of large data files, such a method is not complete. Determining the time 16 Figure 15. Data from seismic survey [6]. interval and the skipping step is not an easy task for a complicated data set. Moreover, for a data set whose patterns vary greatly, a single time frame might not be appropriate. These shortcomings of Fourier transform gave rise to wavelets. Wavelets were first applied in analyzing data from seismic surveys in geophysics and later on in many other signal processing applications [6]. Figure 15 displays sample data from a seismic survey. Using the multiple transform method is certainly not a good approach in this case: there are short-duration and high-frequency mixed with low-frequency bursts; thus, equal-time intervals will not be able to detect both of them. In this section, the Haar transform and Haar wavelets will be introduced. They represent the simplest type of wavelet analysis, and can serve as the prototype for all other wavelet operations [7]. Two components that play primary roles in any wavelet analysis are the scaling function, known as the father wavelet, and the wavelet function, also known as the mother wavelet. Many theorems or results of this section come directly or after a few simple steps from the definitions with some knowledge of linear algebra; thus, some proofs will be omitted. Keep in mind that all signals are plotted on a two dimensional time-axis against a displacement-axis. This section only works with discrete signals. Thus, every signal f is of the form f = (f1 , f2 , . . . , fN ), where N is a positive even integer (similar to FFT, N is ideally a power of 2). The Haar transform decomposes a discrete signal f into two subsignals: one reserves the trend of f , and the other reserves its fluctuation [7]. 3.2.1. Haar transform, 1-level. Definition 1. The first level of the Haar transform is the mapping H1 defined by H 1 f 7−→ (a1 |d1 ), where f is a discrete signal, a1 = (a1 , a2 , . . . , aN/2 ), and d1 = (d1 , d2 , . . . , dN/2 ) such that (6) am = (7) dm = f2m−1 + f2m √ , 2 f2m−1 − f2m √ . 2 WAVELETS 17 Figure 16. The graph of signal f on the left and its Haar transform, 1-level on the right [7]. For example, let f = (4, 6, 10, 12, 8, 6, 5, 5), then we have f a1 d1 4 6 10 √ √ 12 8 11√ 2 5 √2 − 2 − 2 √ 5 √ 6 5 7√ 2 5 2 2 0 Thus, √ √ √ √ √ √ √ H1 (4, 6, 10, 12, 8, 6, 5, 5) 7−→ (5 2, 11 2, 7 2, 5 2| − 2, − 2, 2, 0). The mapping H1 has an inverse; given (a1 |d1 ), f can easily be obtained: f2m−1 = am + d m am − d m √ √ , f2m = . 2 2 If f is extracted from a continuous signal and the spaced sample values of f are small, each value of d1 will be significantly smaller than the magnitude of f . Figure 16 illustrates this idea. While a1 appears to be a compression of f (a1 is the left half of the graph in Figure 16(b)), d1 (the right half of the same figure) is close to 0 in magnitude. This property is called the small fluctuations feature. Consequently, signals with small fluctuations feature can be compressed using a smaller number of bits. For instance, if small fluctuations in d1 are disregarded, we obtain a compression that is 50% of the original or a 2:1 compression. 3.2.2. Conservation and Compaction of Energy [7]. Definition 2. The energy of a signal f is defined by Ef = N X fi2 . i=1 Theorem 3 (Conservation of Energy). The 1-level Haar transform conserves energy, i.e., E(a1 |d1 ) = Ef for all signal f . 18 Proof. We have N/2 Ea1 |d1 = X N/2 a2i + i=1 = i=1 N/2 = X d2i i=1 N/2 h X X f2m−1 + f2m 2  f2m−1 − f2m 2 i √ √ + 2 2 2 2 f2m−1 + f2m i=1 = Ef
Hence the proof is complete.  General Principle - Compaction of Energy: The energy of the trend subsignal a1 accounts for a large percentage of the energy of the transformed signal (a 1 |d1 ). This is consistent with the small fluctuations feature of d1 mentioned in the previous subsection. 3.2.3. Haar transform, multiple levels. The multiple level Haar transform is performed by consecutively applying Formula (1) on the trend subsignal of the previous operation. For instance, a1 is obtained from the 1-level Haar transform applied on the signal f . Computing the second trend a2 and the second fluctuation d2 for a1 , we obtain a 2-level Haar transform of f . We can keep on doing this as long as the resulting trend signal allows. The Conservation of Energy Theorem still holds here: E(an |dn |dn −1|...|d2|d1 ) = Ef . Definition 4. The cumulative energy profile of a signal f is a signal defined by f2 f2 + f2 f2 + f2 + f2  1 2 2 3 , 1 , 1 ,...,1 . Ef Ef Ef Figure 17 is a graph of f , its 2-level Haar transform, and their cumulative energy profiles. We can see that the elements in the cumulative energy profile of the 2level Haar transform signal approach 1 much more rapidly than the original signal, suggesting that the compression by Haar transform conserves most of the energy of the original signal. 3.2.4. Haar Wavelets. This section introduces two major concepts in wavelet analysis: the wavelet and the scaling signal. They are extremely useful in mathematical operations of wavelet analysis, which will be illustrated in the following sections. 1 Definition 5 (Haar Wavelets). 1-level Haar wavelets W11 , W21 , . . ., WN/2 are 1 1 1 1 defined as Wn = (wn1 , wn2 , . . . , wnN ) such as 1 wni =    √1 2 −1 √ 2 0 if i = 2n − 1 if i = 2n otherwise. Definition 6 (Haar scaling signals). 1-level Haar scaling signals V 11 , V21 , . . ., 1 1 1 1 VN/2 are defined as Vn1 = (vn1 , vn2 , . . . , vnN ) such as WAVELETS Figure 17. The graph of a) signal f , b) 2-level Haar transform of signal f , c) the cumulative energy profile of signal f and d) the cumulative energy profile of its 2-level Haar transform[7]. 1 vni =    √1 2 √1 2 0 if i = 2n − 1 if i = 2n otherwise. So, the 1-level Haar wavelets are: W11 W21 1 WN/2 1 1 = ( √ , − √ , 0, 0, . . . , 0), 2 2 1 1 = (0, 0, √ , − √ , 0, 0, . . . , 0), 2 2 .. . 1 1 = (0, 0, . . . , 0, √ , − √ ), 2 2 and the 1-level Haar scaling signals are: V11 V21 1 VN/2 1 1 = ( √ , √ , 0, 0, . . . , 0), 2 2 1 1 = (0, 0, √ , √ , 0, 0, . . . , 0), 2 2 .. . 1 1 = (0, 0, . . . , 0, √ , √ ). 2 2 19 20 A relationship can be established between subsignals of f and the 1-level Haar wavelets and scaling signals, using the familiar scalar product (known as the dot product in Linear Algebra): am dm 1 = f · Vm , 1 = f · Wm . For any type of wavelets, the mother wavelets act as a window sliding along the signal (translation), whereas the scaling signals allow zooming in and out at each point (dilation). m m m We can also define multiple level Haar wavelets Wnm = (wn1 , wn2 , . . . , wnN ) and m m m m scaling signals Vn = (vn1 , vn2 , . . . , vnN ) in a similar manner:  1 √ if 2m−1 n + 1 ≤ i ≤ 2m n m 2 vni = 0 otherwise,  1 √ if 2m−1 n + 1 ≤ i ≤ 3 × 2m−2 n  2 m −1 √ wni = if 3 × 2m−2 n + 1 ≤ i ≤ 2m n  2 0 otherwise. For example, 1 1 1 1 V22 = (0, 0, 0, 0, √ , √ , √ , √ , 0, 0, . . . , 0), 2 2 2 2 1 1 −1 −1 W22 = (0, 0, 0, 0, √ , √ , √ , √ , 0, 0, . . . , 0). 2 2 2 2 Consequently, we can obtain multiple level Haar transform from multiple level Haar wavelets and scaling signals:
m am = f · V1m , f · V2m , . . . , f · VN/2 m ,
m dm = f · W1m , f · W2m , . . . , f · WN/2 m . The multiple level Haar wavelets and scaling signals are essential in Haar wavelet analysis. Similar to the coefficients in Fourier transform, they provide the means to analyze the signals. More details will be provided in later sections. 3.3. Daubechies Wavelets. Haar wavelets are considered the most basic of all. In 1988, Ingrid Daubechies discovered another family of wavelets that were named after her [8]. Unlike Haar wavelets, Daubechies wavelets are continuous. Consequently, they work better with continuous signals. They also have longer supports, i.e. they use more values from the original signals to produce averages and differences. These improvements enable Daubechies wavelets to handle complicated signals more accurately. We will first examine the simplest of the Daubechies family of wavelets: the Daub4 wavelets. Although the scaling and wavelet numbers are different, the idea of the Daubechies wavelets and Daubechies wavelet transform are very similar to that of the Haar wavelets. Since all wavelet analyses are similar in definitions and properties, the following sections on different types of wavelets will not go into particulars. The previous section on Haar wavelet analysis might be useful as a reference. WAVELETS 21 3.3.1. Definitions. The Daub4 wavelets use four coefficients for their scaling signals and wavelets, compared to two in that of Haar wavelets. The scaling coefficients of Daub4 wavelets are √ √ √ √ 1+ 3 3+ 3 3− 3 1− 3 √ , α2 = √ , α3 = √ , α4 = √ , α1 = 4 2 4 2 4 2 4 2 and the wavelet coefficients are: β1 = α 4 , β2 = −α3 , β3 = α 2 , β4 = −α1 . Then, the first level Daub4 scaling signals are Vn1 = (v1 , v2 , ...vN ), for n = 1, 2, ..., N/2, in which [7]  α1 if i = 2n − 1     if i = 2n  α2 α3 if i = (2n + 1) mod N vi =   α if i = (2n + 2) mod N  4   0 otherwise. The first level Daub4 wavelets Wn1 = (w1 , w2 , ...wN ), for n = 1, 2, ..., N/2, are defined similarly [7]:  if i = 2n − 1   β1   if i = 2n  β2 β3 if i = (2n + 1) mod N wi =   β if i = (2n + 2) mod N  4   0 otherwise. According to the definition, we have the 1-level Daub4 wavelets: W11 W21 W31 1 WN/2−1 1 WN/2 = (β1 , β2 , β3 , β4 , 0, 0, . . . , 0) = (0, 0, β1 , β2 , β3 , β4 , 0, 0, . . . , 0) = (0, 0, 0, 0, β1, β2 , β3 , β4 , 0, 0, . . . , 0) .. . = (0, 0, . . . , 0, β1 , β2 , β3 , β4 ) = (β3 , β4 , 0, 0, . . . , 0, β1 , β2 ), and the 1-level Daub4 scaling signals: V11 = (α1 , α2 , α3 , α4 , 0, 0, . . . , 0) V21 V31 = (0, 0, α1 , α2 , α3 , α4 , 0, 0, . . . , 0) 1 VN/2−1 1 VN/2 = (0, 0, 0, 0, α1, α2 , α3 , α4 , 0, 0, . . . , 0) .. . = (0, 0, . . . , 0, α1 , α2 , α3 , α4 ) = (α3 , α4 , 0, 0, . . . , 0, α1 , α2 ). Figure 18 depicts the two wavelets: Haar and Daub4. Their shapes explain why Haar wavelet works better with discontinuous functions while Daub4 wavelet has an advantage while working with continuous ones. 22 Figure 18. (a) Left figure: Haar wavelet, (b) Right figure: Daub4 wavelet Notice that the scaling signals and wavelets are orthogonal to each other. This can be easily proved using the scaling and wavelet numbers’ properties stated in the next section. 3.3.2. The coefficients’ properties. The wavelet and scaling coefficients satisfy conditions that are essential for the properties of wavelets α21 + α22 + α23 + α24 = 1, √ α1 + α2 + α3 + α4 = 2, β12 + β22 + β32 + β42 = 1, β1 + β2 + β3 + β4 = 0. These identities can be used to prove the orthogonality among all scaling signals and wavelets. 3.3.3. Daub4 First and Multiple-Level Transforms. Much like 1-level Haar transD1 form, Daub4 transform is defined as the mapping f 7−→ (a1 |d1 ). Each element in 1 the first trend subsignal a = (a1 , . . . , aN/2 ) is the scalar product ai = f · Vi1 . Similarly, the fluctuation subsignal d1 = (d1 , . . . , dN/2 ) is the scalar product di = f ·Wi1 . 0 Let us define the elementary signals V10 , V20 , ..., VN as V10 = (1, 0, 0, . . . , 0) V20 = (0, 1, 0, 0 . . . , 0) .. . = (0, 0, . . . , 0, 1). 0 VN We recognize that 1 Vm 0 0 0 0 = α1 V2m−1 + α2 V2m + α3 V2m+1 + α4 V2m+2 , 1 Wm 0 0 0 0 = β1 V2m−1 + β2 V2m + β3 V2m+1 + β4 V2m+2 , where the sub index is mod N . Higher level Daub4 transforms are obtained by applying the 1-level Daub4 transform consecutively on the trend subsignal of the previous level transform. The higher level scaling signals and wavelets are defined accordingly: (8) k Vm k Wm k−1 k−1 k−1 k−1 = α1 V2m−1 + α2 V2m + α3 V2m+1 + α4 V2m+2 k−1 k−1 k−1 k−1 = β1 V2m−1 + β2 V2m + β3 V2m+1 + β4 V2m+2 . WAVELETS 23 Figure 19. (a) Signal A. (b) 2-level Daub4 transform. (c) and (d) Magnifications of the signal’s graph in two small squares; the signal is approximately linear [7]. The above formulas can be modified and applied to any kind of wavelets. To summarize, the first and multiple-level Daub4 transforms are achieved from these wavelets and scaling signals:
m am = f · V1m , f · V2m , . . . , f · VN/2 m ,
(9) m dm = f · W1m , f · W2m , . . . , f · WN/2 m , which are similar to those of Haar transforms. 3.3.4. The Daub4 Transform’s Property [7]. Property 1. If a signal f is approximately linear over the support of a k-level Daub4 k wavelet Wm , then the k-level fluctuation value f · Wkm is approximately zero [7]. The support of a k-level wavelet depends on the number k. For example, the 1-level Daub4 wavelet has 4 time-unit support (4 non-zero coefficients); the 2level Daub4 wavelet has 6 time-unit support, etc... Figure 19 illustrates this idea: magnifications of the original signal appear linear; and the fluctuation subsignals seem to be 0. This property is useful in determining whether Daub4 wavelets are adequate for certain applications. Property 2. Similar to Haar transform, Daub4 transforms also conserve energy. The proof for this can be found in [7]. Figure 20 compares the efficiency of Haar transform and Daub4 transform. The upper graphs show that the detail subsignals from Daub4 transform are significantly less than that from Haar transform. As a result, the cumulative energy of Daub4 transform reaches 1 much faster. 24 Figure 20. (a)2-level Haar transform of signal A. (b) 2-level Daub4 transform on the same signal. (c) Cumulative energy profile for the Haar transform in (a). and (d) Cumulative energy profile for the Daub4 transform in (b)[7]. 3.4. Other Daubechies Wavelets. Other Daubechies wavelets are very similar to Daub4 wavelets. Ingrid Daubechies developed two families of wavelets: the DaubJ wavelets (for J = 4, 6, 8, . . . , 20) and the CoifI wavelets (for I = 6, 12, 18, 24, 30) [7]. 3.4.1. DaubJ Wavelets. We provide a general definition for DaubJ wavelets (for J = 4, 6, 8, . . . , 20) as done in [7]. For each J, the scaling numbers αi (for i = 1, 2, ..., J) are computed. The wavelet numbers are then defined accordingly: βi = (−1)i+1 αJ−i . The 1-level Daub4 scaling signals are: Vn1 = (v1 , v2 , ...vN ), for n = 1, 2, ..., N/2, in which  α1 if i = (2n − 1) mod N     α if i = (2n) mod N  2    .. . vi =  αJ−1 if i = (2n + J − 3) mod N     α if i = (2n + J − 2) mod N  J   0 otherwise. WAVELETS 25 The 1-level Daub4 wavelets are:Wn1 = (w1 , w2 , ...wN ), for n = 1, 2, ..., N/2, in which  if i = (2n − 1) mod N   β1   β if i = (2n) mod N  2    .. . wi =  βJ−1 if i = (2n + J − 3) mod N     β if i = (2n + J − 2) mod N    J 0 otherwise. The scaling numbers and wavelet numbers of DaubJ transform still satisfy some identities: α21 + α22 + . . . + α2J−1 + α2J α1 + α2 + . . . + αJ−1 + αJ i i 0 β1 + 1 β2 + . . . + (J − 2)i βJ−1 + (J − 1)i βJ for i = 0, 1, ..., J − 4. = 1, √ = 2, = 0, The multiple level DaubJ transforms are defined similarly to Daub4. The 1-level transform is applied consecutively on the previous level’s trend signal to achieve a higher transform. 3.4.2. Coiflets. Like the DaubJ family, the CoifI family, also known as the “coiflets”, consists of different wavelets defined in a similar manner [7]. Thus, we will examine the representative Coif6, which should give us good understanding of the family in general. The six Coif6 scaling numbers are: √ 1−√ 7 , 16 √ 2 7 14−2 √ , 16 2 α1 = α4 = α2 = α5 = √ 5+√ 7 , 16 √2 1−√ 7 , 16 2 α3 = α6 = √ 14+2 √ 7, 16 √ 2 −3+ √ 7. 16 2 The wavelet numbers are defined based on the scaling numbers: βi = (−1)i+1 αJ−i . Besides common identities, the Coif6 scaling numbers satisfy additional ones: √ 2 α1 + α 2 + α 3 + α 4 + α 5 + α 6 = (−2)i α1 + (−1)i α2 + 0i α3 + 1i α4 + 2i α5 + 3i α6 = 0 for i = 1, 2, β1 + β 2 + β 3 + β 4 + β 5 + β 6 0β1 + 1β2 + 2β3 + 3β4 + 4β5 + 5β6 = 0 = 0. The Coif6 first level scaling signals Vn1 = (v1 , v2 , ...vN ), and wavelets Wn1 = (w1 , w2 , ...wN ), for n = 1, 2, ..., N/2, are also determined slightly differently [7]:  α1 if i = (2n − 3) mod N     α if i = (2n − 2) mod N  2    if i = (2n − 1) mod N  α3 α4 if i = (2n) mod N vi =   α if i = (2n + 1) mod N  5    α if i = (2n + 2) mod N  6   0 otherwise, 26  β1     β2      β3 β4 wi =    β5     β6   0 if i = (2n − 3) mod N if i = (2n − 2) mod N if i = (2n − 1) mod N if i = (2n) mod N if i = (2n + 1) mod N if i = (2n + 2) mod N otherwise. Formulas (8) and (9) may be employed to find higher level wavelets, scaling signals, and to perform Coiflet transforms at different levels. Coiflets share all properties with Haar and DaubJ wavelets. 3.5. Wavelet Applications. In this section, some basic applications of wavelet analysis will be introduced. The Haar wavelet will be examined in all applications as a model for other types of wavelets. Since wavelet analytical mechanisms are similar across different families of wavelets, more complex wavelets will be examined without elaborate explanation. Juxtaposition of two or more different wavelet analyses in one application will help indicate one wavelet’s advantages over the others. 3.5.1. Multiresolution Analysis. Since discrete signals are subjects of wavelet analysis in this paper, all elementary algebraic operations such as addition, subtraction, and scalar multiplication can be performed on any two or more signals. Multiresolution analysis allows the original signal to be built up from lower resolution signals and necessary details. Definition 7 (First Signals). The First Average Signal A1 is defined by a a a a aN/2 aN/2  1 1 2 2 A1 = √ , √ , √ , √ , . . . , √ , √ . 2 2 2 2 2 2 The First Detail Signal D1 is defined by D1 =  d −d d −d dN/2 −dN/2  1 1 2 2 √ , √ , √ , √ ,..., √ , √ . 2 2 2 2 2 2 Recalling the elementary signals defined in the previous section, we have: (10) f= N X fi Vi0 . i=1 The above formula is called the natural expansion of a signal f in terms of the 0 natural basis of signals V10 , V20 , . . ., VN . WAVELETS 27 Figure 21. The graph of signal A built up from 10-level Haar MRA. Ten averaged signals from A10 to A1 are displayed from top to bottom, from left to right [7]. It follows that f a a a a aN/2 aN/2  1 1 2 2 √ , √ , √ , √ ,..., √ , √ 2 2 2 2 2 2  d −d d dN/2 −dN/2  1 1 2 + √ , √ , √ ,..., √ , √ 2 2 2 2 2 = A1 + D1 N/2 N/2 X X di Wi1 ai Vi1 + = = i=1 i=1 = N/2 N/2 X X

f · Wi1 Wi1 . f · Vi1 Vi1 + i=1 i=1 This is the first level of Haar multiresolution analysis (MRA). Since the multiple level Haar transform can be applied consecutively on average subsignals, we can expand further to obtain multiple level Haar MRA: f = Ak + Dk + . . . + D 2 + D1 , in which N/2k A k = X
f · Vik Vik X
f · Wik Wik . i=1 N/2k Dk = i=1
k
k The values f · Vi and f · Wi are called wavelet coefficients. Each of the component signals has lower resolution than f ; however, if a high-enough level of MRA is employed, the original f can be obtained. In Figure 21, the original signal was achieved after 10 levels of Haar MRA. Since the original signal is continuous, Daubechies wavelets yield better results. Figure 22 and Figure 23 show that more complex wavelets approach the original signal after fewer steps of MRA. 28 Figure 22. Daub4 MRA of the same signal. The graph are of 10 averaged signals A10 through A1 [7]. Figure 23. Daub20 MRA of the same signal. The graph are of 10 averaged signals A10 through A1 [7]. 3.5.2. Compression of Audio Signals. There are two basic categories of compression techniques: lossless compression and lossy compression [7]. A lossless compression technique yields a decompression free in error from the original signal, while a decompression resulted from lossy compression suffers a degree of inaccuracy. However, the lossy compression usually succeeds more often at reducing the size of the data set. The wavelet transform is a lossy compression technique. Method of Wavelet Transform Compression [7] Step 1. Perform a wavelet transform of the signal Step 2. Set equal to 0 all values of the wavelet transform which are insignificant, i.e., which lie below some threshold value. Step 3. Transmit only the significant, non-zero values of the transform obtained from Step 2. This should be a much smaller data set than the original signal. Step 4. At the receiving end, perform the inverse wavelet transform of the data transmitted in Step 3, assigning zero values to the insignificant values which were not transmitted. This decompression step produces an approximation of the original signal. Figures 24 and 25 are two examples of compression using wavelet transform method. Because of the nature of Haar wavelets, discrete signals like signal 1 in Figure 24 can be more easily compressed with high degree of accuracy, while continuous signals such as signal 2 in Figure 25 are much harder to compress, WAVELETS 29 Figure 24. a)Original signal 1, b)10-level Haar transform of signal 1, c) energy map of Haar transform, and d) 20:1 compression of signal 1, 100% of energy[7]. Figure 25. a)Signal 2, b)12-level Haar transform of signal 2, c) energy map of Haar transform, and d) 10:1 compression of the original signal, 99.6% of energy of signal 2[7]. and when they are compressed, the decompresion doesn’t yield high precision even though a smaller compression size and higher level Haar transform were used. Although Haar transforms are prominent in compressing piecewise constant signals, Daubechies transforms work much better with continuous signals. Figure 26 30 Figure 26. (a) The original signal. (b) 12-level Coif30 transform. (c) Energy map of the transform and (d) 32:1 compression of the signal [7]. is the 12-level Coif30 transform of signal 2. We see that a 32:1 compression of this signal by Coif30 transform yields better result when decompressed than the 10:1 compression by Haar transform. 3.5.3. Removing Noise from Audio Signals. When transmitted over a distance, signals are usually contaminated with noise, “the undesired change that has altered the values of the original signal” [7]. Noise is often encountered in three types [7]: (1) Random noise. The noise signal is highly oscillatory, its values alternating rapidly between values above and below an average, or mean, value. We will be mostly working with random noise. (2) Pop noise. This type of noise is heard on old analog recordings obtained from phonograph records. The noise is perceived as randomly occuring, isolated “pops.” (3) Localized random noise. Sometimes the noise appears as in type 1, but only over a short segment or segments of the signal. This can occur when there is a short-lived disturbance in the environment during transmission of the signal. The simple model for a contaminated signal is given as follows: contaminated signal = original signal + noise Denoting f as the contaminated signal, s as the original signal and n as noise, we have an equation f = s + n. To filter out random noise, the threshold method of wavelet denoising can be implemented. In this method, only the transform values whose magnitudes are greater than a threshold Ts > 0 will be kept. Equivalently, we can discard all WAVELETS 31 Figure 27. a) Signal B, 210 values. b) 10-level Haar transform of signal B. The two horizontal lines are at values of ±.25 (the denoising threshold). c) Thresholded transform. d) Denoised signal [7]. the transform values whose magnitudes lie below a noise threshold Tn satisfying Tn < T s . The Root mean Square Error (RMSE) is used to measure the effectiveness of noise removal method: s PN 2 i=1 (fi − si ) RM SE = N s PN 2 i=1 (ni ) = N √ En = √ . N Smaller RMSE indicates a better denoising result. This is similar to the least square method in determining error for a set of data. Figure 27 is a denoising example of signal 1. Part (a) of Figure 27 suggests that some random noise was added to the original signal. Part (b) shows the denoising threshold to be used on 10-level Haar transform of the contaminated signal, part (c) shows the Haar transform after the noise was filtered, and part (d) is the denoised signal. The result is fairly consistent with the original signal given in the previous report. The RMSE between signal B and signal 1 is 0.057. After denoising, the RMSE reduces to 0.011. Figure 28 is another example of denoising method applied on signal 2. For this signal, 212 values were used, together with a higher level of Haar transform (12 instead of 10), and a smaller denoising threshold (0.2 instead of 0.25). The RMSE between signal C and signal 2 is 0.057. After denoising, the RMSE is 0.035. 32 Figure 28. a) Signal C, 212 values. b) 12-level Haar transform of signal C. The two horizontal lines are at values of ±.2 (the denoising threshold). c) Thresholded transform. d) Denoised signal [7]. Even though a higher level of Haar transform as well as a smaller denoising threshold were used, denoising of signal C still did not yield good results, which implies that Haar wavelet denoising method does not work particularly well with continuous signals such as signal 2. Figure 29 consists of signal 2, its 12-level Coif30 transform, the thresholded transform and the denoised signal. In comparison with the denoised signal by Haar transform , this thresholded Coif30 transform yields much better result. This example shows that Daubechies wavelets once again prove to be a better approach in coping with continuous signals. The thresholded method, however, poses a problem: what should the threshold value be? One method of choosing the threshold is to rely on the mean µ and standard deviation σ of the probability density function [7]. A further explanation can be sought in [7]. 3.6. Other Applications. All examples introduced in the previous section are simple and more direct applications of wavelet analysis. The same analysis can be applied to more complicated data sets, such as sound signals or two dimensional images, to produce high quality compression images, or to filter the unexpected noise from a recording. Although wavelets’ history can be traced back to Haar analysis in 1910, most of the development has been in the last 20 years, starting with Stromberg (1981) and Morlet (1984) [9] while working on seismology. Similar to calculus’ development, wavelets quickly evolved beyond mathematics and proved their efficacy in other natural sciences, including both theoretical fields such as physics, chemistry, biology, and appplied fields such as computer science, engineering, and econometrics [9]. The technique is still being developed and put WAVELETS 33 Figure 29. (a) The familiar signal. (b) 12-level Coif30 transform, with threshold = ±0.2. (c) Thresholded transform. (d) Denoised signal.[7] into trial in different disciplines, among them are statisticians and stock market researchers. 4. Wavelets vs. Fourier Analysis Throughout the previous section, Fourier analysis was occasionally compared to wavelets. This section provides a concrete example inwhich wavelets outperform Fourier analysis in yielding desired results. The basic Fourier transform gives a global picture of a data set’s spectrum, whereas wavelet transforms offer a more flexible way to examine a signal, a function or an image. In addition, wavelet transforms also provide information on where or when each frequency component is occuring. These advantages are especially embraced in studying non-stationary or inhomogeneous objects [9]. Windowed Fourier transform, as mentioned in section 3.1, is capable of obtaining a localized frequency; however, the window width remains the same along the data set, which suggests a limitation if the data set contains jumps with unexpected high or low frequencies. This problem can be resolved by wavelet analysis as it allows the degree of localization to be adjustable. For example, consider the signal f in Figure 30, defined by: f = sin(2πv1 x)e−π( x−0.2 10 0.1 )
x−0.5 10 + sin(2πv1 x) + 2 cos(2πv2 x) e−π( 0.1 )
x−0.5 10 + 2 sin(2πv2 x) − cos(2πv3 x) e−π( 0.1 ) . 34 10 8 6 4 2 0 −2 −4 −6 −8 −10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 30. The original signal f . 350 300 frequency 250 200 150 100 50 0 0 50 100 150 200 250 k 300 350 400 450 500 Figure 31. The power spectrum of signal f . The signal comprises three distinct components; the first one has one dominant frequency v1 , the second one has two dominant frequencies v1 and v2 , and the third one also has two: v2 and v3 . Figure 31 is the Fourier transform of f , and Figure 32 displays the power spectrum of each individual component. Fourier transform only detects the dominant frequencies without telling us where they come in. This problem can be solved by using wavelet analysis. Figure 33 plots the wavelet coefficients of signal f against time and scale (logarithm of frequencies). Both the dominant frequencies and their time of appearance are included in the figure. 5. Conclusion Compared to other mathematical topics, wavelets are in a rather fledgling stage. However, their applications are widespread in many fields, both theoretical and WAVELETS 35 1000 frequency 800 600 400 200 0 0 200 400 600 k 800 1000 1200 0 200 400 600 k 800 1000 1200 0 200 400 600 k 800 1000 1200 1000 frequency 800 600 400 200 0 1000 frequency 800 600 400 200 0 Figure 32. The power spectra for each component. 0.8 Magnitude Magnitude 0.8 0.6 0.6 500 400 0.4 0.4 300 200 500 0.2 400 2 300 8 7 6 4 200 5 4 Scale 3 2 100 1 Time 0.2 100 6 8 Time Figure 33. Wavelet transform of signal f . The plot is captured at two different angles 36 practical. Arising from the Fourier analysis’ failure to cope with large and complex data files, wavelets rapidly develop to resolve these problems. This paper briefly introduced Fourier series and chose Fast Fourier Transform (FFT) to be the representative method for Fourier analysis. A modification of FFT - the windowing process or short-term Fourier transform - was also examined as a transition before moving to wavelets. Since all wavelets are constructed similarly, the first and simplest, the Haar wavelet, was studied first and in more detail. Two other families of wavelets were mentioned in the paper: DaubJ wavelets and Coiflets. Some simple applications of wavelet analysis mentioned in the paper include multiresolution analysis, the denoising problem and compression of audio signals. These applications also provide a basis to make a comparison between Fourier and wavelet analyses. All examples in this paper, including discrete data and audio signals, are onedimensional. However, wavelet analysis is capable of dealing with higher dimensional data sets, such as pictures. Information on wavelet analysis on 2D data sets can be found in [7] or most books on applications of wavelets. More complex applications of wavelets, such as signal detection and applications in statistics, can also be topics for further study. References [1] [2] [3] [4] [5] [6] [7] [8] [9] WikiPedia, Fourier Series. http://en.wikipedia.org/wiki/Fourier_series WikiPedia, Fourier Series. http://en.wikipedia.org/wiki/Wavelets Stewart, James. Calculus Early Transcendentals, 4th Edition. Fourier Analysis, Course Notes, Modeling II. Sunspots. Retrieved from website http://www.exploratorium.edu/sunspots/activity. html. Boggess, A. and Narcowich, F., A First Course in Wavelets with Fourier Analysis. New Jersey 2001. Walker, J. , A Primer on Wavelets and their Scientific Applications. CRC Press LLC 1999. Riddle, L. Ingrid Daubechies. Biographies of Women Mathematicians. Retrieved on Dec 1, 2005 from website http://www.agnesscott.edu/lriddle/women/daub.htm Abramovich, F., Bailey, T. and Sapatinas, T.. Wavelet Analysis and Its Statistical Applications. The Statistician, Vol. 49, No. 1 (2000), pp 1-29.