-
Notifications
You must be signed in to change notification settings - Fork 156
/
wavelet3.html
180 lines (150 loc) · 6.25 KB
/
wavelet3.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
<HTML>
<HEAD>
<TITLE>Wave:Wavelets 3</TITLE>
<link rel="icon" href="http://atoc.colorado.edu/research/wavelets/favicon.ico" type="image/x-icon" />
<link rel="shortcut icon" href="http://atoc.colorado.edu/research/wavelets/favicon.ico" type="image/x-icon" />
</HEAD>
<BODY BGCOLOR="white">
<H1>Wavelet Analysis</H1>
<H2>
<UL>
<LI><A HREF="wavelet1.html">Introduction</A>
<LI><A HREF="wavelet2.html">Wavelets</A>
<LI><FONT COLOR="green">Algorithms</FONT>
<LI><A HREF="montecarlo.html">Monte Carlo</A>
<LI><A HREF="references.html">References</A>
</UL>
</H2>
<HR><!----------------------------------------------------------------->
<P>
<H2>
Scales
</H2>
      In the <A HREF="wavelet2.html">previous section</A>
we saw how one can use a typical wavelet (the Morlet) to decompose
a time series into time-frequency phase space. We now turn to the actual
computation of the wavelet transform.
<P>
      The first thing to consider is the shape of the wavelet.
For decomposing the NINO3 SST data, we chose the Morlet wavelet because:
<OL>
<LI>it is commonly used,
<LI>it's simple,
<LI>it looks like a wave.
</OL>
There are an infinite number of other mother wavelets that could be chosen
(see Farge 1992 for examples).
<P>
      For the Morlet wavelet transform, where the mother
wavelet is:
<P ALIGN=center>
<IMG SRC="images/wl_eqn2.1.gif" ALT="Eqn 2.1" ALIGN=center WIDTH=228 HEIGHT=28>
</P>
we must first choose the wavenumber <I>w<SUB>0</SUB></I>, which gives
the number of oscillations within the wavelet itself. One condition
of the wavelet transform is that the average of the wavelet itself
must be zero. In practice, if we choose <I>w<SUB>0</SUB></I>=6, then
the errors due to non-zero mean are smaller than the typical computer
round-off errors (Farge 1992).
<P>
      Our next choice is a set of scaling parameters <I>s</I>,
such that we adequately sample all the frequencies present in our
time series. We first choose the smallest resolvable scale,
<I>s</I><SUB>0</SUB>, as some multiple of our time resolution, <I>dt</I>.
For the NINO3 SST data, we have seasonal data, thus <I>dt</I> = 0.25 years.
The smallest wavelet we could possibly resolve is 2<I>dt</I>, thus
we choose <I>s</I><SUB>0</SUB> = 2<I>dt</I> = 0.5 years.
<A NAME="eqn3.1"><A NAME="eqn3.1b">
The larger scales (longer periods) are chosen as power-of-two multiples
of this smallest scale,
<P>
<FONT SIZE=4 COLOR=00FF00>(3.1a)</FONT> . . . . . . . . . . </A>
<IMG SRC="images/wl_eqn3.1.gif" ALIGN=center WIDTH=283 HEIGHT=25>
<P>
<FONT SIZE=4 COLOR=00FF00>(3.1b)</FONT> . . . . . . . . . . </A>
<IMG SRC="images/wl_eqn3.1b.gif" ALIGN=center WIDTH=229 HEIGHT=24>
<P>
The largest scale chosen should be less than 1/2 the length of the
entire time series. The choice of scales for the NINO3 SST data is
shown on the right-hand axis of <A HREF="wavelet2.html#fig3">Figure 3</A>.
In theory, if these scales are chosen wisely, then one
can construct an orthogonal complete basis set. In reality, one usually
over-samples the scales so as to provide information on
freqencies in between the orthogonal scales. Thus, in Figure 3, the
scales have been over-sampled by choosing 10 sub-scales within each
scale.
<P>
<H2>
Algorithms
</H2>
<A NAME="eqn3.2">
<A NAME="eqn3.3">
      It is possible to compute the wavelet transform in
the time domain using <A HREF="wavelet2.html#eqn2.3">Eqn 2.3</A>.
However, it is much simpler to use the fact that the wavelet transform
is the convolution between the two functions <I>x</I> and <I>Psi</I>,
and to carry out the wavelet transform in Fourier space using the
Fast Fourier Transform (FFT).
In the Fourier domain, the wavelet transform is simply,
<P>
<FONT SIZE=4 COLOR=00FF00>(3.2)</FONT> . . . . . . . . . . </A>
<IMG SRC="images/wl_eqn3.2.gif" ALIGN=center WIDTH=276 HEIGHT=61>
<P>
<A NAME="eqn3.4">
where the ^ indicates the Fourier transform (FT), and the Fourier transform
of the time series is given by:
<P>
<FONT SIZE=4 COLOR=00FF00>(3.3)</FONT> . . . . . . . . . . </A>
<IMG SRC="images/wl_eqn3.3.gif" ALIGN=center WIDTH=227 HEIGHT=60>
<P>
To use this formula,
the FT of the wavelet function should be known analytically.
In addition, the wavelets must be normalized as:
<P>
<FONT SIZE=4 COLOR=00FF00>(3.4)</FONT> . . . . . . . . . . </A>
<IMG SRC="images/wl_eqn3.4.gif" ALIGN=center WIDTH=258 HEIGHT=55>
<P>
Unlike the convolution, the FFT method allows the computation of
all <I>n</I> points simultaneously, and can be efficiently coded using
any standard FFT package.
<P>
      The steps to compute the wavelet transform for a time
series are thus:
<OL>
<LI>Choose a mother wavelet,
<LI>Find the Fourier transform of the mother wavelet,
<LI>Find the Fourier transform of the time series,
<LI>Choose a minimum scale <I>s</I><SUB>0</SUB>, and all other scales,
<LI>For each scale, do:
<UL>
<LI>Using Eqn. 3.4 (or whatever is appropriate for your
mother wavelet), compute the daughter wavelet at that scale;
<LI>Normalize the daughter wavelet by dividing by the square-root
of the total wavelet variance (the total of
(<I>Psi</I>)<SUP>2</SUP> should then be one, thus preserving
the variance of the time series);
<LI>Multiply by the FT of your time series;
<LI>Using Eqn. 3.2, inverse transform back to real space;
</UL>
<LI>Make a contour plot.
</OL>
<P>
      One problem with performing the wavelet transform in
Fourier space is that this assumes the time series is periodic. The
result is that signals in the wavelet transform at one end of the
time series will get wrapped around to the other end. This effect
is more pronouced at larger scales as the influence of each wavelet
extends further in time. One way to avoid this is to pad one end of the
time series with zeroes. A clever method is to pad with enough zeroes
to make the length of the time series equal to a power of two, and
thereby speed up the FFT as well.
<P>
<CENTER>
<H3>
-- <A HREF="montecarlo.html">on to Monte Carlo</A> --
</H3>
</CENTER>
<HR><!----------------------------------------------------------------->
<A HREF="./">back to Wavelet Home Page</A>
</BODY>
</HTML>