Abstract
We proposed the combination of signal denoising technology and Hankel transforms algorithm which were both based on Haar wavelet decomposition. Therefore, one can achieve above two purposes simultaneously while the wavelet decomposition is carried out just for once. The principle and its derivation of the Haar wavelet method for Hankel transforms were put forward. Numerical examples and engineering application showed that the precision of Haar wavelet method is about magnitude of and ; it can maintain good accuracy when using fewer wavelet coefficients. Moreover, it has better anti-noise performance and better computational stability than the filter method, so it can be applied to the Hankel transforms with noisy data.
1. Introduction
The Hankel transforms (HT) arise naturally in the discussion of the problems posed in cylindrical coordinates and hence [1], as a result of separation of variables, involve Bessel functions. The HT is frequently used as a tool for solving numerous scientific problems. For example, the stratified model and cylindrical coordinates are widely used in geophysical research, and the HT arises in forward and inverse calculation with zero or first order. The general HT pair with the kernel being is defined aswhere is the th-order Bessel function of the first kind. Analytical evaluations are rare and their numerical computations are difficult because of the oscillatory behavior of the Bessel function and the infinite length of the interval.
To overcome these difficulties, various different techniques are available in the literature. The first is the fast HT as proposed by Siegman in [2]. Here, via an exponential change of variables, the problem is transformed in the space of the logarithmic coordinates and the fast Fourier transform in that space. The disadvantage is requiring sampling over an exponential grid, thereby leading to important errors for functions with an oscillating tail. Moreover, it is sensitive to the smoothness of the functions [3]. The second is the back-projection and projection-slice methods [4], which carry out the HT as a double integral by means of one of the standard integral representations of the Bessel functions. But these methods generally require the efficient implementation of Chebyshev and Abel transforms. The computational complexity is unfortunately [5]. In [6], the authors used Filon quadrature philosophy to evaluate zero-order HT. They separated the integrand into the product of slowly varying component and a rapidly oscillating one. This methods works quite well for computing for , but the error is appreciable for . And the calculation of inverse HT is more difficult as is no longer a smooth function but a rapidly oscillating one. In 1998, Yu et al. [7] gave another method to compute zero-order quasi-discrete HT by approximating the input function by a Fourier-Bessel series over a finite integration interval. It leads to a symmetric transformation matrix. And later in 2004, Guizar-Sicairos and Gutierrez-Vega [8] obtained a powerful scheme to calculate the HT of order by extending the zero-order HT algorithm of Yu to higher orders. Their algorithm is based on the orthogonality properties of Bessel functions. Similarly, Gupta et al. [9] has proposed to transform the integral kernel function using orthogonal exponential expansion. This expansion reduces the integral to a simple algebraic sum because the analytical solution of HT of an exponential function is readily available. The essence of Pravin method is adopting exponential functions for fitting or interpolating bases [10].
More recently, Postnikov [11] and Zykov and Postnikov [12] proposed, for the first time, a novel and powerful method for computing zero- and first-order HT by using Haar bases and piecewise-linear bases. After that, Singh et al. adopted the linear Legendre multi-wavelets [3], rationalized Haar (RH) wavelets [13], and the hybrid method of block-pulse and Legendre polynomials [1] for small random perturbation in data function situation [14]. The above methods can be cast into a general class as expansion of integrand by wavelets.
However, from the engineering perspective, the wavelet theory is exactly originated from earth science which is initially proposed for observed signal analysis, and it is widely used in the signal denoising. Therefore, in this paper, we proposed to combine the signal denoising technology and the HT calculation algorithm which are both based on wavelet decomposition, so as to achieve the above two purposes simultaneously by choosing the decomposed coefficient and its threshold, while the wavelet decomposition is carried out just for once.
The criterion for selecting appropriate wavelet is that the chosen wavelets and their Bessel integral are easy to calculate. As for filtering, most wavelets can easily satisfy. Hence, the simplest and orthogonal Haar wavelet is adopted. Firstly, the principle of Haar wavelet decomposition is presented and the HT algorithm based on wavelet is given. Then, numerical examples and EM sounding application with comparison are carried out to illustrate the feasibility and performance of proposed algorithms.
2. Principle of Haar Wavelet Decomposition
Haar wavelet is the simplest, orthogonal, and compact support wavelet, and its scaling function and generating function are defined asThe corresponding binary scale translations are where and are integer referred separately as scale factor and translation factor. The scale space which is expanded by possess orthogonality just in scale , while the wavelet space which is expanded by possesses orthonormality in all scales. For arbitrary function , it can be decomposed in scale asHere, the is the projection of in scale space , which is corresponding to the smooth approximation of under scale and reflects the overall information. The is the projection of in wavelet space and reflects the detail information under scale . All the scale coefficient and wavelet coefficient can be directly determined by the inner product of and the corresponding basis functions as
Babolian and Shahsavaran [15] have proved that the approximation error in (4) is determined bywhere is the maximum absolute derivative value of and is the number of the terms in (4). So the error bounds of the wavelet approximation have an inverse relationship to the decomposition level, and the larger the means the better the convergence.
3. Hankel Transform Based on Haar Wavelet Decomposition
The wavelet method for HT is essentially a kind of integrand expansion. The selected wavelet function should satisfy the following conditions. Firstly, it should have good approximation ability for integrand, and meanwhile, an efficient coefficient calculation algorithm is available. In addition, its HT should have an analytical or simple expression. In this paper, the chosen function is Haar wavelet. It can acquire any precision as long as enough terms are adopted. And its decomposition coefficients can be efficiently calculated by the Mallat algorithm [16] which is similar to fast Fourier transform. Moreover, the HT of Haar function is corresponding to integral of the Bessel function and therefore has simple analytic expression. The Haar wavelet method will be illustrated in detail in the following.
(1) The Integrand Is Decomposed by Haar Wavelet. The first step is to determine the proper scale according to the precision or signal to noise ratio (SNR). And then, one can determine the truncation length according to the attenuation. Finally, the integrand can be expressed by scale function and wavelet function asHere, the scale and wavelet coefficients can be calculated efficiently by Mallat algorithm. And based on multi-resolution analysis, for the discrete integrand interval and truncation data number , one can derive that, in the scale , the domain length of scale function and wavelet function are both and the total number of the scale function is . But for the number of wavelets, it is for scale , 2m for scale , and so on. The total number of all the scales is 2N.
(2) HT in Wavelet Decomposition Form. Substitute the expansion into the transform formula and apply truncation processing. Then scale and wavelet functions are integral in domain, respectively. We can obtain the following:
(3) HT of Haar Function. The above integral is exactly to calculate the -order Bessel integral which can employ the following integral formula:Specifically, the expression of zero-order integral isThe first-order integral expression is
(4) Calculation of HT. It can be seen clearly that as the scale increases, the smooth approximation error of integrand increases; the detail part corresponding to the wavelet approximation in this scale gets vaguer. So the total approximation error is comparatively large. Hence, in order to improve the approximation accuracy, the wavelet components under scales , , , and so on can be gradually added. As the scale decreases, the amplitude of the wavelet component gradually reduces. However, for the wavelet filtering, the noise signals just correspond to the wavelet components with smaller amplitude. Hence, one can control approximation accuracy and filtering by selecting the scale of the adding wavelet. In fact, one can utilize easily the smooth approximation and its corresponding details, the balance between approximation accuracy and denoising performance can be controlled by adjusting . So after determining the above parameters, for the given variable , the corresponding HT can be calculated by (8).
4. Numerical Experiments
In this section, four typical HT [17] which include zero and first order were investigated. The integrands involve the monotonous and nonmonotonic cases. The HT experiments for noisy data and no noise situation were implemented. And the Haar method was compared with the classical filter method. The HT pairs were listed as
4.1. The Integrand Approximation Error Analysis
In the following experiments, the discrete interval of above integrands was set as ; the truncation length is . So the ratios of maximum to minimum for the integrands were all greater than , which indicated that the outside integrands value is far less than the chosen. Then, setting a single scale , the corresponding scale and wavelet coefficient numbers were both 64. The integrands were approximated by the smooth approximation and the wavelet approximation and the approximation errors were shown in Table 1.
If the scale is set, the scale and wavelet coefficient numbers were both 128, and the approximation errors were shown in Table 2.
The RMSE in Tables 1 and 2 is defined as
From Tables 1 and 2, we can find that the approximation error of Haar function is comparatively large for continuous functions because of its constant in the domain. The approximation error is associated with scale, and in general, the maximum error is reduced by half when the scale is reduced one level.
4.2. Hankel Transforms without Noise
Firstly, setting and the interval was 0.1. In the filter method, two kinds of node schemes were adopted, namely, the 61 and 120 node schemes for zero-order case and the 47 and 140 node schemes for first-order case [18]. The scale for wavelet method is . The number of truncation terms in Bessel function integral (as shown in (10)) is .
And secondly, we extend the region of to with interval 2.0. The maximum absolute errors and their comparisons were shown in Table 3.
We can conclude that the overall accuracy of wavelet method is about magnitude of to ; this is mainly caused by the limited approximation accuracy for integrand due to the constant value characteristics of the scale and wavelet functions. Numerical results also showed that when adding the wavelet detail information of scales and , even directly choosing the original discrete sequence, the accuracy improvement is still limited. It is identical to integrand approximation that the maximum error is reduced by half while the scale is reduced one level. Hence, we can increase the scale for better efficiency while a smaller error loss. Simultaneously, the truncated level can also affect the precision; numerical experiment showed that when increase from 11 to 31, the accuracy can be improved. For example, the maximum absolute error in Table 3 gets smaller as . But the cost is the time consuming.
However, for the filter method, we can find easily that increasing the filter node number may lead to better accuracy. The filter method is only adapted to specific intervals. For example, in the interval , the filter method may be unstable (for integrals 2 and 4), while the Haar wavelet method is stable.
From the time consuming aspect, the filter method has undoubtedly significant advantages. The time consumed is about 4~6 ms for integral 1. But for Haar wavelet method, the coefficient calculation time is about 10~12 ms and the integral calculation time is approximately 1.7~1.8 s (or 440 ms for ). However, the wavelet method has significant computational stability.
Then, consider a well-known HT [1] in optical diffraction theory where and areThe is the optical transfer function of an aberration-free optical system with a circular aperture and is the corresponding spread function. In numerical simulations, the discrete interval of the integrand is set as and the truncation region is with data points. And setting and the interval is 0.05. In the filter method, the high precision scheme of 120 nodes was adopted. While in wavelet method, the truncated level was set and scale . The absolute error curves and comparison were shown in Figure 1. The results showed that the maximum absolute error and RMSE in filter method were and . However, for wavelet method, the corresponding errors were and .
As can be seen from Figure 1, the accuracy of Haar wavelet method is about the magnitude of to , while the filter method is about magnitude of to . So the wavelet method’s accuracy is much higher than the filter’s in this experiment.
So we can conclude that the filter method is much suitable for the integrand with exponential form because the filter method utilizes an exponential replacement of variables, and its coefficients are calculated by exponential integrand optimally. However, the Haar wavelet method is based on integrand decomposition and can be adopted much more widely. Its defect lies in the truncation error, approximation error, and relative larger time consuming.
And finally, the HT of noisy data situation was investigated. The parameters were set the same as [1], so as to compare conveniently. The noisy data function is , where is a uniform random variable with values in . Three different noise terms as , , and were employed. The error curves of Haar method were shown in Figure 2. Hence, the errors were in the range of while the errors of hybrid method of Block-pulse and Legendre polynomials in [1] were in the range of .
4.3. Hankel Transforms of Noisy Data
Owning to the inevitable noise in practical application, the HT with noisy data was simulated and error variation was analyzed. The integrands were chosen as the former. Set and the interval is 0.05 and consider the Gaussian noise with 20 dB, 40 dB, and 60 dB situations. The curves of HT for integral 1 are shown in Figure 3 with 20 dB noise.
It can be obviously seen that the transformed function in the filter method contains significant noise because it directly substitutes the noisy data in integrand. However, for the Haar wavelet method, the small amplitude noise is decomposed in detail space to different scales. In the implementation of transforms, just the smooth approximation and detail part in scale were considered; the noisy signals smaller than scale were not considered. So the effect of noise has not been reflected consequently, the transformed function by the Haar wavelet method approximates the exact solution very well.
Hereafter, the noisy data experiments with 40 dB and 60 dB were carried out. The absolute error curves comparison is shown in Figure 4. As the SNR increases, the influence of noise gets smaller, and the accuracy of the numerical calculation improves. For the filter method, the random noise of integrand directly affects the transformed function, and the error distribution is also random. But for the Haar wavelet method, the selected smooth and detail part in scale for integrand approximation reflects the essence of filtering, so the transformed function is smooth, and its error curve is also relatively regular.
(a)
(b)
In addition, we introduced the root mean square error (RMSE) to illustrate the overall error distribution in the above four integrals. The comparison with noise was shown in Table 4.
As can be seen, the Haar wavelet method possesses good anti-noise performance because of its filtering effect on the integrand. It can maintain high accuracy even with larger noise (for 20 dB). However, the filter method is greater impact of noise. As the SNR increases 20 dB which means noise is reduced to about 10 percent, the RMSE is also reduced to 10 percent. So it is consistent with the theoretical analysis.
Then, different scales as and 4 were investigated. The corresponding scale and wavelet coefficient numbers were 16 and 256, respectively. We took the average of 100 times experiments for a statistical average. The root mean square errors were shown in Table 5. We can find that as the scale increasing, the coefficient numbers decreasing, and, however, the RMSE errors get larger.
5. Haar Wavelet Method in Electromagnetic Sounding Application
In DC sounding, consider the following three-hierarchical-layer model: = 400 Ω·m, = 15 m, Ω·m, = 20 m, = 1000 Ω·m. From recurrence relation, the apparent resistivity conversion function isHT of the apparent resistivity and resistivity conversion function is as follows:We considered using the apparent resistivity to calculate the conversion function in sounding application. Firstly, the apparent resistivity is given in the forward calculation and then the integrand obtained in (16). Finally, calculating the is corresponding to the first-order HT. In filter method, the 47-node scheme is adopted, while the scale is chosen in Haar wavelet method. The ranges of and are correspondingly set as and whose intervals are 1 and 0.002.
When the integrand does not contain noise, the curves of conversion function are shown in Figure 5. The numerical results show that the maximum absolute error, maximum relative error, and RMSE in filter method are 16.7422, 4.1886%, and 9.7080. However, for wavelet method, the corresponding errors are 81.2516, 14.8812%, and 12.6663. The loss of accuracy is mainly caused by the truncation and approximation error.
Then, the noisy data for HT were carried out. Numerical simulations showed that the filter method is unstable until the SNR is greater than 80 dB. However, for Haar wavelet method, it can maintain its accuracy even for large noise as SNR = 20 dB. The curves of 20 dB, 40 dB, and 60 dB Gaussian noise were shown in Figure 6. The maximum absolute error, maximum relative error, and RMSE in filter method are 16.6394, 90.8510%, and 14.6798 in case of SNR = 20 dB. Form Figure 7, we can find that the transform is still stable for 20 dB noise. So the wavelet method has good anti-noise performance and significant computational stability.
6. Conclusion
We have obtained a stable algorithm for HT, and the following conclusions can be drawn:(1)The combination of signal denoising technology and Hankel transforms algorithm was proposed in this paper. The common technology of them is the wavelet decomposition, and specifically, the Haar wavelet is chosen. Therefore, we can achieve above two purposes simultaneously while the wavelet decomposition is carried out just for once.(2)The principle and its derivation of Haar wavelet method for HT were given. Numerical examples and engineering application simulation showed that the Haar wavelet method has better computational stability than the filter method; its precision is about magnitude of and ; it can maintain good accuracy even using a few wavelet coefficients.(3)When considering the effects of noisy data, numerical experiments reveal that the Haar wavelet method only chooses the scale and wavelet coefficient in the same level, which is actually a filtering process, so it can get similar results to the noiseless case even with large noise. However, the filter method is easily affected by noise, leading to accuracy decrease or even unstable calculation. So the wavelet method has good anti-noise performance and significant computational stability. It can be applied to the HT with noisy data.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by the State Key Laboratory of Transmission Equipment and System Safety and Electrical New Technology Foundation Project of 2007DA10512713210.