A Class of Algorithms for Continuous Wavelet Transform Based on the Circulant Matrix

: The Continuous Wavelet Transform (CWT) is an important mathematical tool in signal processing, which is a linear time-invariant operator with causality and stability for a ﬁxed scale and real-life application. A novel and simple proof of the FFT-based fast method of linear convolution is presented by exploiting the structures of circulant matrix. After introducing Equivalent Condition of Time-domain and Frequency-domain Algorithms of CWT, a class of algorithms for continuous wavelet transform are proposed and analyzed in this paper, which can cover the algorithms in JLAB and WaveLab, as well as the other existing methods such as the cwt function in the toolbox of MATLAB. In this framework, two theoretical issues for the computation of CWT are analyzed. Firstly, edge effect is easily handled by using Equivalent Condition of Time-domain and Frequency-domain Algorithms of CWT and higher precision is expected. Secondly, due to the fact that linear convolution expands the support of the signal, which parts of the linear convolution are just the coefﬁcients of CWT is analyzed by exploring the relationship of the ﬁlters of Frequency-domain and Time-domain algorithms, and some generalizations are given. Numerical experiments are presented to further demonstrate our analyses.


Introduction
In recent years, different Time-frequency representations, such as, empirical mode decomposition [1], wavelet transform [2] and its variants, empirical wavelet transform [3,4], synchrosqueezed wavelet transforms [5], have been used for analyzing nonlinear and non-stationary signals.To name only a few, the method of fused empirical mode decomposition and wavelets is applied to detection-location of damage in a truss-type structure [6]; Wavelet transform is used for the pattern recognition for diagnosis, condition monitoring and fault detection [7][8][9].It is also used to design an algorithm for Brain-computer interfacing [10], to detect the exact onset of chipping of the cutting tool from the workpiece profile [11], to determine the length of piles [12]; Synchrosqueezed wavelet transform is used for global and local health condition assessment of structures [13], for modal parameters identification of smart civil structures [14].
This paper gives priority to the computation of continuous wavelet transform (CWT) for Morlet-type wavelets.However, for this type of wavelets, it is not possible to use a multiresolution framework for the computation of CWT [15].Michael Unser [16] used Exponentials or B-spline window to approximate Gabor window, which can achieve O(N) complexity per scale.Another kind of method for the computation of CWT with Morlet-type wavelets concerns directly discretizing the integral expression of Wavelet transform.These methods allow us to use arbitrary values for the scale variable, and require the explicit expression of wavelet function and cannot deal with those wavelets without analytic expressions such as Daubechies wavelets [17].This kind of method include Time-domain methods [17] and Frequency-domain methods [18].
Frequency-domain methods for CWT are also widely used in the freewares such as JLAB (available online http://www.jmlilly.net),Wavelab (available online http://www-stat.stanford.edu/$\sim$wavelab), which can achieve satisfactory precision while the edge effect, defined in [18], may occur and the complexity of those methods is O(N log(N)) where N is the length of the signal.
The function cwt in wavelet toolbox of MATLAB, which has poorer precision than Frequency-domain methods mentioned above, is a representation of the Time-domain method [17].For the details of the comparison of the different methods, we refer to [17].
The computation of linear convolution can be realized by using the circular convolution, while the computation of circular convolution can be computed by using the FFT-based fast method [19].In this paper, a novel and simple method is given to prove the FFT-based fast method of linear convolution by exploiting the structures of circulant matrix.After introducing Equivalent Condition of Time-domain and Frequency-domain Algorithms of CWT, a class of algorithms for continuous wavelet transform are proposed and analyzed, which can cover the algorithms in JLAB and WaveLab, as well as the other existing methods such as the cwt function in the toolbox of MATLAB.In this framework, two theoretical issues for the computation of CWT are analyzed.
Firstly, by using Equivalent Condition of Time-domain and Frequency-domain Algorithms of CWT to design Frequency-domain Algorithm, the edge effect, defined in [18], can be avoided.To be specific, in [18], the time series is padded with sufficient zeroes to bring the total length N up to the next-higher power of two, thus limiting the edge effects and speeding up the Fourier transform.In fact, the same number of zeros is padded for all scales in [18], which may cause some troubles, for example, increasing the amount of data to process, and consequently the computational complexity.In this article, the time series (the data) is padded with 2aT zeros while a is the scale.These two zero padding methods are compared in Remark 2 of this paper.
Secondly, due to the fact that linear convolution expands the support of the signal, which parts of the linear convolution are just the coefficients of CWT is analyzed by exploring the relationship of the filters of Frequency-domain and Time-domain algorithms (see Theorem 4), and some generalizations are given (see Theorem 5).
This paper is organized as follows.Section 2 gives some definitions and theorems concerning circulant matrix and linear convolution.Section 3 analyzes algorithms of continuous wavelet transform.Section 4 presents numerical experiments to demonstrate our results and finally, we end this paper with conclusions and discussions in Section 5.

Primary Definitions and Theorems
Definition 1. Circulant matrix: An N × N circulant matrix C takes the form [20,21] A circulant matrix is fully specified by one vector, c, which appears as the first column of C.

Definition 2. Discrete Fourier transform(DFT):
N is the N-th root of unity.
Equation ( 2) can be written as where F N is a N × N matrix, defined as It is easy to verify that where " * " means conjugate transposition.Then, the inverse discrete Fourier transform (IDFT) is given as Theorem 1. References [20,21] The matrix C defined in (1) can be diagonalized by the DFT matrix F N , namely, where c is the first column of C, i.e., c = [ Proof.It is easy to verify that the normalized eigenvectors of C are given by where ω j = e 2πij N with i = √ −1 is the N-th root of unity.The corresponding eigenvalues are then given by From (2), ( 3) and ( 9), we have Using {v j } N−1 j=0 defined in Equation ( 8) as column vectors to form a then we have Furthermore, we have 4), ( 5) and (11).Then we have Definition 3. Linear convolution [23]: A time-invariant linear operator L can be represented as a linear convolution.To be specific, we denote by δ[n] the discrete Dirac Any signal f [n] can be decomposed as a sum of shifted Diracs: Let Lδ[n] = h[n] be the discrete impulse response.Linearity and time invariance implies that where " " represents linear convolution.

Definition 4. Causality and Stability [23]:
One can verify that the filter is stable if and only if h ∈ l 1 (Z).
Proposition 1. Assume that the length of f and h are finite.To be specific, If f = [ f 0 , f 1 , . . ., f n−β−1 ] T , and h = [h 0 , h 1 , . . ., h β ] T , then the linear convolution of f and h defined in ( 14) can be written as where H ∈ R n×(n−β) , and Proof.We refer the reader to [19] for details.
Suppose we wish to compute the polynomial product c(x) = a(x) • b(x), the ordinary product expression for the coefficients of c(x) involves a linear convolution.

Definition 5. Circular convolution [19]:
The circular convolution of a signal {x k=0 is defined as a matrix vector multiplication as follows. y where C is defined in equation ( 1), and " " represents circular convolution.
Proposition 2. The computation of circular convolution can be realized by using FFT-based fast method.
In other word, Equation ( 17) can be written as Proof.The matrix involved in ( 17) is a circulant matrix, thus ( 18) is obtained by using Theorem 1.
The equivalent condition of circular convolution and linear convolution will be given in the following theorem.
then linear convolution g = f h of f and h defined in ( 14) can be computed as circular convolution f h of f and h, namely, where F N , F N −1 are defined by ( 4), ( 5) respectively; ". * "means componentwise product of two vectors.
Proof.From Proposition (1), we get g = H f , where H ∈ R n×(n−β) is defined in (16).Now extend H to a circulant matrix [21] thus the first column of C is h.Therefore, by using Definition 5 and Proposition 2.
From Theorem 2, the algorithm of linear convolution of g = f h is given as follows (See Algorithm 1).
Algorithm 1: FFT-based fast method for linear convolution. Input: Appending an array of zeros to the end of h to obtain h so that the length of h is n.
3. Appending an array of zeros to the end of f to obtain f so that the length of f is n Additionally, the asymptotic time complexity is O(nlog 2 (n)) for Algorithm 1.Note the condition n = length( f ) + length(h) − 1 plays a pivotal role in the equivalence of linear convolution and circular convolution.This condition will also be used in the following part of this paper.

Definition 6. Equivalent Condition of Linear Convolution and Circular Convolution
Equivalent Condition of Linear Convolution and Circular Convolution.

A Class of Algorithms for Continuous Wavelet Transform
where ψ(t) is a continuous function called the mother wavelet and the overline represents the operation of complex conjugate.
For real-life applications, the length of f (t) is finite and the support of ψ(t) is compact.Without loss of generality, assume the sampling rates of the signal and the wavelet are equal to 1, then Equation ( 21) can be written as , where m ∈ Z.
Assume the sampled signal is { f [i]} N−1 i=0 and the support of the mother wavelet then the support of ψa (t) is [0, 2aT].Now, we can get the following result.
Theorem 3. Assume the length of f (t) is finite and the support of ψ(t) is compact, and the sampling rates of the signal and the wavelet are equal to 1, then CWT of { f [i]} N−1 i=0 can be written as where we assume aT is an integer for the convenience of analysis.Furthermore, the corresponding filter is causal and stable.
Proof.By using Proposition 1 and Equations ( 23) and ( 24) is deduced from (22).Furthermore, define then the discrete filter {h a [i]} i∈Z is causal and stable since h a [i] = ψa (i) = 0 if i < 0 and h a ∈ l 2 (Z).Furthermore, Equation ( 24) can be written as can be implemented with linear convolution [17].From Algorithm 1, we can get the algorithm of CWT (see Algorithm 2).
Algorithm 2: Time-domain fast algorithm for CWT.
, where the overline represents the operation of complex conjugate;

Frequency-Domain Algorithm for CWT
The algorithm can be further optimized by taking advantage of the analytic expression of ψ(ω), the Fourier transform of the mother wavelet ψ(t).In fact, the wavelet transform defined in (21) can also be written as a frequency integration by applying the Fourier Parseval formula [23].
Assume the sampled signal is { f [n]} N−1 n=0 .Discretizing (27) and considering the periodic property of discrete Fourier transform yields It is seen by comparing (28) and ( 18) that (28) represents a circular convolution.By considering Equivalent Condition of Linear Convolution and Circular Convolution, It is reasonable to let In this case, the values of { f ( 2π M k)} M−1 k=0 can be computed as the discrete Fourier transform of As for the computation of { ψ( 2π M ka)} M−1 k=0 , we can make use of the analytic expression of ψ(ω).For example, the Morlet wavelet is defined as [24] where σ 2 is shape parameter, and η is center frequency.Furthermore, the Morlet wavelet is approximately analytic and therefore ψ( 2π M ka) ≈ 0 for k < 0. So by the periodic property of discrete Fourier transform.Equation (32) will be used to design algorithm of CWT (see the 4th-6th steps of Algorithms 3).29) is called Equivalent Condition of Time-domain and Frequency-domain Algorithms of CWT.

Remark 1.
In freeware such as Wavelab, the parameter M in Equation (29) takes value N, where N is the length of signal.In this case, { f ( 2π M k)} M−1 k=0 can be obtained as the discrete Fourier transform of the data.Furthermore, the length of the result of ifft is just N.However, the method will produce artificial periodicity, which is called edge effect by Torrence [18], in the wavelet coefficients if signal is not periodic.In order to limit the edge effect, in [18], the data is padded with sufficient zeros before doing the CWT, then the first N coefficients of the corresponding convolution are just the coefficients of the CWT.Nevertheless, the reason for this is not answered in previous papers in the author's knowledge.Put another way.Why the last step of Algorithm 3 is dissimilar from the last step of Algorithm 2? The answer will be given in the following Theorem 4.
Lemma 1.The filter used in Equation ( 27) is Proof.We refer the reader to [23] for details.
Theorem 4. The wavelet coefficients of the CWT of { f [n]} N−1 n=0 is the first N coefficients of (28).While the length of the total coefficients of (28) is M.
Theorem 5. Assume that (Without loss of generality, aT is assumed to be an integer.) Proof.If m = 0, Equation ( 34) is simply Equation ( 28), the conclusion is deduced from Theorem 4. Since , the conclusion for the m = 2 case can also be proved in the same way.Now, Frequency-domain algorithm of CWT with Morlet wavelet as the mother wavelet is given as follows.
Algorithm 3: Frequency-domain fast algorithm for CWT (In real application, the case m = 0 is enough for the calculation of CWT). Input: , the Fourier transform of the Morlet wavelet ψ(t) T, where [−T, T] is the support of ψ(t).Output: 2. Define ϕ m (ω) = e −imaTω ψ(aω), m = 0, 1, 2, where the overline represents the operation of complex conjugate; 3. Let M = 2aT + N; Remark 2. Note that the previous data preparation method takes M = 2 q , where q = min n∈Z {2 n > N} [17,18].However, this method may fail for some real life data.We propose M = N + 2aT, Equivalent Condition of Time-domain and Frequency-domain Algorithms of CWT, then the edge effect, defined in [18], that may occur in JLAB, Wavelab and [25] can be avoided.See Figures 6 and 7 for details.

Numerical Experiments
The experiments are conducted on two types of data, one for synthetic data, another for real-life data.Entropy can be used to measure the sparsity of wavelet coefficients [17].In order to define entropy, {|W f (a i , b j )| 2 } i,j are rearranged as {c k } k=1,...,M , then are normalized to obtain: The wavelet entropy is calculated as with the convention 0 log 0 = 0 by definition.Experiments for Data 1: Data 1 is a synthetic signal with length N. The first half contains sinusoidal signal superimposition with three different frequencies,namely a 4 sin(30πt) + a 3 sin(60πt) + a 1 sin(120πt); The latter half contains the 60 Hz sinusoidal signal with amplitude a 2 , namely a 2 sin(120πt), where a 4 = 1, a 3 = 1.2, a 1 = 1.4, a 2 = 0.6, and the sampling frequency is 400 Hz.Data 1 with length N = 1024 is presented in Figure 1a. Figure 1b, the absolute values of CWT coefficients of Figure 1a, computed with the awt function in Wavelab, where the shape parameter σ 2 = 1 and the center frequency η = 8, manifests edge effects.Figure 1c without edge effects is computed by using the case 0 of Algorithm 3. Table 1 gives the computational results of Data 1 with wavelet parameter σ 2 = 1, η = 8 by using different methods."cwt" means the cwt function in wavelet toolbox of Matlab."Zhao's method" is an improved version of cwt [26]."Direct" means the computation method of linear convolution by using Equation (15).If the length of signal is far larger than the length of the filter, an "Overlap-add" procedure for the calculation of linear convolution is faster than Direct method and Algorithm 1 [23].The "wavelet entropy" shown in Table 1 can measure the sparsity of wavelet coefficients [17]."Err_i" means the relative error of a i and the corresponding maximum amplitude of CWT coefficients.
From Table 1, we know that the precision, wavelet entropy of "Direct", Algorithms 2 and 3, "Overlap-add" are almost the same and are more optimal to cwt function in toolbox of Matlab.In fact, the precision of two methods in [17], "FFT based method" and the proposed method is almost the same.This phenomenon indicates some equivalence of these two methods.As a matter of fact, these two methods can be categorized as Algorithms 2 and 3 respectively.It is seen from Table 2 that the wavelet entropy for Algorithm 2 is the smallest of all the methods.Therefore, Algorithm 2 is the most suitable method for this temperature data.
Figures 6 and 7 are calculated with different data preparation methods.To be specific, Figure 6 is computed with the frequency-domain method of CWT with data preparation method given in previously published papers, namely, M = 2 q , q = min n∈Z {2 n > N} while Figure 7 is computed with the case 0 of Algorithm 3 with M = N + 2aT.Note the different structures of contours in the upper edge of these two figures.The contour is not closed in the upper edge of Figure 6, which means that there will be a valley and a peak in the left and right part of this edge.So the edge effect is obvious in Figure 6.   .This contour, computed with frequency-domain method of CWT with data preparation method given in previously published papers, namely, M = 2 q , q = min n∈Z {2 n > N}, manifests edge effect.

Conclusions
The continuous wavelet transform of a signal with finite length for a fixed scale is considered as a linear time-invariant operator.Furthermore, the filter with causal and stability is constructed.The algorithm of linear convolution constitutes a unifying framework to the continuous wavelet transform methods previously published in [17].The precision of these methods based on this framework is almost the same, no matter what method is used, and higher than other methods that use the approximation of wavelet, for example, the cwt function in wavelet toolbox of Matlab.The edge effect is also easily handled by using Equivalent Condition of Time-domain and Frequency-domain Algorithms of CWT.
The algorithms of CWT consist of two methods, time-domain method and frequency-domain method.For time-domain method, by constructing the causal filter ψa (t), we know the wavelet coefficients are just the middle part of the corresponding linear convolution.For frequency-domain method, by exploring the relationship of ψa (t) and ψa (t), we know the wavelet coefficients are just the first N coefficients of the corresponding convolution.Furthermore, by constructing the different filters, the wavelet coefficients can be the first N coefficients, or the middle N coefficients, or the last N coefficients of the corresponding convolution for frequency-domain method.
There are three methods for the calculation of linear convolution.The first one is to directly implement the definition of linear convolution.The second one is known as the FFT-based method, for example, Algorithm 1. Lastly, if the length of the signal is far larger than the length of the filter, an overlap-add procedure for the calculation of linear convolution is faster than the previous two methods [23].How to combine the Frequency-domain method and an overlap-add procedure for the computation of CWT is a question for future research.

3. 1 .
Time-Domain Algorithm for CWT Definition 7. Continuous wavelet transform: The continuous wavelet transform of a function f (t) at a scale a ∈ R + and translational value b ∈ R is expressed by the following integral [17] W f (a, b)

Figure 1 .
Figure 1.(a) is Data 1 with N = 1024.(b) the absolute values of CWT coefficients of (a), computed with the awt function in Wavelab, where the shape parameter σ 2 = 1 and the center frequency η = 8, manifests edge effects.(c) without edge effects is computed with Algorithm 3.

Figure 2
is computed by using the case 1 of Algorithm 3. The wavelet coefficients of Figure2bcan exactly characterize the time-frequency local properties of data of Figure1a.At the same time, Figure2a,c fail to do so.In fact, it is known from the case 1 of Algorithm 3, the wavelet coefficients should be chosen as the middle part of the coefficients of the convolution.However, the wavelet coefficients of Figure2a,c are respectively chosen as the first N and the last N part of the coefficients of the convolution.

Figure 2 .
Figure 2.Figure 2 is computed by using the case 1 of Algorithm 3. The wavelet coefficients of (b) can exactly characterize the time-frequency local properties of data of (a).At the same time, (a,c) fail to do so.

Figure 2 Experiments for Data 2 :
Figure 2.Figure 2 is computed by using the case 1 of Algorithm 3. The wavelet coefficients of (b) can exactly characterize the time-frequency local properties of data of (a).At the same time, (a,c) fail to do so.

Figure 5 .Figure 6
Figure 5.This contour is computed by using CWT function in wavelet toolbox of Matlab with the Morlet wavelet parameter σ 2 = 1, η = 5 for Data 2.

Figure 7 .
Figure 7.This contour, computed with the case 0 of Algorithm 3 with M = N + 2aT, manifests no edge effect.

Table 1 .
Comparisons between the different CWT methods for Data 1.

Table 2 .
The wavelet entropy corresponding to different wavelet transform methods.