Wavelet Transforms in Time Series Analysis

1 Wavelet Transforms in Time SeriesAnalysisAndrew TangbornGlobal Modeling and Assimilation Office, Goddard Space Flight fourier What is a Wavelet ?3. Continuous and Discrete Wavelet Transforms4. Construction of Wavelets through dilation Example - Haar wavelets6. Daubechies Compactly Supported Data compression, efficient Soft Continuous transform - Morlet Wavelet10. Applications to approximating error correlationsFourier Transforms A good way to understand how wavelets work and why they are useful is bycomparing them with fourier Transforms . The fourier transform converts a time Series into the frequency domain:Continuous Transformof a function f(x): f( ) = Z f(x)e i xdxwhere f( ) represents thestrengthof the function at frequency , where Transformof a function f(x): f(k) = Z f(x)e ikxdxwherekis a discrete discrete dataf(xj),j= 1,..,N fk=NXj=1fje( i2 (k 1)(j 1)/N) The Fast fourier transform (FFT) iso(NlogN) : Single Frequency Signalf(t) = sin(2 t)012345678910 1 Coefficient Imaginary CoefficientFourier transform Discrete fourier transform (DFT) locates the single frequency and re-flection.

2 Complex expansion is not exactly sine Series - causes some Frequency Signalf(t) = sin(2 t) + 2 sin(4 t)012345678910 3 2 10123 Signal012345678910 1 Coefficients Imaginary CoefficientsFourier transform Both signal frequencies are represented by the fourier Signalf(t) = sin(2 t) + 2 sin(4 t) < t < .55f(t) = sin(2 t) otherwise012345678910 3 2 10123 Signal012345678910 coefficients Imaginary CoefficientsFourier transform First frequency is found, but higher intermittent frequency appears asmany frequencies, and is not clearly is a Wavelet ? A function that is localized in time and frequency, generally with a zeromean. It is also a tool for decomposing a signal by the fourier transform :A signal is only decomposed into its frequency information is extracted about location and time. What happens when applying a fourier transform to a signal that hasa time varying frequency?

3 02468101214161820 1 fourier transform will only give some information onwhich frequencies are present, but will give noinformation on when they Representation of DecompositionOriginal SignalTime >Frequency > The signal is represented by an amplitude that is changing intime. There is no explicit information on TransformTime > Frequency > The fourier transform results in a representation that depends only on frequency. It gives no information on TransformTime >Frequency > The Wavelet transform contains information on both the timelocation and fre-quency of a typical (but not required) properties of wavelets Orthogonality - Both Wavelet transform matrix and Wavelet functions can for creating basis functions for computation. Zero Mean (admissibility condition) - Forces Wavelet functions towiggle(oscillatebetween positive and negative). Compact Support - Efficient at representing localized data and are Wavelets Defined?

4 Families of basis functions that are based on(1) dilations: (x) (2x)(2) translations: (x) (x+ 1)of a given general mother Wavelet (x). How do we use (x)?The general form is: jk(x) = 2j/2 (2jx k)wherej: dilation indexk: translation index2k/2needed for normalization How do we get (x)?Dilation EquationsConstruction of Wavelets We consider here only orthogonally/compactly supportedwavelets- Orthogonality means: Z jk(x) j k (x)dx= kk jj Wavelets are constructed from scaling functions, (x): (x) come from the dilation equation: (x) =Xkck (2x k)ck:Finiteset of filter coefficients General features:- Fewer non-zerock s meanmorecompact andlesssmooth functions- More non-zerock s meanlesscompact andmoresmooth functionsRestrictions on the Filter Coefficients Normality: Z (x)dx= 1 Z Xkck (2x k)dx= 1 Xkck Z (2x k)dx= 1 Xkck Z 12 (2x k)d(2x k) = 1 Xkck Z 12 ( )d( ) = 1 Xkck(12)(1) = 1orXkck= 2 Simple Examples- Smallest number ofckis 1: just get (x) = , zero Haar Scaling Function:c0= 1,c1= Scaling FunctionThe scaling function equation is: (x) = (2x) + (2x 1)The only function that satisfies this is.

5 (x) = 1 if 0 x 1 (x) = 0 otherwise 1 Translation and Dilation of (x): (2x) = 1 if 0 x 1/2 (2x) = 0 elsewhere (2x 1) = 1 if 1/2 x 1 (2x 1) = 0 elsewhereso the sum of the two functions is then: 1 Wavelet Wavelets are constructed by taking differences of scalingfunctions (x) =Xk( 1)kc1 k (2x k)differencing is caused by the ( 1)k:so the basic Haar Wavelet is: 1 1 the family comes from dilating and translating: jk(x) = 2j/2 (2jx k)so that thej= 1 wavelets are: 1 1 of Haar wavelets Translation no overlap 1 1 Dilation cancellation 1 1 Compactly Supported WaveletsThe following plots are wavelets created using larger numbers of filter coefficients, but having allthe properties of orthogonal wavelets. For example,D4has coefficients:ck=14(1 + 3),14(3 + 3),14(3 3),14(1 3)(Reference: Daubechies, 1988)020406080100120 and Discrete Wavelet TransformsContinuous(Twavf)(a,b) =|a| 1/2 Zdtf(t) (t 1b)wherea: translation parameterb: dilation parameterDiscreteTwavm,n(f) =a m/2oZdtf(t) (a mot nbo)m: dilation parametern: translation parameterao,bodepend on the Wavelet usedFast Wavelet transform (Reference: S.)

6 Mallat, 1989) Uses the discrete data:hf0f1f2f3f4f5f6f7i Pyramid Algorithm o(N) !! - Start at finest scale and calculate differencesj=0j=1j=2f1f2f3f4f5f6f7f8a0,0 b0,0a1,0a1,1b1,0b1,1a2,0a2,1a2,2a2,3b2,0 b2,1b2,2b2,3and averages- Use Averages at next coarser scale to get new set of differences (bj,k) and averages(aj,k)And the coefficients are simply the differences (bj,k) and the average for the coarsestscale (a0,0):ha0,0b0,0b1,0b1,1b2,0b2,1b2,2b2,3 iFor Haar Wavelet :aj 1,k=c0aj,k+c1aj,k+1bj 1,k=c0aj,k c1aj,k+1- The differences are the coefficients at each scale. Averages used for next scale. How do we store all this? (eg, Numerical Recipes algorithm).Intermittent Signal - Wavelet Transform012345678910 3 2 10123 Signal020406080100120 5 4 3 2 101234j=2j=3j=4j=5 Higher Frequency intermittent signal shows up in thej= 5,6scales, and only in themiddle portion of the domain. Localized frequency data is Compression - Efficient Representation012345678910 2 1 Localized Data, 128 2 1 term Wavelet Reconstruction012345678910 1 term fourier Reconstruction Localized function is represented with more accuracy using16 Wavelet coefficientsbecause only significant coefficients need to be retained.

7 fourier coefficients are by Soft ThresholdingBasic Idea: Wavelet Representation highly compresses coherent data into justa few of each coefficient is relatively large. White noise contains energy at all time scales and time is spread to many (if not all) Wavelet coefficients are relatively small in amplitude. How do we remove the noise contribution in the coefficients?Ans: Shrink all of them just a little. How much do we shrink each coefficient?t=q2log(n) / nwheren= number of data points and =noise standard 1 Signal with white noise of known 1 coefficients Shrunk toward zero. Soft Thresholding takes advantage of the fact that white noise is rep-resented equally by all coefficients. Data compression means that coherent signal is representedby just afew coefficients with relatively large values. White noise (or non-whiteas well) is spread over many coefficients, adding just a small amount tothe magnitude of Continuous Morlet Wavelet TransformThe Morlet mother Wavelet is a complex exponential ( fourier ) with aGaussian envelop which ensures localization: (t) =exp(i 0t)exp( t2/2 2)where 0is the frequency and is a measure of the spread or that while the footprint is infinite, the exponential decay createsan effective footprint which is relatively and dilations of the Morlet Wavelet : (b,a)(t) =1aexp"i 0 t ba!

8 #exp t ba!2/2 2 The Morlet Wavelet using matlabIn matlab, the Morlet mother Wavelet can be constructed using thecommand:[psi,x] = Morlet(-8,8,128);on 128 grid points, and domain of [-8,8]. 8 6 4 202468 1 The Morlet transform with MatlabGiven a time Series :y(t) =(sin(2 t) +sin(32 t) t .6sin(2 t)otherwiseon 128 grid points, the Morlet transform can be calculated and plottedusing the command:coef=cwt(y,[1 2 4 8 16 32 64 128], morl , plot ) the top level of the coefficient plot shows bright (or large) values forscale 128 (largest scale). The next two layers are not zero, indicatingleakagebetween differenttime 1 (x)Absolute Values of Ca,b Coefficients for a = 1 2 4 8 16 ..time (or space) bscales a20406080100120 1 2 4 8 16 32 64128 Example: Southern Oscillation Index Time Series 1951-2005 The SOI is the monthly pressure fluctuations in air pressure between Tahiti andDarwin.)

9 It is generally a noisy time Series and benefits from some smoothing. Softthresholding is a way to remove the noise in the signal without removing important in-formation. The multivariate ENSO signal comes from sea level pressure, surface wind,sea surface tempemperature, surface air temperature and total amount of 8 6 4 20246 Time in monthsSOI signalOriginal SOI signal 1951 20050100200300400500600700 8 6 4 20246 Time in MonthsSOI signal Denoised SOI SignalMultivariate ENSO Signal from CDCM orlet Wavelet Decomposition of SOI time series0100200300400500600700 8 6 4 20246 MonthSOIM onthly Southern Oscillation Index, 1951 2006 Absolute Values of Ca,b Coefficients for a = 1 2 3 4 5 ..Time in monthsTime Scale in Months100200300400500600 110016623229836443052962850100150200 Morlet Coefficients for SOItime in monthsTime scale in months100200300400500600 1 19 37 55 73 9110912714516350100150200 Application to Approximation of Error Correlations Error correlations are essential for carrying out atmospheric data assimilation.

10 They tell us how errors in model outputs are spatially related, and how to spread informationfrom observational data into the model. Error correlations are generally the most computationallyand memory intensive parts of a dataassimilation system. Below is the error correlation of a chemical constituent assimilation system (a), and waveletapproximations created by retaining 10 % (b), 10% (c) and 2% (d) of the Wavelet time scales from climate signals:Reference: Seasonal-to-Interannual variability of Ethiopia/Hornof Africa Monsoon. Part 1: Associations of Wavelet Filteredlarge-Scale Atmospheric Circulation and Global Sea SurfaceTemperature, Segelet et al. J. of Climate, can we use information separated out by time-scaleto make predictions?Reference: Webster and Hoyos, BAMS, Vol 85, 2004. Prediction of Asian monsoons on 15-35 day timescale. Morphology of the Monsoon Intraseasonal Oscillation (MISO) Madden-Julian Oscillation (MJO):- Eastward propagating convection- Largest variance in 20-40 day spectral bandFacts about South Asian Summers(1) Convection stronger in eastern Indian Ocean.