1. technology and computing

Journal of Information & Computational Science 11:17 (2014) 6369–6385
Available at http://www.joics.com
November 20, 2014
High Precision Computation of Morlet Wavelet Transform
for Multi-period Analysis of Climate Data ⋆
Hua Yi ∗, Zhiquan Chen ∗, Yanhua Cao ∗
School of Mathematics and Physics, Jinggangshan University, Ji’an 343009, China
Abstract
The algorithms of continuous wavelet transform based on Morlet wavelets are analyzed. The algorithm
analysis, advantage and disadvantage, and remedial method of FFT-based method are given. And then
the algorithm proposed in this paper is free of shortcomings of FFT-based method and its precision
is the best among existing linear-convolution-based methods. Thirdly, two different normalizations
for the computation of continuous wavelet transform are proposed in order to obtain high-precision
characterization of data. Finally, we stress the importance of suitable algorithms by analyzing the
multi-period structure of climate data.
Keywords: Morlet Wavelet; Continuous Wavelet Transform; Entropy; FFT; Convolution; Matlab
1
Introduction
The Continuous Wavelet Transform (CWT) is widely used as a common tool for analyzing localized variations of power within a time series, and therefore the program realization of CWT
for signal analysis is becoming increasing interest. This article focuses on the computation of
CWT with Morlet wavelets based on Matlab (matrix laboratory) for two reasons. Firstly, Matlab is a high-level language and interactive environment and widely used in academic and research institutions. In addition, the algorithms of CWT were developed by many researchers
and the source codes about CWT, written in Matlab language, are freely available in internet
besides cwt function in wavelet toolbox of Matlab. For instance, a Matlab freeware for data
analysis, called JLAB, witten by Lilly, is available at http://www.jmlilly.net. Torrence and Compo (1998) have given a practical step-by-step guide to wavelet analysis, with examples taken
from time series of the El Ni˜
no-Southern Oscillation (ENSO), and developed a code at the URL
http://paos.colorado.edu/research/wavelets/ for computing the continuous wavelet transform of
⋆
Supported jointly by the National Natural Science Foundation of China (No. 11347210), the Natural Science
Foundation of Jiangxi Province (No. 20132BAB201014), and the Doctoral Starting up Foundation of Jinggangshan
University (No. JZB1304).
∗
Corresponding author.
Email addresses: [email protected] (Hua Yi), [email protected] (Zhiquan Chen),
[email protected] (Yanhua Cao).
1548–7741 / Copyright © 2014 Binary Information Press
DOI: 10.12733/jics20104715
6370
H. Yi et al. / Journal of Information & Computational Science 11:17 (2014) 6369–6385
time series [11]. Another freeware, which contains the source program of CWT, known as WaveLab, can be retrieved at: http://www-stat.stanford.edu/∼wavelab, and so on.
Secondly, the complex Morlet wavelets, with the property of optimal localization in time and
frequency domains, have been widely used for modeling multiscale nonstationary processes in
space and time. In 1993, Meyers studied the dispersion of Yanai waves with Morlet wavelets [8].
The wavelet transform with Morlet wavelets has also been used to identify periodic oscillations
of real-life signals [2]-[4], [13]-[15].
There are several different ways to implement a continuous wavelet transform. The simplest
method is numerical integration for the transform using classical quadrature (such as a trapezoidal
rule). The drawback of this technique is that it is time consuming. Unser proposed a fast
noniterative algorithm in which the input signal and the analyzing wavelet are both represented
by polynomial splines for the evaluation of continuous wavelet transforms. However, this method
can only deal with integer scales [12].
To fully exploit the advantages of continuous wavelet transform, FFT-based algorithms, which
allow us to use arbitrary values for the scale variable, are widespreadly employed in aforementioned
Matlab freeware toolboxes, such as JLAB, WaveLab, and Torrence’s codes etc. However, the
method will produce an artificial periodicity in the wavelet coefficients if signal is not periodic
since this method assumes that data is cyclic. Therefore, a data processing method is proposed
to overcome this problem after the principle of FFT-based method is analyzed and the issues will
be addressed in Section 3.2.
At the same time, the function cwt in wavelet toolbox of Matlab, will not suffer this pseudo
periodic problem because the algorithm involves linear convolution instead of circular convolution.
Yet the precision for analyzing data is less than FFT-based method. And the errors caused by
differential coefficient was analyzed [17]; the algorithm was modified by Zhao [17]. The algorithm
is further studied in this paper so that the precision for analyzing data is enhanced, which is
another goal of this paper and will be addressed in Section 3.1.
The third goal, the focus of Section 2, is to propose two normalizations of mother wavelet
and time-frequency atoms ψa,b (t) when conducting continuous wavelet transform so that the
characteristics of signal are accurately depicted with wavelet coefficients.
In Section 4, a physical signal—temperature data is analyzed in order to illustrate the importance of suitable continuous wavelet transform method, which is the fourth goal of this paper.
2
2.1
Morlet Wavelets and Continuous Wavelet Transform
Morlet Continuous Wavelet Transform (MWT)
The Morlet wavelet is a plane wave modulated by a Gaussian function, and is defined as
−t2
ψ(t) = C · e 2σ2 eiηt ,
(1)
where η is the center frequency, σ 2 is the shape parameter to control time and frequency resolution, and C is a constant which should take different values according to different definition of
continuous wavelet transform. The Morlet Wavelet Transform (MWT) of signal f (t) at time b
H. Yi et al. / Journal of Information & Computational Science 11:17 (2014) 6369–6385
6371
and scale a is calculated [6]:
∫
Wf (a, b) =
+∞
−∞
(t − b)
1
f (t) √ ψ ∗
dt,
a
a
where ∗ denotes complex conjugation. Then the constant C should be
(2)
1
1
(σ 2 π) 4
and the mother
wavelet takes the form [6]
ψ(t) =
1
−t2
1
(σ 2 π) 4
· e 2σ2 eiηt
(3)
so that ∥ψa,b (t)∥2 = ∥ψ(t)∥2 = 1, where ψa,b (t) are the family of time-frequency atoms obtained
by scaling ψ by a and translating it by b, that is
1 (t − b)
ψa,b (t) = √ ψ
.
a
a
(4)
.
2.2
Modified Morlet Continuous Wavelet Transform (MMWT)
Different normalization may be found in the literature [1, 16]. If we define the family ψa,b as
ψa,b
1 (t − b)
= ψ
,
a
a
(5)
we can get ∥ψa,b ∥1 = ∥ψ(t)∥1 [1, 15]. Furthermore, if we change the constant C so that the mother
wavelet is defined as
√
−t2
2
iηt
2σ 2 · e
ψ(t) =
·
e
,
(6)
πσ 2
we have ∥ψa,b ∥1 = ∥ψ(t)∥1 = 2 [15]. In that case, we can define the Modified Morlet Wavelet
Transform (MMWT) as
∫ +∞
1 (t − b)
dt.
(7)
Wf (a, b) = ⟨f, ψa,b ⟩ =
f (t) ψ ∗
a
a
−∞
It should be mentioned that the complex Morlet wavelet function known as cmor, is defined in
2
wavelet toolbox as cmor(t) = (π · F b)−0.5 e2iπ·F c·t e−t /F b , where F b is a bandwidth parameter, and
F c is a wavelet center frequency. In order to ensure that ∥cmor(t)∥1 = 2, we modify the cmor in
this paper so that
2
cmor(t) = 2 · (π · F b)−0.5 e2iπ·F c·t e−t /F b .
(8)
For more information on the comparison between MWT and MMWT, we refer the reader to [15].
The algorithms for MWT and MMWT are almost the same except for the constant C of
mother wavelet and the normalizations of time-frequency atoms ψa,b (t). On the other hand,
MMWT maintains the same amplitude between the original frequency component and the wavelet
coefficients, which qualify itself as a superior tool for analyzing data [15, 16]. Therefore, algorithms
for MMWT will be exclusively addressed in the following part of this article.
6372
3
H. Yi et al. / Journal of Information & Computational Science 11:17 (2014) 6369–6385
Algorithms for MMWT
3.1
Time-domain-based Algorithms
The algorithm of cwt in wavelet toolbox, Zhao’s codes, and the algorithm proposed in this article
are all based on the convolution of signal and wavelets in time domain. The proposed algorithm
is first introduced, and then the other two are given.
3.1.1
The Proposed Algorithm
Assume the support of the mother wavelet ψ(t) is [−T, T ]. If define
1 (t)
ψa (t) = ψ
a a
(9)
1 ∗ ( −t )
∗
e
,
ψa (t + a · T ) = ψ a (t) = ψa (−t) = ψ
a
a
(10)
and
the support of ψea (t) = a1 ψ ∗ ( a·Ta−t ) is [0, 2a·T ]. Then the samples of ψea (t) is the impulse response of
a causal system. The wavelet transform, which is basically a linear filter whose response function
is given by the wavelet function, can be rewritten as a convolution:
∫
Wf (a, b) =
+∞
−∞
1 (t − b)
f (t) ψ ∗
dt = f ∗ ψ a (b) = f ∗ ψea (b + a · T ).
a
a
(11)
Denote the sampled signal by {fd (i), i = 0, 1, · · · , N − 1}, with sample frequency 1. As for the
samples of wavelet ψ a (t), we make use of the analytic expression of ψ(t) defined in (6) or cmor
defined in (8) and sample ψ a (t) in [−aT, aT ] with sample frequency 1 . And so we obtain
{
}
ψ a (ti ) : ti = −aT + i, i = 0, 1, · · · , 2a · T ,
(12)
the discrete version of ψ a (t)1 . Then we have the following result corresponding to the discrete
version of (11).
Proposition 1: Let
ha (i) = ψ a (ti − aT ) = ψea (ti ), i = 0, · · · , 2a · T.
(13)
Wf (a, i) = fd ∗ ha (i + aT ), i = 0, 1, · · · , N − 1,
(14)
Then
where ∗ denotes the linear convolution operator.
1
Without loss of generality, assume the sampling rates of data to be analyzed and the wavelet are equal to 1
and a · T is an integer for the convenience of analysis.
H. Yi et al. / Journal of Information & Computational Science 11:17 (2014) 6369–6385
6373
Proof: For i = 0, 1, · · · , N − 1, discretizing (11) yields
Wf (a, i) = fd ∗ ψ a (ti ) =
N
−1
∑
fd (k)ψ a (ti−k )
k=0
=
=
N
−1
∑
k=0
N
−1
∑
fd (k)ψea (aT + ti−k ) =
N
−1
∑
fd (k)ψea (−aT + aT + i − k)
k=0
fd (k)ψea (taT +i−k ) = fd ∗ ψea (taT +i ) = fd ∗ ha (i + aT )
(15)
k=0
In brief, {fd ∗ ha (i), i = 0, 1, · · · , N + 2aT − 1} can be implemented with linear convolution,
for example, the conv function in Matlab. Then {Wf (a, i), i = 0, 1, · · · , N − 1} are just the
{fd ∗ ha (i), i = aT, aT + 1, · · · , aT + N − 1}.
The computational complexity for the convolution-based algorithm depends on how the linear
convolution operation is implemented. In real applications, one would like to use the convolution
theorem together with the FFT in order to perform linear convolutions in the frequency-domain
−1
2aT
with reduced complexity [9]. For the data {fd (i)}N
i=0 and the filter {ha (i)}i=0 . The classical
frequency-domain approach requires to pad both signals with enough zeros to at least a size of
N + 2aT . The complexity can be expressed as a function of N + 2aT . Then, we have that three
FFTs are needed to transform the signals to the frequency-domain and back to the time-domain.
log2 (N + 2aT ) complex
For each of these, the FFT (or IFFT) requires approximately N +2aT
2
multiplications. The point-wise multiplication of the spectra requires N + 2aT multiplications.
)
Hence, the total number of multiplications performed is N +2aT + 3(N +2aT
log2 (N +2aT ) for each
2
scale a. For more information on some other method to perform the linear convolution adapting
to the feature of the data, we refer the reader to [7, 10].
The main source program of this method is given in Appendix A2.
3.1.2
The Algorithm of cwt in the Toolbox of Matlab
According to (7),
Wf (a, b) =
∑∫
k
k+1
k
1 (t − b)
dt.
f (t) ψ ∗
a
a
If f (t) = f (k), for t ∈ [k, k + 1], then
∫ k
∫
(
) )
( ∫ k+1 ( t − b )
1∑
∗ t−b
∗
Wf (a, b) =
dt −
ψ
dt = f ∗ ∆ ψ a (b).
f (k)
ψ
a k
a
a
−∞
−∞
(16)
(17)
According to (17), the convolution in algorithm of cwt in the toolbox of Matlab involves a filter
by first integrating, and then calculating the difference for the filter defined in (10). There is no
error introduced until now because theoretically the composition of integration and differential
will be an identity transform. However, practically the redundant operation will bring in some
error. Zhao et al. proposed to discard directly these redundant operation [17]. The precision of
Zhao’s codes is better than cwt’s.
6374
3.1.3
H. Yi et al. / Journal of Information & Computational Science 11:17 (2014) 6369–6385
The Algorithm of Zhao
In [17], three steps were used to obtain the filter defined in (13). Firstly, by using wavef un
function of Matlab toolbox, the samples of the mother wavelet ψ(t) in [−T, T ] is computed as
{
}N
2T
ψ(tj ), tj = −T +
(j − 1)
.
(18)
N −1
j=1
N , controlled by a precision parameter, may be 210 or 215 or others. In other words, N also
depends on the sampling interval
dw =
2T
N −1
(19)
of the mother wavelet. In this case, except for round-off error, there is no other error introduced
for the samples because these samples are computed with the analytic expression of the Morlet
wavelet.
Secondly, in order to obtain
{
}i=2aT +1
1 t′i ′
′
ψa (ti ) = ψ( ), ti = −aT + (i − 1)
,
a a
i=1
(20)
the discrete version of ψa (t), Zhao et al. adopted the same method with that of cwt function
in Matlab. That is, the discrete samples of ψa (t) are the dilation and translation version of the
samples of ψ(t). For a fixed i of (20), search j of (18) so that
1 ( t′ ) 1
ψa (t′i ) = ψ i ≈ ψ(tj ).
(21)
a
a
a
Insert the expressions of t′i and tj into (21), we have
1 (
i − 1) 1
ψ −T +
≈ ψ(−T + dw (j − 1)).
a
a
a
Then
j=
⌊i−1⌋
a · dw
+1
suffices an approximation for (22)2 . From (23), it is easy to deduce that
′
ti
− t j < dw .
a
(22)
(23)
(24)
By controlling the size of dw , any wanted precision can be obtained for (21) because the Morlet
wavelet function is continuous.
2aT +1
Thirdly, {ha (i)}i=1
is obtained by flipping left and right, taking conjugation3 , taking trans+1
.
lation for {ψa (t′i )}2aT
i=1
The code in Matlab, searching j for a fixed i, is j = [1 + f loor([0 : a ∗ 2T ]/(a ∗ dw ))]
Conjugation operation had been ignored in [17] because only real Morlet wavelet whose wavelet type is 4 was
considered in [17]. The real Morlet wavelet is just the real part of the cmor whose wavelet type is 5
2
3
H. Yi et al. / Journal of Information & Computational Science 11:17 (2014) 6369–6385
6375
Now, the difference between the proposed method and Zhao’s codes is that the discrete samples
of ψea (ti ) are obtained by different means. In Zhao’s code, the discrete samples of the filters are
the dilation and translation version of samples of ψ(t), except the samples of ψ(t) are computed
by analytic expression of Morlet wavelet. In the proposed method, all samples of response function are from the analytic expression of Morlet wavelet, which is the reason of that accuracy is
enhanced.
In Zhao’s method, for each scale a, in order to obtain the {ha (i)}, it is necessary to compute j
for each i using (23). These operations requires extra 2aT + 1 multiplications. In addition, N +
)
2aT + 3(N +2aT
log2 (N + 2aT ) multiplications is needed for the computation of linear convolution.
2
Hence, the total number of multiplications performed for Zhao’s method is N + 4aT + 1 +
3(N +2aT )
log2 (N + 2aT ) for each scale a.
2
3.2
3.2.1
Frequency-domain-based Algorithm
Algorithm Analysis
The wavelet transform defined in (7) can also be written as a frequency integration by applying
the Fourier Parseval formula:
∫ +∞
1 ˆ ˆ
1
Wf (a, b) =
⟨f , ψa,b ⟩ =
fˆ(ω)ψˆ∗ (aω)eibω dω,
(25)
2π
2π −∞
where
ˆ
ψˆa,b (ω) = e−ibω ψ(aω),
(26)
and ψˆ is the Fourier transform of ψ.
In practice, we must deal with the sampled versions of the signal. Denote the sampled signal
by {fd (n) = f (nT ), n = 0, · · · , N − 1}, where T is the sampling period. Discretising (25) yields
M
1
Wf (a, nT ) =
M
2
∑
k=− M
2
+1
( 2π ) ( 2π
) 2π
fˆ
k ψˆ∗
· k · a ei M ·k·n ,
M
M
n = 0, 1, · · · , N − 1.
(27)
In freeware such as Wavelab, M usually takes value N so that the number of {fˆ( 2π
k)} which can
M
be obtained as the discrete Fourier transform of {f (nT )} is just N . As for the computation of
ka)}, we can make use of the analytic expression of
{ψˆ∗ ( 2π
M
σ
ˆ
ψ(ω)
= 2 · e−
2 (ω−η)2
2
(28)
if ψ(t) is defined as (6). Therefore, the CWT can be computed by first finding the Fourier
transforms of the signal and the normalized wavelet. The inverse Fourier transform of the product
k)} and the scaled wavelet {ψˆ∗ ( 2π
ka)} yields one constant-scale slice of the transform.
of {fˆ( 2π
M
M
The main source program of FFT-based algorithm for continuous wavelet transform is given in
Appendix A1.
6376
3.2.2
H. Yi et al. / Journal of Information & Computational Science 11:17 (2014) 6369–6385
The Advantage and Disadvantage of this Method
FFT-based method is popularized due to its low computation complexity and its high precision,
and so it is widely adopted in aforementioned free softwares. However, the method will produce an
artificial periodicity in the wavelet coefficients if signal is not periodic since this method assumes
the data is cyclic if we do nothing about the preprocessing of data. This is demonstrated in the
examples herein, and the method for dealing with this problem is proposed.
For pedagogical reasons, we will concentrate on a Matlab generated data to illustrate this
point. As shown in Fig. 1, the data 1 is a synthetic signal of length 1024. The first half contains
sinusoidal signal superimposition with three different frequencies, and various amplitude, namely
0.8 sin(30πt) + sin(60πt) + 1.2 sin(120πt); The latter half contains the 60 Hz sinusoidal signal with
amplitude 0.6, namely 0.6 sin(120πt). Here the sampling frequency is 400 Hz.
MultiSignal
2.5
2.0
1.5
1.0
0.5
0
−0.5
−1.0
−1.5
−2.0
−2.5
0
200
400
600
800
1000
1200
Fig. 1: A synthetic signal of length 1024. The first half contains sinusoidal signal superimposition with
three different frequencies, namely 0.8 sin(30πt) + sin(60πt) + 1.2 sin(120πt); The latter half contains the
60 Hz sinusoidal signal with amplitude 0.6, namely 0.6 sin(120πt). Here the sampling frequency is 400
Hz
ABSOLUTE−COEFS
The absolute values of the MMWT coefficients of data 1 are presented in Fig. 2. The three
periodic components embedded in the signal are unveiled respectively by their corresponding
MMWT real coefficients in Fig. 3. The outward bulging of the coefficients at the end of the
signal are shown in Fig. 2 and Fig. 3. The line of wavelet coefficients bends to meet artificial
periodic boundary conditions, falsely indicating an increase of amplitude at the end of the signal.
This demonstrates the generic problem when using FFT-based algorithm on non periodic data—
the FFT-based algorithm, in which FFT (fast Fourier transform) is adopted, can induce the
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0
0
200
400
600
Time
800
1000
40 30 20
1200 80 70 60 50 cy/Hz
Frequen
10
Fig. 2: The absolute values of the MMWT coefficients of data 1 computed with FFT-based method
without data preprocessing, where the shape parameter σ 2 = 1 and the center frequency η = 8
6377
H. Yi et al. / Journal of Information & Computational Science 11:17 (2014) 6369–6385
2
0
−2
The MMWT coefficients of 60 Hz
0
200
400
600
800
1000
1200
1000
1200
1000
1200
The MMWT coefficients of 30 Hz
2
0
−2
0
200
400
600
800
The MMWT coefficients of 15 Hz
1
0
−1
0
200
400
600
800
Fig. 3: The MMWT real coefficients of the three periodic components of the signal computed with
FFT-based method without data preprocessing, where the shape parameter σ 2 = 0.8 and the center
frequency η = 10
periodicity into the transform, greatly distorting the information at the end of the signal. This
is a critical problem when the initial and final signal characteristics are very different.
There are several different data preparation methods and here we propose one that is to pad
the end of the time series with zeroes before doing the MMWT and then remove them afterward
[11]. In this study, the time series {fd (n) = f (nT ), n = 0, · · · , N − 1} is padded with sufficient
zeroes so that
{
f (nT ), x = 0, 1, · · · , N − 1
f˜d (n) =
(29)
0,
n = N, N + 1, · · · , M,
where M = 2q for some q ∈ Z and M is the smallest number that is greater than n. It is important
to note that zero padding does not provide any additional information about the spectrum of the
data but results in a “better display” of the Fourier transform of the data [9]. Then
n = 0, 1, · · · , N − 1.
Wfd (a, nT ) = Wf˜d (a, nT ) ,
(30)
Fig. 4 is computed with this method. Comparing Fig. 3 and Fig. 4, we see this method does not
The MMWT coefficients figure of the first period, 60 Hz
2
0
−2
0
200
1000
1200
2
0
−2
0
400
600
800
The MMWT coefficients of the second period, 30 Hz
200
1000
1200
1
0
−1
0
400
600
800
The MMWT coefficients of the third period, 15 Hz
200
1000
1200
400
600
800
Fig. 4: The MMWT real coefficients of the three periodic components of the signal computed with FFTbased method with data preprocessing, where the shape parameter σ 2 = 0.8 and the center frequency
η = 10
6378
H. Yi et al. / Journal of Information & Computational Science 11:17 (2014) 6369–6385
introduce significant distortions of the end regions and eliminates the problem of non periodic
data.
3.2.3
The Computational Complexity for the Frequency-domain-based Method
For the data with length N , pad the end of the data with zeroes so that the length of the data is
M and M ≥ N + 2aT . The FFT of the data requires M2 log2 (M ) multiplications. The spectra of
σ 2 ( 2π ka−η)2
M
2
wavelet , {ψˆ∗ ( 2π
ka) = 2·e−
, k = 0, 1, · · · , M2 , − M2 +1, · · · , −1} can be precomputed and
M
stored in memory. Furthermore, only M2 + 1 discrete spectra need to be stored because the Morlet
wavelet is approximately analytic and therefore ψˆ∗ ( 2π
ka) ≈ 0 for k < 0. Therefore the point-wise
M
M
multiplication of the spectra requires 2 + 1 multiplications. Then the IFFT requires M2 log2 (M )
multiplications. Hence, the total number of multiplications performed is M2 + 1 + M log2 (M ) for
each scale for FFT-based method. In a word, the asymptotic time complexity is O(M log2 (M ))
for all the four methods discussed in this paper.
4
3.3
Comparisons between the Different Algorithms
ABSOLUTE−COEFS
The advantage of convolution-based method (include three methods involved in this paper) compared with FFT-based method is that it will not result in edge effect without previously processing
the data because the algorithm involves linear convolution in time domain instead of FFT, which
assumes the data is cyclic. Fig. 5, computed with the proposed method without preprocessing of
data, unlike Fig. 2, is free of artificial periodic problem.
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0
600 800
1000 1200 80
Time
70
60
50
30
40
Frequency/Hz
10
Fig. 5: Absolute values of MMWT coefficients of data 1 computed with the proposed method without
preprocessing of data, where the shape parameter σ 2 = 1 and the center frequency η = 8
Another merit of the proposed method compared with cwt and Zhao’s is that it accurately
characterizes the amplitude of the three periodic component embedded in the data 1. For the
MMWT real coefficients map of 60 Hz main period shown in Fig. 6, computed with the proposed
method with parameter σ 2 = 0.8, η = 10, the maximum amplitude of MMWT real coefficients
is 1.2 in the first 512 points, the maximum amplitude 0.6 in the second 512 points, while the
corresponding local amplitude of the original signal is respectively 1.2, 0.6. And we can get the
relative errors are 0.2675 × 10−4 and 0.0000 × 10−4 respectively. Unlike the 60 Hz signal, lasting
in all 1024 samples, the 30 Hz, the 2nd period, and the 15 Hz, the 3rd period, shown in Fig. 6,
4
Note the first element in the wavelet spectra corresponds to k = 0 because the first element in the result of
FFT of the data also corresponds to k = 0
6379
H. Yi et al. / Journal of Information & Computational Science 11:17 (2014) 6369–6385
The MMWT real coefficients of 60 Hz
2
0
−2
0
200
400
600
800
1000
1200
1000
1200
1000
1200
The MMWT real coefficients of 30 Hz
2
0
−2
0
200
400
600
800
The MMPS real coefficients of 15 Hz
1
0
−1
0
200
400
600
800
Fig. 6: The real MMWT coefficients of data 1 computed with the proposed method without preprocessing
of data, where the shape parameter σ 2 = 0.8 and the center frequency η = 10
only occur in the first 512 samples, the relative errors between the measured values and the true
values are 0.2589 × 10−4 and 0.0000 × 10−4 respectively. If computed with FFT-based method
with preprocessing of data, the relative errors are the same with the proposed method. When
computed with another parameters, the relative errors for different methods are displayed in
Table 1. At the same time, it should be mentioned that the errors in column “proposed method”
are identical with column “FFT based method” in Table 1, and this is not typographical errors.
This phenomenon indicates some equivalences between these two methods. In fact, both of the
two methods can supply error-free linear convolution except for round-off error. In addition, the
number of multiplications for FFT-based method is less than that of other methods.
As another example illustrating the different performance of these methods, we recompute the
scalograms and wavelet ridges of two hyperbolic chirps shown in Figure 4.14 of [6], with the four
different continuous wavelet transform methods, and the result is shown in Fig. 7. It is obvious
that the proposed method and FFT-based method outperform the other two, while the precision
of Zhao’s is slightly better than that of cwt method.
Table 1: The relative errors between the amplitude of the MMWT real coefficients and the corresponding
amplitude of periodic components embedded in the data 1 with different methods
σ2 = 1
Proposed method
cwt method
Zhao’s method
FFT-based method
relative error1
0.2196 × 10−5
3.254 × 10−2
0.9530 × 10−2
0.2196 × 10−5
relative error2
0.0000 × 10−5
3.652 × 10−2
0.3688 × 10−2
0.0000 × 10−5
relative error3
0.2108 × 10−5
0.3960 × 10−2
1.722 × 10−2
0.2108 × 10−5
relative error4
0.0001 × 10−5
0.6820 × 10−2
2.971 × 10−2
0.0001 × 10−5
Proposed method
cwt method
Zhao’s method
FFT-based method
relative error1
0.0042 × 10−6
3.305 × 10−2
0.8670 × 10−2
0.0042 × 10−6
relative error2
0.0000 × 10
−6
−2
3.678 × 10
−2
0.2115 × 10
0.0000 × 10−6
relative error3
0.0041 × 10−6
0.3358 × 10−2
0.8435 × 10−2
0.0041 × 10−6
relative error4
0.7074 × 10−6
0.6217 × 10−2
0.4829 × 10−2
0.7074 × 10−6
η = 10
σ 2 = 1.5
η = 10
6380
H. Yi et al. / Journal of Information & Computational Science 11:17 (2014) 6369–6385
400
200
Figure based on cwt
Frequency
Frequency
Figure based on cwt
400
200
200
400
200
0
0 100 200 300 400 500 600 700 800 900 1000
200
0 100 200 300 400 500 600 700 800 900 1000
Figure based on the proposed method
400
200
0
Frequency
Frequency
0
0 100 200 300 400 500 600 700 800 900 1000
Figure based on FFT−based algorithm
0 100 200 300 400 500 600 700 800 900 1000
Figure based on the algorithm of ZHAO
400
0
Frequency
Frequency
0
0 100 200 300 400 500 600 700 800 900 1000
Figure based on the proposed method
400
200
0
Frequency
Frequency
0
0 100 200 300 400 500 600 700 800 900 1000
Figure based on the algorithm of ZHAO
400
0 100 200 300 400 500 600 700 800 900 1000
Figure based on FFT−based algorithm
400
200
0
0 100 200 300 400 500 600 700 800 900 1000
Fig. 7: The scalograms and wavelet ridges of two hyperbolic chirps shown in Figure 4.14 of [6] are
calculated with four different continuous wavelet transform methods, where the Morlet wavelet with the
shape parameter σ 2 = 1 and the center frequency η = 5 is used
4
Application: The Necessity of Suitable Method
The different analytical effect of temperature data with different algorithms are presented in this
section, which illustrates the necessity of choosing suitable algorithm. The temperature data
are from Jan. 1, 1951 to Dec. 31, 2010, and the experimental site is No. 50978 weather station,
Heilongjiang, China.
Fig. 8 and Fig.9 are contours of real part of MMWT of summer average temperature data for
50978 station, the two figures are computed with same parameters except for the continuous
wavelet transform algorithms, the Fig. 8 with Zhao’s method, the Fig. 9 with proposed method
(the corresponding figures computed by cwt method and FFT-based are not shown here for lack
of space but are available from the authors upon request). In these figures, the coefficients with
negative values, which mean that the temperature is in cold stage, are plotted with dotted, blue
line, while the coefficients with positive values, which mean that the temperature is in warm
stage, are plotted with dashed, red line. And so the negative and positive oscillation of local time
structures of summer temperature are visualized clearly from both figures.
However, the contour in Fig. 9 are more distinct than those in Fig. 8, which means that there
are less peaks and valleys in Fig. 9 compared with those in Fig. 8.
6381
H. Yi et al. / Journal of Information & Computational Science 11:17 (2014) 6369–6385
60
Period/year
50
28
17.9594
10.7561
5.3781
2.689
0
1951
1960
1970
1980
Time/year
1990
2000
2010
Fig. 8: The contour of real part of wavelet coefficients of summer average temperature of 50978 station
computed with Zhao’s method
60
Period/year
50
28
17.8302
10.9126
5.3781
2.6125
0
1951
1960
1970
1980
Time/year
1990
2000
2010
Fig. 9: The contour of real part of wavelet coefficients of summer average temperature of 50978 station
computed with the proposed method
The above phenomena, which are perceived by our visual sensors, can also be illustrated from
another quantitative way. That is, the proposed method can determine the multi-period of
the data more definite than Zhao’s methods. The multi-period, which concentrate the essential
dynamic of data, shown in Table 2, are calculated by using wavelet variance [3, 4, 11]. The
performance of proposed method superior to that of Zhao’s and FFT-based can be seen from
two aspects. Firstly, the number of periods computed by Zhao’s and by FFT-based outnumber
that of proposed method. With a few periods, we can give a more accurate characterization of
6382
H. Yi et al. / Journal of Information & Computational Science 11:17 (2014) 6369–6385
data. Secondly, the period computed by proposed method are almost the same when the center
frequency η, the wavelet parameter, changes. On the contrary, the number of period calculated
by other methods will change when the wavelet parameter changes, which may lead to confusion
when we determine what period exist in the data.
Table 2: The multi-period (unit: year) of 60 year summer temperature data differentiated by different
continuous wavelet transform methods
σ 2 = 1, η = 5
σ 2 = 1, η = 6
Zhao’s
2.1189, 2.2776, 2.4305, 2.5750
2.1497, 2.3956, 2.5381, 2.6697, 2.8081, 2.9324
method
2.6890, 5.3781, 10.7561, 17.9594
5.4170, 11.1515, 17.8302, 44.6060, 51.9091, 57.4307
proposed
2.6125, 5.3781, 10.9126, 17.8302
2.5565, 5.4170, 10.9916, 17.8302
F-based
2.6505, 5.3394, 10.2260, 16.2327
2.5750, 5.3781, 10.0794, 16.2327, 21.0512, 31.7698
The above analysis indicates that the wavelet transform coefficients computed by the proposed method are more sparse than those of all the other methods. On the other hand, sparsity of wavelet coefficients may be measured with the entropy of those wavelet coefficients
[5]. The entropy here is termed wavelet entropy. In order to define entropy, we rearrange the
{|Wf (ai , bj )|2 }{i,j} computed by (7) as {ck }k=1,...,M , then normalize it to obtain:
ck
dk = ∑M
k=1 ck
.
(31)
The wavelet entropy is calculated as
En = −
M
∑
dk log dk ,
(32)
k=1
with the convention 0 log 0 = 0 by definition. Next, from the Table 3, we see the wavelet entropy
of the proposed method is the smallest of all the others, which singles out the proposed method
as the superior method to the others.
Table 3: The wavelet entropy corresponding to different wavelet transform methods
σ 2 = 1, η = 5
2
σ = 1, η = 6
5
cwt
Zhao’s
FFT-based
Proposed
9.3327
9.1273
9.2079
9.1263
9.3318
9.1300
9.1985
9.1291
Conclusion
The wavelet transform is basically a linear filter whose response function is given by the wavelet
function. So the algorithms of continuous wavelet transform consist of two methods, frequencydomain calculation and time-domain calculation. FFT-based method, the typical representation
of frequency-domain computation, with its low computation complexity and its high precision, is
H. Yi et al. / Journal of Information & Computational Science 11:17 (2014) 6369–6385
6383
widely adopted in free softwares. But, M in (27) taking value N , the length of data, implicitly
assumes the data is cyclic, and so the artificial periodic problem emerges. Meanwhile, different
values of M implies different data preprocessing method.
The proposed method, a representation of time-domain computation, is free of artificial periodic
problem. And the computation precision of the proposed method and FFT-based method is
superior to other existing methods because the two methods calculate the response function by
the analytic expression of wavelet, the proposed method in time-domain, and the FFT-based
method in frequency-domain.
The proposed method, with the most sparse coefficients, is the most suitable method to analyzing the multi-period structure of temperature data.
The asymptotic time complexity is O(M log2 (M )) for all the four methods discussed in this
paper.
It should be noted that the proposed method can be adapted to other wavelets having analytic
expressions although this paper exclusively discusses Morlet wavelets. And it is not surprising
that the precision of the proposed method is superior to that of cwt function in wavelet toolbox of
Matlab and Zhao’s considering the fact that cwt and Zhao’s can deal with those wavelets without
analytic expressions such as Daubechies wavelets.
Appendix
A1: The Main Source Program of FFT-based Algorithm for MMWT
The main source program of FFT-based algorithm is given as follows:
n= l e n g t h ( S i g ) ; % t h e number o f s i g n a l
f f t x = f f t ( S i g ( : ) ) ; % F o u r i e r Transform o f s i g n a l S i g
L e n S c a l e s=l e n g t h ( S c a l e s ) ; % t h e number o f S c a l e s
omega = [ ( 0 : ( n / 2 ) ) ((( −n /2)+1): −1) ] . ∗ ( 2 ∗ p i /n ) ;
Coefs = z e r o s ( LenScales , n ) ;
f o r k = 1 : L e n S c a l e s % computation o f c o e f f i c i e n t s f o r each s c a l e
a= S c a l e s ( k );% a i s t h e s c a l e v a l u e
omega = omega ( : ) ;
P s i = 2∗ exp(−sigma2 / 2 ∗ ( a . ∗ omega − e t a ) . ˆ 2 ) ; % sigma2
i s th e shape parameter , and e t a i s t h e c e n t e r f r e q u e n c y .
P s i = P s i . ∗ ( omega>0);% P s i ( omega )=0 , when omega<0
C o e f s =(k , : ) = i f f t ( f f t x . ∗ P s i ) ;
end
A2: The Main Source Program of the Proposed Method for MMWT
The main source program(in this time, the Morlet wavelet is given as cmor function defined in
(8) in stead of (6)) of proposed method is given as follows:
[
√
√
√ ]
wsupport= 8 √F2 b ; % the effective support of Morlet defined in (8) is − 8 √F2 b , 8 √F2 b
6384
H. Yi et al. / Journal of Information & Computational Science 11:17 (2014) 6369–6385
L e n S c a l e s=l e n g t h ( S c a l e s ) ;
% t h e number o f S c a l e s
lenSig = length ( Sig ) ;
% t h e number o f s i g n a l
C o e f s = z e r o s ( L e n S c a l e s , l e n S i g );% t h e c o e f f i c i e n t s matrix
f o r k = 1: LenScales
% computation o f c o e f f i c i e n t s f o r each s c a l e
a = Scales (k );
% a i s t he s c a l e v a l u e
t = −round ( a∗ wsupport ) : 1 : round ( a∗ wsupport ) ;
% t h e e f f e c t i v e s u p p o r t o f Morlet with s a l c e a i s
[−a∗ wsupport , a∗ wsupport ] ,
and t h e sample f r e q u e n c y o f w a v e l e t i s 1 Hz
wvalue= f l i p l r ( 2 ∗ ( p i ∗Fb)ˆ( −1/2)∗ exp ( 1 i ∗2∗ p i ∗Fc . ∗ t /a ) .
∗ exp(− t . ˆ 2 / ( Fb∗a ˆ 2 ) ) ) ;
wvalue=c o n j ( wvalue ) ;
%computation o f t he w a v e l e t f u n c t i o n with t he a n a l y t i c e x p r e s s i o n .
temp = conv ( Sig , wvalue ) / a;% computation o f c o n v o l u t i o n
C o e f s ( k , :)= temp ( round ( a∗ wsupport )+1:
l e n g t h ( temp)−round ( a∗ wsupport ) ) ;
end
A3: The Main Source Program of Zhao’s Method for MWT
The main source program of Zhao’s method is given as follows:
signal = signal ( : ) ’ ; len = length ( signal ) ; nbscales =
l e n g t h ( s c a l e s ) ; c o e f s = z e r o s ( n b s c a l e s , l e n ) ; f o r k=1: n b s c a l e s
a=s c a l e s ( k ) ;
[ p s i , x v a l ] = wavefun (wname ) ;
x v a l = xval−x v a l ( 1 ) ;
dx = x v a l ( 2 ) ;
xmax = x v a l ( end ) ;
j =[1+ f l o o r ( [ 0 : a∗xmax ] / ( a∗dx ) ) ] ;
i f l e n g t h ( j )==1, j =[1 1 ] ; end
f= f l i p l r ( p s i ( j ) ) ;
f=c o n j ( f ) ;
c o e f s ( k , : ) = wkeep ( conv ( s i g n a l , f ) , l e n ) / s q r t ( a ) ; end
References
[1]
[2]
[3]
[4]
Ren´e Carmona, Wen-liang Hwang, Bruno Torresani, Practical Time-Frequency Analysis: Gabor
and Wavelet Transforms, with an Implementation in S, Volume 9, Academic Press, 1998
M. Issac, G. Renuka, C. Venugopal, Wavelet analysis of long period oscillations in geomagnetic
field over the magnetic equator, Journal of Atmospheric and Solar-Terrestrial Physics, 66(11),
2004, 919-925
David Labat, Recent advances in wavelet analyses: Part 1 - A review of concepts, Journal of
Hydrology, 314(1), 2005, 275-288
David Labat, Josyane Ronchail, Jean Loup Guyot, Recent advances in wavelet analyses: Part
2 - Amazon, parana, orinoco and congo discharges time scale variability, Journal of Hydrology,
314(1), 2005, 289-311
H. Yi et al. / Journal of Information & Computational Science 11:17 (2014) 6369–6385
6385
[5]
J. Lin, L. Qu, Feature extraction based on morlet wavelet and its application for mechanical fault
diagnosis, Journal of Sound and Vibration, 234(1), 2000, 135-148
[6]
Stephane Mallat, A Wavelet Tour of Signal Processing: The Sparse Way, Access Online via Elsevier, 2008
[7]
Jorge Martinez, Richard Heusdens, Richard C. Hendriks, A generalized poisson summation formula
and its application to fast linear convolution, Signal Processing Letters, IEEE, 18(9), 2011, 501-504
[8]
Steven D. Meyers, Brian G. Kelly, James J. O’Brien, An introduction to wavelet analysis in oceanography and meteorology: With application to the dispersion of yanai waves, Monthly Weather
Review, 121(10), 1993, 2858-2866
[9]
John G. Proakis, Digital Signal Processing: Principles, Algorithms, and Applications, 4th Ed.,
Pearson Education India, 2007
[10] Ivan W. Selesnick, C. Sidney Burrus, Fast Convolution and Filtering, The Digital Signal Processing
Handbook, 1999, 8-2
[11] C. Torrence, G. P. Compo, A practical guide to wavelet analysis, Bulletin of the American Meteorological Society, 79(1), 1998, 61-78
[12] Michael Unser, Akram Aldroubi, Steven J. Schiff, Fast implementation of the continuous wavelet
transform with integer scales, IEEE Transactions on Signal Processing, 42(12), 1994, 3519-3523
[13] R. Werner, The latitudinal ozone variability study using wavelet analysis, Journal of Atmospheric
and Solar-Terrestrial Physics, 70(2), 2008, 261-267
[14] H. Yi, Q. Fan, An algorithm for the determination of multi-period structure of time series, IEEE
International Conference on Information and Automation (ICIA), 2010, 1684-1689
[15] H. Yi, H. Shu, The improvement of the Morlet wavelet for multi-period analysis of climate data,
Comptes Rendus Geoscience, 344(10), 2012, 483-497
[16] H. Yi, H. Shu, T. Zhang, Q. Fan, Applications of Morlet wavelets in time-frequency localization
of signals, Mathematica Applicata, 23(2), 2010, 395-400
[17] Y. Zhao, X. Yuan, Y. Wei, Realization of continuous wavelet transform of sequences by matlab,
Journal of Sichuan University (Natural Science Edition), 43(2), 2006, 325-329