# A Stable Hankel Transforms Algorithm Based on Haar Wavelet Decomposition for Noisy Data

Authors: Huaiqing Zhang; Chuanbao Zhu; Xiangfeng Su; Xin Nie
Journal: Mathematical Problems in Engineering
 (2015)

---

## 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 of1.0E-4 and 1.0E-5; 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.


---

## Body

## 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 Jν is defined as(1)Fνr=∫0∞λfλJνλrdλ,fλ=∫0∞rFνrJνλrdr,where Jν 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 О(N2) [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 F0(r) for r⩾1, but the error is appreciable for r<1. And the calculation of inverse HT is more difficult as F0(r) 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 n 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 as(2)ϕt=1,t∈0,1,0,t∉0,1,ψt=1,t∈0,0.5,-1,t∈0.5,1,0,t∉0,1.The corresponding binary scale translations are (3)ϕj,kt=2-j/2ϕ2-jt-k,ψj,kt=2-j/2ψ2-jt-k,where j and k are integer referred separately as scale factor and translation factor. The scale space Vj which is expanded by Φj,k(t) possess orthogonality just in scale j, while the wavelet space Wj which is expanded by Ψj,k(t) possesses orthonormality in all scales. For arbitrary function f(t), it can be decomposed in scale J as(4)ft≈fmt=∑k∈ZcJ,kϕJ,kt+∑j=1J∑k∈Zdj,kψj,kt=PJft+D1ft+⋯+DJft.Here, the PJf(t) is the projection of f(t) in scale space VJ, which is corresponding to the smooth approximation of f(t) under scale J and reflects the overall information. The Djf(t) is the projection of f(t) in wavelet space Wj and reflects the detail information under scale j. All the scale coefficient cJ,k and wavelet coefficient dj,k can be directly determined by the inner product of f(t) and the corresponding basis functions as(5)cJ,k=ft,ϕJ,kt,dj,k=ft,ψj,kt.Babolian and Shahsavaran [15] have proved that the approximation error in (4) is determined by(6)ft-fmt2=M23·1m2,ft-fmt=O1m,where M is the maximum absolute derivative value of f(t) and m 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 m 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 J 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 λf(λ) can be expressed by scale function and wavelet function as(7)λfλ=∑k∈ZcJ,kϕJ,kλ+∑j=1J∑k∈Zdj,kψj,kλ.Here, 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 N, one can derive that, in the scale J, the domain length of scale function and wavelet function are both ΔJ=2JΔ and the total number of the scale function is m=N/2J. But for the number of wavelets, it is m for scale J, 2m for scale J-1, 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:(8)Fνr=∫0∞λfλJνλrdλ=∑k=0m-1cJ,k∫0∞ϕJ,kλ·Jνλrdλ+∑j=1J∑k=0N/2j-1dj,k∫0∞ψj,kλJνλrdλ=t=λr∑k=0m-1cJ,k1r∫0∞ϕJ,ktrJνtdt+∑j=1J∑k=0N/2j-1dj,k1r·∫0∞ψj,ktrJνtdt≈∑k=0m-1cJ,k1r∫rkΔJrk+1ΔJJνtdt+∑j=1J∑k=0N/2j-1dj,k·1r∫rkΔjrk+1/2ΔjJνtdt-∫rk+1/2Δjrk+1ΔjJνtdt.(3) HT of Haar Function. The above integral is exactly to calculate the ν-order Bessel integral which can employ the following integral formula:(9)∫0zJνtdt=2∑k=0∞Jν+2k+1z.Specifically, the expression of zero-order integral is(10)∫0zJ0tdt=2J1z+J3z+J5z+⋯=2limM→∞⁡∑m=0MJ2m+1z.The first-order integral expression is(11)∫0zJ1tdt=1-J0z.(4) Calculation of HT. It can be seen clearly that as the scale J 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 J-1, J-2, J-3, 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 J. So after determining the above parameters, for the given variable r, the corresponding HT Fν(r) 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(12)∫0∞λexp-λ2J0λrdλ=12exp-r24Integral1,∫0∞exp-4λJ0λrdλ=116+r2Integral2,∫0∞λ2exp-λ2J1λrdλ=r4exp-r24Integral3,∫0∞exp-3λJ1λrdλ=9+r2-3r9+r2Integral4.


### 4.1. The Integrand Approximation Error Analysis

In the following experiments, the discrete interval of above integrands was set ash=0.001; the truncation length is N=4096. So the ratios of maximum to minimum for the integrands were all greater than 1.0E+5, which indicated that the outside integrands value is far less than the chosen. Then, setting a single scale J=6, 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.Table 1
Approximation errors of integrands with Haar wavelet (J=6).
Errors
Integrand 1
Integrand 2
Integrand 3
Integrand 4
Max of integrand
4.2888
E
-
1
9.9601
E
-
1
3.6788
E
-
1
9.9700
E
-
1
Max of absolute errors
1.5491
E
-
2
5.9239
E
-
2
9.0990
E
-
3
4.4934
E
-
2
RMSE
3.1257
E
-
3
6.4337
E
-
3
2.3884
E
-
3
5.5765
E
-
3If the scaleJ=5 is set, the scale and wavelet coefficient numbers were both 128, and the approximation errors were shown in Table 2.Table 2
Approximation errors of integrands with Haar wavelet (J=5).
Errors
Integrand 1
Integrand 2
Integrand 3
Integrand 4
Max of absolute errors
7.4988
E
-
3
2.9272
E
-
2
4.4037
E
-
3
2.2089
E
-
2
RMSE
1.5606
E
-
3
3.2141
E
-
3
1.1926
E
-
3
2.7851
E
-
3The RMSE in Tables1 and 2 is defined as(13)ERMSE=1N∑i=1Nerri2.From Tables1 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, settingr∈[0.1,3.0] 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 J=6. The number of truncation terms in Bessel function integral (as shown in (10)) is M=11.And secondly, we extend the region ofr to [0.1,100.0] with interval 2.0. The maximum absolute errors and their comparisons were shown in Table 3.Table 3
The maximum absolute errors comparison between Haar wavelet (J=6) and filter method.
Error
Max of absolute errors in[0.1, 3.0]
Max of absolute errors in[0.1, 100]
Haar
61/47 point filter
120/140 point filter
Haar
61/47 point filter
Integral 1
2.4233
E
-
4
4.1481
E
-
6
1.6097
E
-
7
3.0357
E
-
2
2.5726
E
-
6
Integral 2
4.9964
E
-
4
1.4520
E
-
8
1.7021
E
-
11
8.0027
E
-
3
1.4520
E
-
8
Integral 3
9.4047
E
-
5
6.9370
E
-
6
3.8131
E
-
7
1.1050
E
-
4
3.6599
E
-
6
Integral 4
7.1477
E
-
5
6.4569
E
-
11
2.9526
E
-
13
7.3041
E
-
5
3.3231
E
-
11We can conclude that the overall accuracy of wavelet method is about magnitude of1.0E-4 to 1.0E-5; 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 J=5 and J=4, 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 J for better efficiency while a smaller error loss. Simultaneously, the truncated level M can also affect the precision; numerical experiment showed that when M increase from 11 to 31, the accuracy can be improved. For example, the maximum absolute error 3.0357E-2 in Table 3 gets smaller as 1.1379E-2. 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 intervalr<0, 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 forJ=8). However, the wavelet method has significant computational stability.Then, consider a well-known HT [1] in optical diffraction theory where f(λ) and Fν(r) are(14)fλ=2πarccosλ-λ1-λ2,0≤λ≤1,F0r=2J12r/2r2,0≤r<∞.The f(λ) is the optical transfer function of an aberration-free optical system with a circular aperture and F0(r) is the corresponding spread function. In numerical simulations, the discrete interval of the integrand g(λ)=λf(λ) is set as h=1.0/8192 and the truncation region is [h,1.0] with N=8192 data points. And setting r∈[0.01,20.0] 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 M=31 and scale J=6. 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 3.7499E-1 and 8.8241E-2. However, for wavelet method, the corresponding errors were 1.4209E-5 and 6.8513E-6.Figure 1
Absolute error curves comparison of Haar method and filter method with 120 nodes.As can be seen from Figure1, the accuracy of Haar wavelet method is about the magnitude of 1.0E-5 to 1.0E-7, while the filter method is about magnitude of 1.0E-1 to 1.0E-4. 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 gε(λi)=g(λi)+εθi, where θi is a uniform random variable with values in [-1,1]. Three different noise terms ε as ε0=0.0, ε1=0.002, and ε2=0.005 were employed. The error curves of Haar method were shown in Figure 2. Hence, the errors were in the range of [-1.0E-5,2.0E-5] while the errors of hybrid method of Block-pulse and Legendre polynomials in [1] were in the range of [-1.5E-4,5.0E-5].Figure 2
Error curves comparison of Haar method with different noise terms.


### 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. Setr∈[0.05,3.0] 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.Figure 3
Comparison of exact, Haar method and filter method in example 1 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 scaleJ were considered; the noisy signals smaller than scale J 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 Figure4. 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 J for integrand approximation reflects the essence of filtering, so the transformed function is smooth, and its error curve is also relatively regular.Figure 4
Absolute error comparison between Haar wavelet and filter method with noise of 40 dB and 60 dB.
(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 Table4.Table 4
The root mean square errors comparison between Haar wavelet (J=6) and filter method.
RMSE
20 dB
40 dB
60 dB
Haar
61/47 nods
Haar
61/47 nodes
Haar
61/47 nodes
Integral 1
4.6623
E
-
4
1.2357
E
-
2
1.8164
E
-
4
1.1841
E
-
3
1.6494
E
-
4
1.1110
E
-
4
Integral 2
8.7687
E
-
4
6.3719
E
-
3
4.6967
E
-
4
6.3924
E
-
4
4.8124
E
-
4
6.0250
E
-
5
Integral 3
2.1133
E
-
4
7.8903
E
-
3
5.8432
E
-
4
7.0230
E
-
4
5.5511
E
-
4
6.9836
E
-
4
Integral 4
5.4067
E
-
4
2.8169
E
-
3
6.7550
E
-
5
2.2125
E
-
4
5.4316
E
-
5
2.2098
E
-
5As 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 asJ=8 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.Table 5
The root mean square errors comparison between Haar wavelet scales.
RMSE
20 dB
40 dB
60 dB
J
=
8
J
=
4
J
=
8
J
=
4
J
=
8
J
=
4
Integral 1
3.5801
E
-
3
9.3433
E
-
4
4.9504
E
-
4
2.4083
E
-
4
2.9313
E
-
4
1.7519
E
-
4
Integral 2
3.3828
E
-
3
7.3215
E
-
4
6.4996
E
-
4
4.0977
E
-
4
9.3518
E
-
4
4.6637
E
-
4
Integral 3
2.0787
E
-
3
8.2472
E
-
4
5.3414
E
-
4
1.0870
E
-
4
4.2525
E
-
4
6.6308
E
-
5
Integral 4
3.0913
E
-
3
9.2945
E
-
4
7.5151
E
-
4
1.5243
E
-
4
7.4011
E
-
4
1.0482
E
-
4


## 4.1. The Integrand Approximation Error Analysis

In the following experiments, the discrete interval of above integrands was set ash=0.001; the truncation length is N=4096. So the ratios of maximum to minimum for the integrands were all greater than 1.0E+5, which indicated that the outside integrands value is far less than the chosen. Then, setting a single scale J=6, 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.Table 1
Approximation errors of integrands with Haar wavelet (J=6).
Errors
Integrand 1
Integrand 2
Integrand 3
Integrand 4
Max of integrand
4.2888
E
-
1
9.9601
E
-
1
3.6788
E
-
1
9.9700
E
-
1
Max of absolute errors
1.5491
E
-
2
5.9239
E
-
2
9.0990
E
-
3
4.4934
E
-
2
RMSE
3.1257
E
-
3
6.4337
E
-
3
2.3884
E
-
3
5.5765
E
-
3If the scaleJ=5 is set, the scale and wavelet coefficient numbers were both 128, and the approximation errors were shown in Table 2.Table 2
Approximation errors of integrands with Haar wavelet (J=5).
Errors
Integrand 1
Integrand 2
Integrand 3
Integrand 4
Max of absolute errors
7.4988
E
-
3
2.9272
E
-
2
4.4037
E
-
3
2.2089
E
-
2
RMSE
1.5606
E
-
3
3.2141
E
-
3
1.1926
E
-
3
2.7851
E
-
3The RMSE in Tables1 and 2 is defined as(13)ERMSE=1N∑i=1Nerri2.From Tables1 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, settingr∈[0.1,3.0] 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 J=6. The number of truncation terms in Bessel function integral (as shown in (10)) is M=11.And secondly, we extend the region ofr to [0.1,100.0] with interval 2.0. The maximum absolute errors and their comparisons were shown in Table 3.Table 3
The maximum absolute errors comparison between Haar wavelet (J=6) and filter method.
Error
Max of absolute errors in[0.1, 3.0]
Max of absolute errors in[0.1, 100]
Haar
61/47 point filter
120/140 point filter
Haar
61/47 point filter
Integral 1
2.4233
E
-
4
4.1481
E
-
6
1.6097
E
-
7
3.0357
E
-
2
2.5726
E
-
6
Integral 2
4.9964
E
-
4
1.4520
E
-
8
1.7021
E
-
11
8.0027
E
-
3
1.4520
E
-
8
Integral 3
9.4047
E
-
5
6.9370
E
-
6
3.8131
E
-
7
1.1050
E
-
4
3.6599
E
-
6
Integral 4
7.1477
E
-
5
6.4569
E
-
11
2.9526
E
-
13
7.3041
E
-
5
3.3231
E
-
11We can conclude that the overall accuracy of wavelet method is about magnitude of1.0E-4 to 1.0E-5; 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 J=5 and J=4, 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 J for better efficiency while a smaller error loss. Simultaneously, the truncated level M can also affect the precision; numerical experiment showed that when M increase from 11 to 31, the accuracy can be improved. For example, the maximum absolute error 3.0357E-2 in Table 3 gets smaller as 1.1379E-2. 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 intervalr<0, 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 forJ=8). However, the wavelet method has significant computational stability.Then, consider a well-known HT [1] in optical diffraction theory where f(λ) and Fν(r) are(14)fλ=2πarccosλ-λ1-λ2,0≤λ≤1,F0r=2J12r/2r2,0≤r<∞.The f(λ) is the optical transfer function of an aberration-free optical system with a circular aperture and F0(r) is the corresponding spread function. In numerical simulations, the discrete interval of the integrand g(λ)=λf(λ) is set as h=1.0/8192 and the truncation region is [h,1.0] with N=8192 data points. And setting r∈[0.01,20.0] 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 M=31 and scale J=6. 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 3.7499E-1 and 8.8241E-2. However, for wavelet method, the corresponding errors were 1.4209E-5 and 6.8513E-6.Figure 1
Absolute error curves comparison of Haar method and filter method with 120 nodes.As can be seen from Figure1, the accuracy of Haar wavelet method is about the magnitude of 1.0E-5 to 1.0E-7, while the filter method is about magnitude of 1.0E-1 to 1.0E-4. 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 gε(λi)=g(λi)+εθi, where θi is a uniform random variable with values in [-1,1]. Three different noise terms ε as ε0=0.0, ε1=0.002, and ε2=0.005 were employed. The error curves of Haar method were shown in Figure 2. Hence, the errors were in the range of [-1.0E-5,2.0E-5] while the errors of hybrid method of Block-pulse and Legendre polynomials in [1] were in the range of [-1.5E-4,5.0E-5].Figure 2
Error curves comparison of Haar method with different noise terms.


## 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. Setr∈[0.05,3.0] 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.Figure 3
Comparison of exact, Haar method and filter method in example 1 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 scaleJ were considered; the noisy signals smaller than scale J 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 Figure4. 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 J for integrand approximation reflects the essence of filtering, so the transformed function is smooth, and its error curve is also relatively regular.Figure 4
Absolute error comparison between Haar wavelet and filter method with noise of 40 dB and 60 dB.
(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 Table4.Table 4
The root mean square errors comparison between Haar wavelet (J=6) and filter method.
RMSE
20 dB
40 dB
60 dB
Haar
61/47 nods
Haar
61/47 nodes
Haar
61/47 nodes
Integral 1
4.6623
E
-
4
1.2357
E
-
2
1.8164
E
-
4
1.1841
E
-
3
1.6494
E
-
4
1.1110
E
-
4
Integral 2
8.7687
E
-
4
6.3719
E
-
3
4.6967
E
-
4
6.3924
E
-
4
4.8124
E
-
4
6.0250
E
-
5
Integral 3
2.1133
E
-
4
7.8903
E
-
3
5.8432
E
-
4
7.0230
E
-
4
5.5511
E
-
4
6.9836
E
-
4
Integral 4
5.4067
E
-
4
2.8169
E
-
3
6.7550
E
-
5
2.2125
E
-
4
5.4316
E
-
5
2.2098
E
-
5As 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 asJ=8 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.Table 5
The root mean square errors comparison between Haar wavelet scales.
RMSE
20 dB
40 dB
60 dB
J
=
8
J
=
4
J
=
8
J
=
4
J
=
8
J
=
4
Integral 1
3.5801
E
-
3
9.3433
E
-
4
4.9504
E
-
4
2.4083
E
-
4
2.9313
E
-
4
1.7519
E
-
4
Integral 2
3.3828
E
-
3
7.3215
E
-
4
6.4996
E
-
4
4.0977
E
-
4
9.3518
E
-
4
4.6637
E
-
4
Integral 3
2.0787
E
-
3
8.2472
E
-
4
5.3414
E
-
4
1.0870
E
-
4
4.2525
E
-
4
6.6308
E
-
5
Integral 4
3.0913
E
-
3
9.2945
E
-
4
7.5151
E
-
4
1.5243
E
-
4
7.4011
E
-
4
1.0482
E
-
4


## 5. Haar Wavelet Method in Electromagnetic Sounding Application

In DC sounding, consider the following three-hierarchical-layer model:ρ1 = 400 Ω·m, h1 = 15 m, ρ2=50 Ω·m, h2 = 20 m, ρ3 = 1000 Ω·m. From recurrence relation, the apparent resistivity conversion function T1(λ) is(15)T1λ=ρ11+K12e-2λh1+K23e-2λh1+h2+K12K23e-2λh21-K12e-2λh1-K23e-2λh1+h2+K12K23e-2λh2,Kij=ρj-ρiρj+ρi.HT of the apparent resistivity and resistivity conversion function is as follows:(16)ρsrr2=∫0∞T1λJ1λrλdλ,T1λ=∫0∞ρsrrJ1λrdr.We considered using the apparent resistivity ρs(r) to calculate the conversion function T1(λ) in sounding application. Firstly, the apparent resistivity ρs(r) is given in the forward calculation and then the integrand obtained ρs(r)/r in (16). Finally, calculating the T1(λ) is corresponding to the first-order HT. In filter method, the 47-node scheme is adopted, while the scale J=3 is chosen in Haar wavelet method. The ranges of r and λ are correspondingly set as [1,2048] and [0.002,0.256] whose intervals are 1 and 0.002.When the integrandρs(r)/r does not contain noise, the curves of conversion function T1(λ) 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.Figure 5
Comparison ofT1(λ) curves.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. TheT1(λ) 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.Figure 6
T
1
(
λ
) calculated by noisy data with wavelet.Figure 7
T
1
(
λ
) comparison with noisy data.


## 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 of1.0E-4 and 1.0E-5; 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.



---

*Source: 191086-2015-06-21.xml*