The Adaptive Analysis of Shock Signals on the Basis of Improved Morlet Wavelet Clusters

Abstract

Morlet wavelets do not satisfy the permissibility condition of wavelet analysis, and there are therefore no inverse transformations for Morlet wavelet transforms. In this paper, we put forward the Yang and Pan transform (YPT), which is an adaptive discrete analysis method for shock signals. First, we improved the Morlet wavelet so that the centre and radius of the frequency window can be easily adjusted in the frequency domain. Second, we proposed the extremum frequency concept and analysed the extremum situation of the improved Morlet wavelet. Third, combining the improved Morlet wavelet and extremum frequency, we advanced the theory of the YPT, which does not need to satisfy the permissibility condition. We then continued by using a smoothing operator that can smooth the potentially distorted signal reconstructed after being analysed by the YPT and filtered by using the threshold filtering theory. This operator proved to be simple and efficient. Finally, a noisy signal was reconstructed after being analysed and filtered using the YPT and threshold filtering, respectively, to verify the validity of the theory, and the YPT was compared with the discrete wavelet transform (DWT). As a supplement to the theory in engineering, the shock signals about a gun automatic mechanism were also analysed using the theory in this paper. Good results were obtained, thereby demonstrating that the YPT can be helpful to further extract the features of shock signals in pattern recognition and fault diagnosis.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Wavelet transforms are a significant method for analysing signals, and Morlet wavelets are often used as the kernel functions of wavelet transforms. In 1982, Morlet first used them to analyse seismic signals, employing a complex wavelet using a Gauss envelope and a special case of Gabor wavelets as well [1]. The continuous wavelet transform of a signal by Morlet wavelets can achieve arbitrary high resolution in the time or frequency domains [2, 3]. There is a 90-degree phase shift between its real and imaginary parts, which makes it easy to obtain the instantaneous frequency and phase of the signal, and it therefore has a wide range of applications. For example, Morlet wavelets have been applied to signal filtering and denoising [4], mechanical fault diagnosis [5–8], the analysis of medical signals [9], research on river runoff in natural environments [10], research on rainwater evaporation [11], problems of polymer pollution [12], the atmospheric system [13], the influence of cosmic rays on organisms [14], and the motion of celestial bodies [15].

However, the Morlet wavelet does not satisfy the permissibility condition of wavelets, and wavelet reconstruction therefore cannot be realised. To obtain the inverse transformation of the Morlet wavelet, many scholars have improved the Morlet wavelet and obtained corresponding reconstruction transforms. Grossmann et al. [16] added proper correction terms to the Morlet wavelet to satisfy the permissibility condition and properly set the parameters for the tuning signal so that the rounding error of the computer had the same order of magnitude as the correction terms so that the correction terms could be omitted [17]. Ji and Yan [17], undertaking research at Northwestern Polytechnical University in China, improved the Morlet wavelet. They fixed some parameters in the Morlet wavelet and put forward the Morlet and Ji transform (MJT) as well. The MJT can reconstruct a signal without needing to meet the admissible condition. It is a continuous transformation of a signal and has a good time-frequency localisation property. Partially referencing the MJT and on the basis of the Morlet wavelet, we improved the Morlet wavelet in another way, further advancing a method for adaptively and discretely analysing signals.

In this paper, we first summarise the MJT and then propose another method to improve the Morlet wavelet cluster. Subsequently, the time and frequency windows of the Morlet wavelet, which were improved by Ji and Yan [17], are analysed, as is the Morlet wavelet improved by us. Third, the concept of extremum frequency is proposed in this paper, and the extremum properties of the Morlet wavelet cluster as improved by us are then analysed. Based on the improved Morlet wavelet cluster and extremum frequency, this paper then presents a new method that adaptively and discretely analyses a vibration signal and simultaneously yields the YPT, which can completely reconstruct the original signal. Additionally, we continue by putting forward a smooth operator that can smooth the potentially distorted signal reconstructed after analysis using the YPT and filtered employing threshold filtering theory. Finally, through the inverse transformation of the YPT, a filter is constructed to verify the correctness and practicability of the YPT by combining a signal with high noise. At the same time, the YPT is compared with the DWT. As a supplement to the theory in engineering, the shock signals about a gun automatic mechanism are also analysed by using the theory in this paper, which can provide help for further extracting the features of shock signals in pattern recognition and fault diagnosis.

2. MJT and Improved Morlet Wavelet Cluster

This section first introduces the content of the MJT, and we then improve the Morlet wavelet cluster in the MJT to obtain our improved Morlet wavelet cluster and analyse the properties of the time and frequency windows for the latter improved wavelet cluster.

2.1. MJT

The analytical expression of the Morlet wavelet is$$\begin{array}{}\text{(1)}& mt=e^{-t^2/2}{e}^{j{\omega }_{0}t},\hspace{1em}{\omega }_0\ge 5.\end{array}$$

The Fourier transform of (1) is$$\begin{array}{c}\text{(2)}& \widehat{m}\omega ={{\displaystyle \int }}_{-\infty }^{+\infty }{e}^{j{\omega }_{0}t}{e}^{-t^2/2}{e}^{j\omega t}dt=\sqrt{2\pi }{e}^{-{\omega -{\omega }_0}^2/2}.\end{array}$$

Given a wavelet $\psi t\in L^1\cap L^2$, the permissibility condition of the wavelet transform is$$\begin{array}{}\text{(3)}& C_{\psi }={{\displaystyle \int }}_{-\infty }^{+\infty }\frac{{\widehat{\psi }\omega }^2}{\omega }d\omega <+\infty .\end{array}$$

We know $\widehat{\psi }0=0$ by (3), whereas $\widehat{m}0\ne 0$; thus, the Morlet wavelet does not satisfy the permissibility condition of the wavelet transform. Following Ji and Yan [17], let $\psi t=1/\sqrt{2\pi }{e}^{-t^2/2+j2\pi t}$, and the modified Morlet wavelet cluster is then obtained by stretching and translating $\psi t$. Their improved Morlet wavelet cluster is$$\begin{array}{}\text{(4)}& {\psi }_{f,\tau }t=\frac{f}{\sqrt{2\pi }}{e}^{-f^2{t-\tau }^2/2+j2\pi ft-\tau },\hspace{1em}f\in R,\tau \in R.\end{array}$$

Assuming signal $ht\in L^2$, the expression for the MJT is [17]$$\begin{array}{c}\text{(5)}& Jf,\tau =ht,{\psi }_{f,\tau }t={{\displaystyle \int }}_{-\infty }^{+\infty }ht{\overline{\psi }}_{f,\tau }tdt\\ ={{\displaystyle \int }}_{-\infty }^{+\infty }ht\frac{f}{\sqrt{2\pi }}{e}^{-f^2{t-\tau }^2/2-j2\pi ft-\tau }dt\\ ={{\displaystyle \int }}_{-\infty }^{+\infty }ht\frac{f}{\sqrt{2\pi }}{e}^{-f^2{\tau -t}^2/2+j2\pi f\tau -t}dt\\ =h\tau \ast {\psi }_{f,0}\tau .\end{array}$$

According to (5), the Fourier transformation of $Jf,\tau$ with respect to $\tau$ can be obtained:$$\begin{array}{c}\text{(6)}& \widehat{J}f={{\displaystyle \int }}_{-\infty }^{+\infty }Jf,\tau e^{-j2\pi f\tau }d\tau ={{\displaystyle \int }}_{-\infty }^{+\infty }{{\displaystyle \int }}_{-\infty }^{+\infty }ht\frac{f}{\sqrt{2\pi }}{e}^{-f^2{\tau -t}^2/2+j2\pi f\tau -t}dt{e}^{-j2\pi f\tau }d\tau \\ ={{\displaystyle \int }}_{-\infty }^{+\infty }{{\displaystyle \int }}_{-\infty }^{+\infty }{e}^{-f^2{\tau -t}^2/2+j2\pi f\tau -t}{e}^{-j2\pi f\tau }d\tau ht\frac{f}{\sqrt{2\pi }}dt\\ ={{\displaystyle \int }}_{-\infty }^{+\infty }ht{e}^{-j2\pi ft}dt\\ =Hf.\end{array}$$

Therefore, according to (6), the inverse transformation of $ht$ can be obtained:$$\begin{array}{c}\text{(7)}& ht={{\displaystyle \int }}_{-\infty }^{+\infty }Hf{e}^{j2\pi ft}df={{\displaystyle \int }}_{-\infty }^{+\infty }{{\displaystyle \int }}_{-\infty }^{+\infty }Jf,\tau e^{-j2\pi f\tau }d\tau e^{j2\pi ft}df.\end{array}$$

Formulas (5) and (7) are, respectively, the positive transformation and inverse transformation formulas of the MJT. According to (5), the signal $ht$ can be transformed into the $f\text{–}\tau$ domain for analysis. The MJT is a continuous transformation that has good time-frequency localisation properties [17], and it is a redundant transformation that requires large amounts of calculations and storage space.

2.2. Improved Morlet Wavelet Cluster

In this section, we improve the Morlet wavelet cluster in (4) and analyse the properties of the time and frequency windows for the improved Morlet wavelet cluster. The improved Morlet wavelet cluster, denoted as ${\phi }_{f,\delta }t$, is expressed as (8). Formula (9) is the Fourier transform of (8):$$\begin{array}{c}\text{(8)}& {\phi }_{f,\delta }t=\frac{f+\delta }{\sqrt{2\pi }}{e}^{-{f+\delta }^2{t}^2/2+j2\pi ft},\text{(9)}& {\widehat{\phi }}_{f,\delta }\gamma ={{\displaystyle \int }}_{-\infty }^{+\infty }{\phi }_{f,\delta }t{e}^{-j2\pi \gamma t}dt={{\displaystyle \int }}_{-\infty }^{+\infty }\frac{f+\delta }{\sqrt{2\pi }}{e}^{-{f+\delta }^2{t}^2/2}{e}^{j2\pi f-\gamma t}dt\\ =\frac{\mathrm{sgn}f+\delta }{\sqrt{\pi }}{{\displaystyle \int }}_{-\infty }^{+\infty }{e}^{-t^2}{e}^{j2\sqrt{2}\pi /f+\delta f-\gamma t}dt\\ =\mathrm{sgn}f+\delta e^{-2{\pi }^2{f-\gamma }^2/{f+\delta }^2}.\end{array}$$

Cauchy’s theorem on the two-dimensional complex plane can be used to compute the definite integral in (9), and it can also be obtained directly from the Poisson integral formula. ${\phi }_{f,0}t={\psi }_{f,0}t$ can be obtained by using (4) and (8). We will now analyse the properties of the time and frequency windows for ${\phi }_{f,\delta }t$, the centre of which in the time domain is denoted as $t^{\ast }$, which can be obtained by the following equation:$$\begin{array}{c}\text{(10)}& t^{\ast }=\frac{{{\displaystyle \int }}_{-\infty }^{+\infty }t{{\phi }_{f,\delta }t}^{2}dt}{{{\displaystyle \int }}_{-\infty }^{+\infty }{{\phi }_{f,\delta }t}^{2}dt}=\frac{{{\displaystyle \int }}_{-\infty }^{+\infty }t{f+\delta /\sqrt{2\pi }{e}^{-{f+\delta }^2{t}^2/2+j2\pi ft}}^{2}dt}{{{\displaystyle \int }}_{-\infty }^{+\infty }{f+\delta /\sqrt{2\pi }{e}^{-{f+\delta }^2{t}^2/2+j2\pi ft}}^{2}dt}\hspace{1em}=0.\end{array}$$

In the process of calculating the definite integral in (10), we obtain ${{\displaystyle \int }}_{-\infty }^{+\infty }t{{\phi }_{f,\delta }t}^{2}dt=0$ and ${{\displaystyle \int }}_{-\infty }^{+\infty }{{\phi }_{f,\delta }t}^{2}dt=f+\delta /2\sqrt{\pi }\mathrm{sgn}f+\delta$ by utilising the properties of the Gamma function. Of course, as far as (10) is concerned, we do not have to calculate the latter in detail if we obtain the former. We can conclude that $t^{\ast }$ is independent of $f$ and $\tau$ via (10); therefore, the centre of the time domain for ${\psi }_{f,0}t$ is also 0. Furthermore, the radius of the time window, which is denoted as $\Delta t$, can be obtained by the following equation:$$\begin{array}{c}\text{(11)}& \Delta t={\frac{{{\displaystyle \int }}_{-\infty }^{+\infty }{t-t^{\ast }}^2{{\phi }_{f,\delta }t}^{2}dt}{{{\displaystyle \int }}_{-\infty }^{+\infty }{{\phi }_{f,\delta }t}^{2}dt}}^{1/2}={\frac{{{\displaystyle \int }}_{-\infty }^{+\infty }{t}^2{f+\delta }^2/2\pi e^{-{f+\delta }^2{t}^2}dt}{{{\displaystyle \int }}_{-\infty }^{+\infty }{{\phi }_{f,\delta }t}^{2}dt}}^{1/2}={\frac{\mathrm{sgn}f+\delta }{2\pi f+\delta }\frac{\Gamma 3/2}{{{\displaystyle \int }}_{-\infty }^{+\infty }{{\phi }_{f,\delta }t}^{2}dt}}^{1/2}=\frac{\mathrm{sgn}f+\delta }{\sqrt{2}f+\delta }.\end{array}$$

Obviously, the radius of the time window for ${\psi }_{f,0}t$ is $1/\sqrt{2}f$, and $f$ together with $\delta$ jointly affects the radius of the time window for ${\phi }_{f,\delta }t$. When the absolute value of $f$ is large, the radii of the time windows for ${\psi }_{f,0}t$ and ${\phi }_{f,\delta }t$ decrease, but the latter can be adjusted by using the parameter $\delta$. Similarly, the centre ${\gamma }^{\ast }$ and radius $\Delta \gamma$ of the frequency window for ${\phi }_{f,\delta }t$ in the frequency domain are obtained as$$\begin{array}{}\text{(12)}& {\gamma }^{\ast }=\frac{{{\displaystyle \int }}_{-\infty }^{+\infty }\gamma {{\widehat{\phi }}_{f,\delta }\gamma }^{2}d\gamma }{{{\displaystyle \int }}_{-\infty }^{+\infty }{{\widehat{\phi }}_{f,\delta }\gamma }^{2}d\gamma }=f,\text{(13)}& \Delta \gamma ={\frac{{{\displaystyle \int }}_{-\infty }^{+\infty }{\gamma -{\gamma }^{\ast }}^2{{\widehat{\phi }}_{f,\delta }\gamma }^{2}d\gamma }{{{\displaystyle \int }}_{-\infty }^{+\infty }{{\widehat{\phi }}_{f,\delta }\gamma }^{2}d\gamma }}^{1/2}=\mathrm{sgn}f+\delta \frac{f+\delta }{2\sqrt{2}\pi }.\end{array}$$

We can conclude $\Delta t\times \Delta \gamma =1/4\pi$ via (11) and (13), which as a constant cannot be affected by $\delta$, but it can be separately adjusted by $\delta$, and the width of the frequency window can be linearly adjusted by $f$ and $\delta$; what is more, $f$ determines the centre of the frequency window.

3. Extremum Frequency

In this section, the concept of extremum frequency, which reflects the average rate of change for a sequence and provides a method for adaptively analysing signals, is first proposed. Based on this concept, the extremum and extremum frequency of the improved Morlet wavelet cluster are analysed, which paves the way for the adaptive analysis of signals in the next section.

3.1. The Extremum Frequency of a Sequence

Assuming a sequence $xnn=1,2,\dots ,M$, if $xn-m_{\text{l}}-1<xn-m_{\text{l}}=\dots ,xn-1=xn=xn+1,\dots =xn+m_{\text{h}}>xn+m_{\text{h}}+1$ and $m_{\text{l}}\in N,m_{\text{h}}\in N$, we regard $xn-m_{\text{l}},\dots ,xn-1,xn,xn+1,\dots ,xn+m_{\text{h}}$ as one maximum aggregate of $xnn=1,2,\dots ,M$. Similarly, if $xn-m_{\text{l}}-1>xn-m_{\text{l}}=\dots ,xn-1=xn=xn+1,\dots =xn+m_{\text{h}}<xn+m_{\text{h}}+1$ and $m_{\text{l}}\in N,m_{\text{h}}\in N$, we regard $xn-m_{\text{l}},\dots ,xn-1,xn,xn+1,\dots ,xn+m_{\text{h}}$ as one minimum aggregate of $xnn=1,2,\dots ,M$. To go a step further, assuming that the sequence $xnn=1,2,\dots ,M$ contains $P_h$ amounts of maximum aggregates and $P_l$ amounts of minimum aggregates as well, the extremum frequency of $xnn=1,2,\dots ,M$, denoted as $f_{\text{e}}$, can be defined as$$\begin{array}{}\text{(14)}& f_{\text{e}}\triangleq \frac{P_h+P_l}{M}.\end{array}$$

3.2. The Extremum of the Improved Morlet Wavelet Cluster

In this section, we let $f+\delta \ne 0$ and analyse the real part of ${\phi }_{f,\delta }t$, which is equal to $f+\delta /\sqrt{2\pi }{e}^{-{f+\delta }^2{t}^2/2+j2\pi ft}$ and denoted as ${\phi }_{\text{R}}t$:$$\begin{array}{c}\text{(15)}& {\phi }_{\mathrm{R}}t=\frac{f+\delta }{\sqrt{2\pi }}{e}^{-{f+\delta }^2{t}^2/2}\cos 2\pi ft,{\phi }_{\text{R}}^{\prime }t=\frac{f+\delta }{-\sqrt{2\pi }}{e}^{-{f+\delta }^2{t}^2/2}{f+\delta }^{2}t\cos 2\pi ft+2\pi f\sin 2\pi ft.\end{array}$$

Without loss of generality, let $f+\delta >0$. Some properties regarding the maximum of ${\phi }_{f,\delta }t$ are now described as follows:

(1) There is only one point ${\xi }_n\in -1/4f+n1/2f,1/4f+n1/2f$ that makes ${{\phi }^{\prime }}_{\text{R}}{\xi }_n=0$ in the interval $-1/4f+n1/2f,1/4f+n1/2f,n\in Z$.

Proof.$$\begin{array}{}\text{(16)}& {{\phi }^{\prime }}_{\text{R}}t=\frac{f+\delta }{-\sqrt{2\pi }}{e}^{-{f+\delta }^2{t}^2/2}\cos 2\pi ft{f+\delta }^{2}t+2\pi f\tan 2\pi ft.\end{array}$$

Let $\rho t={f+\delta }^{2}t+2\pi f\tan 2\pi ft$, then ${{\phi }^{\prime }}_{\text{R}}{\xi }_n=0\iff \rho t=0$:$$\begin{array}{}\text{(17)}& {\rho }^{\prime }t={f+\delta }^2+{2\pi f\sec 2\pi ft}^2>0.\end{array}$$

When $f>0$,$$\begin{array}{c}\text{(18)}& \underset{t\to {-1/4f+n1/2f}^{+}}{\lim }\rho t=-\infty ,\underset{t\to {1/4f+n1/2f}^{-}}{\lim }\rho t=+\infty .\end{array}$$

$\rho t$ is continuous in the interval $-1/4f+n1/2f,1/4f+n1/2f$, and there is therefore only one point, ${\xi }_n\in -1/4f+n1/2f,1/4f+n1/2f$, that makes $\rho {\xi }_n=0$, that is, only one point, ${\xi }_n\in -1/4f+n1/2f,1/4f+n1/2f$, that makes ${{\phi }^{\prime }}_{\text{R}}{\xi }_n=0$.

When $f<0$, similarly, there exists only one point ${\xi }_n\in -1/4f+n1/2f,1/4f+n1/2f$ that makes ${{\phi }^{\prime }}_{\text{R}}{\xi }_n=0$.

(2) ${\xi }_n-1/4f+n1/2f<{\xi }_n<1/4f+n1/2f$ is the only extreme point in the interval $-1/4f+n1/2f,1/4f+n1/2f,n\in Z$.

Proof. When $t\in -1/4f+n1/2f,1/4f+n1/2f$ and $f>0$,$$\begin{array}{}\text{(19)}& {{\phi }^{\prime }}_{\text{R}}t=\frac{f+\delta }{-\sqrt{2\pi }}{e}^{-{f+\delta }^2{t}^2/2}\cos 2\pi ft\rho t.\end{array}$$

(i) Let $n$ be an even number. If $t\in -1/4f+n1/2f,{\xi }_n$, then ${{\phi }^{\prime }}_{\text{R}}t>0$; if $t\in {\xi }_n,1/4f+n1/2f$, then ${{\phi }^{\prime }}_{\text{R}}t<0$. That is to say, ${\xi }_n$ is the only maximum point in the interval $-1/4f+n1/2f,1/4f+n1/2f$.

(ii) Let $n$ be an odd number; similarly, ${\xi }_n$ is the only minimum point in the interval $-1/4f+n1/2f,1/4f+n1/2f$.

When $t\in -1/4f+n1/2f,1/4f+n1/2f$ and $f<0$, ${\xi }_n$ is similarly the only extreme point in the interval $-1/4f+n1/2f,1/4f+n1/2f$.

Because ${{\phi }^{\prime }}_{\text{R}}-1/4f+n1/2f\ne 0$, ${{\phi }^{\prime }}_{\text{R}}1/4f+n1/2f\ne 0$, and with the information given above, the property of (2) has been proved.

3.3. The Discrete Sampling of the Improved Morlet Wavelet Cluster

We can obtain a sequence, which is denoted as $unn=1,2,\dots ,N_0$, if sampling ${\phi }_{\text{R}}t$ from the start time $t_{\text{s}}$ to the end time $t_{\text{e}}$ and by sampling interval $t_0$. Figure 1 shows a schematic diagram for sampling ${\phi }_{\text{R}}t$.

[figure(s) omitted; refer to PDF]

According to the parameters for sampling ${\phi }_{\text{R}}t$, there exists the following formula:$$\begin{array}{}\text{(20)}& N_0=\frac{t_{\text{e}}-t_{\text{s}}}{t_0}+1.\end{array}$$

According to the parameters for sampling ${\phi }_{\text{R}}t$ and the analysis of the extremum for the improved Morlet wavelet cluster in Section 3.2, we can obtain the numbers of the extreme aggregates for $unn=1,2,\dots ,N_0$. Let the numbers be $P$, which can be expressed as$$\begin{array}{c}\text{(21)}& P=2f{t}_{\text{e}}-t_{\text{s}}\text{or\hspace{0.17em}\hspace{0.17em}}P=2f{t}_{\text{e}}-t_{\text{s}}\pm 1.\end{array}$$

Combining (14) and (20) with (21), the extremum frequency of $unn=1,2,\dots ,N_0$, denoted as $f_{\text{em}}$, can be obtained:$$\begin{array}{c}\text{(22)}& f_{\text{em}}=\frac{P}{N_0}=\frac{t_{0}2f{t}_{\text{e}}-t_{\text{s}}}{t_{\text{e}}-t_{\text{s}}+t_0}\text{or\hspace{0.17em}\hspace{0.17em}}{f}_{\text{em}}=\frac{P}{N_0}=\frac{t_{0}2f{t}_{\text{e}}-t_{\text{s}}\pm t_0}{t_{\text{e}}-t_{\text{s}}+t_0}.\end{array}$$

Furthermore,$$\begin{array}{}\text{(23)}& \underset{t_{\text{e}}-t_{\text{s}}\gg t_0}{\underset{f{t}_{\text{e}}-t_{\text{s}}\to +\infty }{\lim }}{f}_{\text{em}}=2f{t}_0.\end{array}$$

The effect of analysing a sequence whose extremum frequency is $f_{\text{e}}$ may be better when we utilise the improved Morlet wavelet, if we make $f_{\text{em}}=f_{\text{e}}$. Because the Morlet wavelet, on the condition of $f_{\text{em}}=f_{\text{e}}$, will have better adaptability to the sequence, the following formula is therefore established:$$\begin{array}{}\text{(24)}& f=\frac{f_{\text{e}}}{2t_0}.\end{array}$$

For a continuous signal $ht$ defined in the interval $t_{\text{s}},t_{\text{e}}$ and containing $P_{\text{c}}$ amounts of extreme points, we can determine that the extremum frequency of the sequence obtained by sampling $ht$ is $f_{\text{ec}}=P_{\text{c}}{t}_0/t_{\text{e}}-t_{\text{s}}+t_0$ if setting the sampling interval equal to $t_0$. When we utilise the improved Morlet wavelet to analyse $ht$, we set $f_{\text{em}}=f_{\text{ec}}$. Formula (25) is therefore established:$$\begin{array}{}\text{(25)}& \underset{t_{\text{e}}-t_{\text{s}}\gg t_0}{\underset{f{t}_{\text{e}}-t_{\text{s}}\to +\infty }{\lim }}{f}_{\text{em}}=\underset{t_{\text{e}}-t_{\text{s}}\gg t_0}{\underset{f{t}_{\text{e}}-t_{\text{s}}\to +\infty }{\lim }}{f}_{\text{ec}}.\end{array}$$

It is equal to$$\begin{array}{}\text{(26)}& f=\frac{P_{\text{c}}}{2t_{\text{e}}-t_{\text{s}}}.\end{array}$$

On the basis of (26), the improved Morlet wavelet, whose parameter $f$ is equal to $P_{\text{c}}/2t_{\text{e}}-t_{\text{s}}$, will have better adaptability when analysing the signal $ht$. According to (23)–(25), we obtain $f=f_{\text{e}}/2t_0=P_{\text{c}}/2t_{\text{e}}-t_{\text{s}}$ on the conditions $t_{\text{e}}-t_{\text{s}}\to +\infty$ and $t_{\text{e}}-t_{\text{s}}\gg t_0$.

4. Analysing a Vibration Signal

4.1. Two Adjacent Morlet Wavelets

In this section, we analyse two adjacent Morlet wavelets improved by us. The centres of the frequency windows separately are $f_{k-1}$ and $f_k$, and the radii are $\Delta {\gamma }_{k-1}$ and $\Delta {\gamma }_k$, respectively, where $k\ge 2$ and $k\in N$. We made their radii just connect with each other in the frequency domain, the situations of which are shown in the following two cases and seen in Figures 2 and 3.

[figure(s) omitted; refer to PDF]

Let $f+\delta \ge 0$. According to Figures 2 and 3, the following formula must be true if we make the radii of the two adjacent Morlet wavelets just connect with each other in the frequency domain:$$\begin{array}{}\text{(27)}& \Delta {\gamma }_k=\Delta {\gamma }_{k-1}-f_{k-1}-f_k.\end{array}$$

It is equal to$$\begin{array}{}\text{(28)}& {\delta }_k=f_{k-1}+{\delta }_{k-1}-2\sqrt{2}\pi f_{k-1}-f_k-f_k.\end{array}$$

4.2. YPT

Assume that a signal $h_{0}t\in L^{2}R$, defined in the interval $t_{\text{s}},t_{\text{e}}$ and containing $P_{\text{c}}$ amounts of extreme points, is sampled to form a discrete sequence $vnn=1,2,\dots ,N_v$ by setting the sampling interval equal to $t_0$, and the extremum frequency of the sequence is $f_{\text{e}}$. According to Section 3.3, we know we should make $f=f_{\text{e}}/2t_0=P_{\text{c}}/2t_{\text{e}}-t_{\text{s}}$, when analysing the signal $h_{0}t$ by utilising the improved Morlet wavelet. We now transform the signal $h_{0}t$ by using the following equations, which are collectively named the “YPT”:$$\begin{array}{c}\text{(29)}& Y_{pk}t;f_k,{\delta }_{pk}=h_{k-1}t\ast {\phi }_{f_k,{\delta }_{pk}}t,Y_{nk}t;-f_k,{\delta }_{nk}=h_{k-1}t\ast {\phi }_{-f_k,{\delta }_{nk}}t,\\ h_{k}t=h_{k-1}t-Y_{pk}t;f_k,{\delta }_{pk}-Y_{nk}t;-f_k,{\delta }_{nk},\end{array}$$where $k\in N^{+}$ in the formulas in (29). On their basis, to go a step further, we take the Fourier transforms of the formulas in (29):$$\begin{array}{c}\text{(30)}& {\widehat{Y}}_{pk}\gamma ;f_k,{\delta }_{pk}=\mathcal{F}{Y}_{pk}t;f_k,{\delta }_{pk},{\widehat{Y}}_{nk}\gamma ;-f_k,{\delta }_{nk}=\mathcal{F}{Y}_{nk}t;-f_k,{\delta }_{nk},\\ {\widehat{h}}_k\gamma =\mathcal{F}{h}_{k}t.\end{array}$$

The equations in (29) dispersedly and adaptively transform the signal $h_{0}t$ into the time domain, and on the contrary, the formulas in (30) dispersedly and adaptively transform the signal $h_{0}t$ into the frequency domain. Taken together, both constitute the positive transformation of the YPT. $f_k$ and ${\delta }_k$ are adaptive to $h_{k-1}t$. We can also obtain the inverse transformation of the YPT according to the formulas in (29) and (30), which can be described by the following equations:$$\begin{array}{}\text{(31)}& h_{0}t={\displaystyle \sum }_{i=1}^k{Y}_{pi}t;f_i,{\delta }_{pi}+Y_{ni}t;-f_i,{\delta }_{ni}+h_{k}t,\text{(32)}& h_{0}t={\mathrm{ℱ}}^{-1}{\displaystyle \sum }_{i=1}^k{\widehat{Y}}_{pi}\gamma ;f_i,{\delta }_{pi}+{\widehat{Y}}_{ni}\gamma ;-f_i,{\delta }_{ni}+{\widehat{h}}_k\gamma .\end{array}$$

For the convenience of our expression in this paper, $Y_{pk}t;f_k,{\delta }_{pk}$ and $Y_{nk}t;-f_k,{\delta }_{nk}$ are uniformly named $Yt;f,\delta$; similarly, $h_{k}t$ is designated as $ht$, ${\widehat{Y}}_{pk}\gamma ;f_k,{\delta }_{pk}$ and ${\widehat{Y}}_{nk}\gamma ;-f_k,{\delta }_{nk}$ are designated as $\widehat{Y}\gamma ;f,\delta$, and ${\widehat{h}}_k\gamma$ is designated as $\widehat{h}\gamma$. According to (31), we can directly reconstruct the original signal $h_{0}t$ by adding signals $Yt;f,\delta$ and $ht$, having been analysed and processed, in the time domain. In addition, we can continue to further transform signals $Yt;f,\delta$ and $ht$ into signals $\widehat{Y}\gamma ;f,\delta$ and $\widehat{h}\gamma$ in the frequency domain and then continue to process signals $\widehat{Y}\gamma ;f,\delta$ and $\widehat{h}\gamma$ and reconstruct the original signal $h_{0}t$ through (32). On a similar principle, according to (32), we can reconstruct the original signal $h_{0}t$ through signals $\widehat{Y}\gamma ;f,\delta$ and $\widehat{h}\gamma$, having been analysed and processed, in the frequency domain. In addition, we can continue to transform signals $\widehat{Y}\gamma ;f,\delta$ and $\widehat{h}\gamma$ into signals $Yt;f,\delta$ and $ht$ in the time domain and then continue to process signals $Yt;f,\delta$ and $ht$ and reconstruct the original signal $h_{0}t$ through (31). Therefore, (31) and (32) provide us 4 types of YPT inverse transforms, namely, 2 types of YPT inverse transforms from the frequency domain to the time domain and 2 types of YPT inverse transforms from the time domain to the time domain.

4.3. Smoothing Operator

A reconstructed vibration signal can be seriously distorted in some local areas after being analysed and processed. It is usually necessary to properly smooth the signal for the specific transformation, which is helpful to reflect the change law of the original signal. In this section, we propose a smooth operator $S$ that can smooth the sequence corresponding to the signal. Before introducing the operator $S$, three operators, $A$, $B$, and $L$, are introduced first. Assuming that there is a finite sequence $X=x_1,x_2,\dots ,x_n$, the operator $A$ can perform the following operation on the sequence $X$:$$\begin{array}{}\text{(33)}& AX,k\triangleq x_1,\dots ,x_{k-1},\frac{x_{k-1}+x_k}{2},x_{k+1},\dots ,x_n,\hspace{1em}k=2,3,\dots ,n.\end{array}$$

In particular,$$\begin{array}{c}\text{(34)}& AX,1\triangleq \frac{x_n+x_1}{2},x_2,\dots ,x_n,AX,0\triangleq X.\end{array}$$

The operator $B$ performs the following operation:$$\begin{array}{}\text{(35)}& BX\triangleq AAA\dots AX,1\dots ,n-1,n.\end{array}$$

The operator $L$ performs the following operation:$$\begin{array}{}\text{(36)}& LX,k\triangleq x_{n-k+1},x_{n-k+2},\dots ,x_n,x_1,x_2,\dots ,x_{n-k},\hspace{1em}k\in Z.\end{array}$$

Obviously, $LX,n=LX,-n=X$. On the basis of operators $A$, $B$, and $L$, we can define the operator $S$:$$\begin{array}{}\text{(37)}& SX,w\triangleq \underset{w}{\underbrace{LBB\dots BX,-w}},\hspace{1em}w\in N^{+}.\end{array}$$

In particular, $SX,0\triangleq X$.

Through a large number of experiments on sequences obtained by the YPT and the filtering theory that we utilised in this paper, we found that the operator $S$ can effectively smooth sequences. However, we cannot excessively smooth a signal.

4.4. The Process of Analysing a Vibration Signal

From the analyses in Sections 4.1 and 4.2, if we analyse the signal $h_{0}t$ through the YPT, the kernel is determining the extremum frequency of the signal and controlling the radius of the frequency window. When $k\ge 1$, $f_k$ corresponds to $h_{k-1}t$. If $k\ge 2$, ${\delta }_k$ must meet the condition ${\delta }_k=f_{k-1}+{\delta }_{k-1}-2\sqrt{2}\pi f_{k-1}-f_k-f_k$, which has been discussed through (28) in Section 4.1, and we can let ${\delta }_1$ be 0. The radius of the frequency window will be zero if $f_k+{\delta }_k=0$ in the process of analysing the signal. However, this matter will make the process of the YPT stagnate. We therefore let ${\delta }_k=0$ again if $f_k+{\delta }_k=0$ to continue transforming the signal until the transformation reaches the required number of the transformation. We can denote this required number as $K_0$. The flow chart of the YPT can be divided into two parts, the first of which concerns analysing the signal, and the second concerns reconstructing the signal; both parts are shown in Figures 4 and 5, respectively.

[figure(s) omitted; refer to PDF]

The parameter model determines how to reconstruct a signal. For this paper, we chose $\text{model}=4$ to reconstruct a signal.

5. Analysis and Reconstruction of a Signal

In this section, the correctness and practicability of the YPT are examined from two aspects: signal analysis and reconstruction. Based on the above research, we analysed a signal according to the flow chart in Figure 4 and reconstructed it according to the flow chart in Figure 5, using $\text{model}=4$. We first analysed a signal containing high amounts of white noise, and we then reconstructed the signal after the process of filtering it. Threshold filtering was adopted as filtering theory to remove the noise in the signal, while the threshold filtering established by us is neither hard threshold filtering nor conventional soft threshold filtering but is rather a theory established on the probability that an interval signal belongs to noise. In summary, first, the probability function that judges the probability that an interval signal belongs to noise is established on the basis of the characteristics that the expectation of the noise is zero and that the noise has a uniform distribution; the threshold is then determined by the probability function, and at last, the interval signals belonging to noise are filtered according to the threshold function. In this paper, this filtering theory is named “probability filtering” and is not discussed in detail.

5.1. Analysis of a Signal

Noise was added to the ideal signal $ht=\sin t+\sin \sqrt{2}t+\pi /3,t\in -10\pi ,10\pi$, by utilising the random function $\mathrm{rand}t$ in the Matlab numerical calculation software environment. The maximum noise amplitude was 3, and the noisy signal is therefore expressed as $ht=\sin t+\sin \sqrt{2}t+\pi /3+3\mathrm{rand}t$. The noisy signal was sampled into 2048 data points, which are plotted in Figure 6.

[figure(s) omitted; refer to PDF]

From Figure 6, it is difficult to observe the change regulation of the noisy signal when we set the maximum noise amplitude to 3. After calculation, the signal-to-noise ratio (SNR) of the noisy signal was 3.13 dB. This noise can therefore be regarded as high noise. The signal $ht$ was transformed 6 times using the YPT, and $K_0$ was therefore equal to 6. For the convenience of calculation, we equally divided the transformed signal $h{K}_0$ into two parts when we stored the transformed signals in memory. There were a total of 14 transformed signals resulting from the YPT. To conveniently observe the transformed signals in the time domain, the real parts of the signals in three-dimensional space are plotted in Figure 7.

[figure(s) omitted; refer to PDF]

Because the noisy signal is real, the two transformed signals whose serial numbers are opposite are conjugates and have the same real part. According to Figure 7, the signals with serial numbers 1 and −1 are noise, and on the contrary, the signals with serial numbers 2 and −2 are useful signals and have higher amplitudes when compared with those of the others. Similarly, the signals with serial numbers 3 and −3 are also useful signals, although they have lower amplitudes when compared with those of the former. To observe them clearly, the real parts of the signals on the two-dimensional plane are plotted in Figure 8. The other signals approximate to straight lines in three-dimensional space because of their lower amplitudes, but they are not straight lines. From the YPT results, we can conclude that the YPT can separate noise and useful signals, for which the essential aspect is that noise is different from the useful signal at the extremum frequency.

[figure(s) omitted; refer to PDF]

We can obtain the transformed signals in the frequency domain if continuing to transform the transformed signals in the time domain into the frequency domain using the YPT. The real parts of the transformed signals in three-dimensional space are plotted in Figure 9. From Figure 9, we can clearly see that the signals with serial numbers 1 and −1 contain large amounts of noise, which accords with the characteristics that the expectation of the noise is zero and is uniformly distributed in the frequency domain. For the signals with serial numbers 2 and −2, each of them has a positive impact amplitude and a lower negative impact amplitude in the frequency domain.

[figure(s) omitted; refer to PDF]

The real parts of the signals with serial numbers 3 and −3 are separately plotted on two-dimensional planes in Figures 10 and 11. The signals in Figures 10 and 11 and those in Figure 9 have similar properties. For the convenience of expression in this paper, the signals whose serial numbers are opposite to each other are, respectively, named “positive and negative signals.” Their real parts have symmetry if the signal $ht$ is real.

[figure(s) omitted; refer to PDF]

5.2. Reconstruction of a Signal

For this section, we adopted $\text{model}=4$ to reconstruct the signal $ht$ from Section 5.1. First, the signals in the frequency domain were filtered employing probability filtering to remove a part of the noise, and we then transformed the filtered signals in the frequency domain into the time domain. After that, we continued to again filter the signals in the time domain using probability filtering to remove a part of the noise in the different space. We finally reconstructed the original signal $ht$ using the filtered signals in the time domain. The process of reconstructing the original signal employed the YPT inverse transform, the filtering theory, and the smoothing operator. In Figure 12, the real parts of the new filtered signals in the frequency domain are shown in three-dimensional space.

[figure(s) omitted; refer to PDF]

Compared with those shown in Figure 9, the signals with serial numbers 1 and −1 in Figure 12 significantly reduced the noise and retained a part of the useful signals, as was the case for the signals with serial numbers 2 and −2. However, for the serial number 2 signal, the positive impact amplitude was reduced slightly because the positive impact amplitude was affected by the noise and the different filtering order when the threshold function readjusted the values of the interval signals. To compare with the signals in Figures 10 and 11 and to also observe the effects of the YPT and filtering theory, the real parts of the new filtered signals with serial numbers 3 and −3 are separately plotted on two-dimensional planes in Figures 13 and 14. From Figures 13 and 14, they have properties similar to those of the new filtered signals with serial numbers 2 and −2.

[figure(s) omitted; refer to PDF]

What follows is to transform the new filtered signals in the frequency domain into the time domain, and we can then continue to filter them in the time domain by probability filtering. The new filtered signals in three-dimensional space are plotted in the time domain in Figure 15. Compared with the signals in Figure 7, each of the signals in Figure 15 generally has less noise. On the contrary, the new filtered signal with serial number 2 lost a large amount of energy in the filtering process, which exposes the lack of this filtering theory but does not influence the validity and practicability of the YPT.

[figure(s) omitted; refer to PDF]

Finally, according to (31), we reconstructed the original signal $ht$ through the new filtered signals in Figure 15. We then properly smoothed the signal using the smoothing operator $SX,w$, and the parameter $w$ was 5. The new signal processed by the YPT and ideal signal without noise correspond to the B curve and the C curve in Figure 16, respectively. From Figure 16, we can determine that the amplitude variation of the new signal was not as high as the ideal signal. Combining Figures 6 and 16, we can conclude that the effect of reconstructing a signal is ideal when utilising the YPT, probability filtering, and the smoothing operator. Compared with the noisy signal in Figure 6, the noise in the new signal was very low, and it clearly reflected its own law of change.

[figure(s) omitted; refer to PDF]

To further verify the validity of the YPT, it is compared with the DWT. We also decomposed the noisy signal into 4 layers by the DWT, and the wavelets for the DWT are 3 order Daubechies wavelets. For the noisy signal decomposed by the DWT, after removing the first three layers of detail signals, we can then obtain the best reconstructed signal by the inverse discrete wavelet transform (IDWT), which corresponds to the A curve in Figure 16. After calculation, the SNR of the noisy signal, A curve, and B curve is 3.13 dB, 10.61 dB, and 16.39 dB, respectively, and as far as denoising and the noisy signal in this paper are concerned, the YPT is therefore better than the DWT.

6. The Application of the YPT in Engineering

In this section, two types of serious fault cracks on the latch sheets of a gun automatic mechanism were set by us, and the shock signals that are obtained by sampling the data, respectively, when the gun automatic mechanism works in normal and two fault patterns were analysed by utilising the method in this paper. The firing frequency of the gun automatic mechanism is $10\text{Hz}$, and the main parameters of the sensors that we used in this paper are shown in Table 1. To obtain more details in the time domain, we set the sampling frequency to the maximum sampling frequency ($204.8\text{kHz}$) of the signal collector, and the sampling time is $5\text{s}$. On the contrary, for a shock signal, the sampling frequency should be as high as possible to utilise the signal near the resonant frequency of the sensor, although the output of the sensor may be distorted at high frequency. The environment of the experiment is shown in Figure 17, and the main measuring point is shown in Figure 17(d).

Table 1

The parameters of the sensor.

Performance	Value
Sensitivity	1.0 mV/g (±15%)
Measurement range	±5000 g pk
Frequency range	0.4 to 7500 Hz (±10%)
Electrical filter cutoff frequency	≥7.5 kHz (−10%)
Resonant frequency	≥50 kHz
Broadband resolution	0.02 g·rms (1 to 10 kHz)

[figure(s) omitted; refer to PDF]

Due to the bad working environment or high overload of the equipment, a lot of noise is usually contained in the shock signals. The noise reduction for the signals therefore needs to be done to eliminate the effect of the noise if we do pattern recognition or fault diagnosis for the equipment. For the gun automatic mechanism, there are two main sources of noise, one of which is at ammunition feeding mechanism that is numbered (1) in Figure 17(a) and another is at the base of the equipment, numbered (2) in Figure 17(a). The former is produced by the irregular friction and collision between the bullet chain and the body of the equipment, and the latter is produced by self-excited vibration on the base that is with the complex structure or the random vibration of the base such as a gun automatic mechanism that works on a tank or fighter plane. These two kinds of noise have the property that the frequency of them is lower when compared with that of the efficient signals that are produced by the free vibrations of each part of the gun automatic mechanism.

The actual effective time period of original signals was intercepted to be analysed. The original signals and the signals filtered by the YPT are plotted in Figure 18. Figures 18(a)–18(c) are the original signals, and Figures 18(d)–18(f) are the signals that correspond to Figures 18(a)–18(c) and are filtered by the method, respectively. To further observe the effects of the method, we transformed them into the frequency domain, which is shown in Figure 19. From Figure 19, we can conclude that contrast to original signals, the amplitude of the part in low frequency is significantly reduced, relatively to the part in high frequency, and the shape of the part in high frequency almost does not change. The results of the analysis show that the noise in shock signals is filtered effectively. The YPT therefore shows good performance in the analysis of the shock signal in engineering, which can provide help for further extracting features for shock signals in pattern recognition and fault diagnosis.

[figure(s) omitted; refer to PDF]

7. Summary and Outlook

Because the Morlet wavelet does not satisfy the permissibility condition of the wavelet transform, it has no inverse wavelet transform. On the contrary, the Morlet wavelet has good time-frequency characteristics. In this paper, it was easy to adaptively adjust the centre and radius of the frequency window for the Morlet wavelet after its improvement. The essence of the wavelet transform is that it is a type of filtering. For this paper, the essence of transforming a signal into the time and frequency domains is also a filter. The YPT in this paper can separate noise from the signal, the essential aspect of which is that noise and useful signals have different extreme frequencies. On the principle of extreme frequency, the original signal is decomposed into different components using the YPT.

In this paper, the noisy signal and the shock signals in engineering could also be complex or real. The extremum frequency is determined by the real part of the signal when analysing a complex signal. There will be two ways for us to conduct research if further analysing a multidimensional signal. For the first, we can analyse each dimension of the multidimensional signal using the YPT alone. For the other, we can extend the dimension of the improved Morlet wavelet so that it has the ability to analyse multidimensional signals. For the latter method, how to determine the centre and radius of the Morlet wavelet in the frequency domain needs to be researched more deeply. Matching a multidimensional signal to the Morlet wavelet in phase is worthy of study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors gratefully acknowledge the support from the National Natural Science Foundation of China (Grant nos. 51675491 and 51175480).

References

[1] J. Morlet, G. Arens, E. Fourgeau, "Wave propagation and sampling theory-part II: sampling theory and complex waves," Geophysics, vol. 47 no. 2, pp. 222-236, DOI: 10.1190/1.1441329, 1982.

[2] M. J. Shensa, "Discrete inverses for nonorthogonal wavelet transforms," IEEE Transactions on Signal Processing, vol. 44 no. 6, pp. 798-807, DOI: 10.1109/78.492533, 1996.

[3] M. J. Shensa, "The discrete wavelet transform: wedding the à trous and Mallat algorithms," IEEE Transactions on Signal Processing, vol. 40 no. 10, pp. 2464-2482, DOI: 10.1109/78.157290, 1992.

[4] H. X. Chen, M. J. Zuo, X. D. Wang, M. R. Hoseini, "An adaptive Morlet wavelet filter for time-of-flight estimation in ultrasonic damage assessment," Measurement, vol. 43, pp. 570-585, DOI: 10.1016/j.measurement.2010.01.001, 2010.

[5] N. Saravanan, K. I. Ramachandran, "A case study on classification of features by fast single-shot multiclass PSVM using Morlet wavelet for fault diagnosis of spur bevel gear box," Expert Systems with Applications, vol. 36, pp. 10854-10862, DOI: 10.1016/j.eswa.2009.01.053, 2009.

[6] B. P. Tang, W. Y. Liu, T. Song, "Wind turbine fault diagnosis based on Morlet wavelet transformation and Wigner-Ville distribution," Renewable Energy, vol. 35, pp. 2862-2866, DOI: 10.1016/j.renene.2010.05.012, 2010.

[7] H. R. Karimi, W. Pawlus, K. G. Robbersmyr, "Signal reconstruction, modeling and simulation of a vehicle full-scale crash test based on Morlet wavelets," Neurocomputing, vol. 93, pp. 88-99, DOI: 10.1016/j.neucom.2012.04.010, 2012.

[8] M. Ayad, D. Chikouche, N. Boukazzoula, M. Rezki, "Search of a robust defect signature in gear systems across adaptive Morlet wavelet of vibration signals," IET Signal Process, vol. 8 no. 9, pp. 918-926, DOI: 10.1049/iet-spr.2013.0439, 2014.

[9] A. Suvichakorn, H. Ratiney, A. Bucur, S. Cavassila, J. P. Antoine, "Toward a quantitative analysis of in vivo proton magnetic resonance spectroscopic signals using the continuous Morlet wavelet transform," Measurement Science and Technology, vol. 20 no. 10,DOI: 10.1088/0957-0233/20/10/104029, 2009.

[10] C. P. Kuang, P. Su, J. Gu, "Multi-time scale analysis of runoff at the Yangtze estuary based on the Morlet wavelet transform method," Journal of Mountain Science, vol. 11 no. 6, pp. 1499-1506, DOI: 10.1007/s11629-014-3049-6, 2014.

[11] M. A. Goulart, L. Sanches, M. T. Vilani, O. B. Pinto Júnior, "Analysis of evapotranspiration by Morlet wavelet in area of Vochysia divergens Pohl in Pantanal," Revista Brasileira De Engenharia Agrícola E Ambiental, vol. 19 no. 2, pp. 93-98, DOI: 10.1590/1807-1929/agriambi.v19n2p93-98, 2015.

[12] A. Suvichakorn, H. Ratiney, A. Bucur, S. Cavassila, J. P. Antoine, "Quantification method using the Morlet wavelet for magnetic resonance spectroscopic signals with macromolecular contamination," Proceedings of the 30th Annual International IEEE EMBS Conference, pp. 2681-2984, .

[13] S. Nicolay, G. Mabille, X. Fettweis, M. Erpicum, "Brief communication “a statistical validation for the cycles found in air temperature data using a Morlet wavelet-based method," Nonlinear Processes in Geophysics, vol. 17, pp. 269-272, DOI: 10.5194/npg-17-269-2010, 2010.

[14] N. Zarrouk, R. Bennaceur, "Extrapolating cosmic ray variations and impacts on life: Morlet wavelet analysis," International Journal of Astrobiology, vol. 8 no. 3, pp. 169-174, DOI: 10.1017/s1473550409990085, 2009.

[15] L. T. Liu, H. Hsu, E. W. Grafarend, "Normal Morlet wavelet transform and its application to the earth’s polar motion," Journal of Geophysical Research, vol. 112 no. B8,DOI: 10.1029/2006jb004895, 2007.

[16] A. Grossmann, R. Kronland-Martinet, J. Morlet, Reading and Understanding Continuous Wavelet Transforms, 1990.

[17] Z. H. Ji, S. G. Yan, "Improved Morlet wavelet and its complete transform," Computer Engineering and Applications, vol. 52 no. 14, pp. 12-18, 2016.