Mathematical Problems in Engineering
Volume 2008 (2008), Article ID 687318, 7 pages
doi:10.1155/2008/687318
Research Article

On the Discrete Harmonic Wavelet Transform

Carlo Cattani1 and Aleksey Kudreyko2

1Department of Pharmaceutical Sciences (DiFarma), University of Salerno, Via Ponte Don Melillo, 84084 Fisciano (SA), Italy
2Department of Mathematics and Computer Science, University of Salerno, Via Ponte Don Melillo, 84084 Fisciano (SA), Italy

Received 29 May 2008; Accepted 26 July 2008

Academic Editor: Cristian Toma

Copyright © 2008 Carlo Cattani and Aleksey Kudreyko. 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.

Abstract

The discrete harmonic wavelet transform has been reviewed and applied towards given functions. The absolute error of reconstruction of the functions has been computed.

1. Introduction

The discrete harmonic wavelet transform was developed by Newland in 1993 [1, 2]. Similar to the ordinary discrete wavelet transform, the classical harmonic wavelet transform can also perform multiresolution analysis of a function. In addition, it has a fast algorithm based on fast Fourier transform for numerical implementation. A distinct advantage of harmonic wavelets is that they are disjoint in frequency domain (see Figure 1) and the Fourier transform of the successive levels decreases in propagation of their bandwidth (1.1). 𝜔 = ⎧ ⎪ ⎨ ⎪ ⎩ 1 𝜓 1 2 𝜋 2 𝑗 f o r 2 𝜋 2 𝑗 ≤ 𝜔 < 4 𝜋 2 𝑗 , 0 e l s e w h e r e . ( 1 . 1 ) Calculating its inverse Fourier transform, we obtain 𝜓 𝑗 𝑘 𝑥 = 𝑒 4 𝜋 𝑖 2 𝑗 𝑥 − 𝑘 − 𝑒 2 𝜋 𝑖 2 𝑗 𝑥 − 𝑘 2 2 𝜋 𝑖 𝑗 𝑥 − 𝑘 , ( 1 . 2 ) where 𝑗 = 0 , … , ∞ and 𝑘 = − ∞ , … , ∞ . This function represents a class of pulsed functions due to its compact support in the space domain.

687318.fig.001
Figure 1: Values of the Fourier transform of harmonic wavelets of different levels.

2. Discretisation of a Real Function

The goal of the wavelet transform is to decompose any arbitrary given function 𝑓 ( 𝑥 ) into an infinite summation of wavelets at different scales according to the expansion 𝑓 ( 𝑥 ) = ∞ ∞ 𝑗 = − ∞ 𝑘 = − ∞ 𝑎 𝑗 , 𝑘 𝜓 𝑗 𝑘 ( 𝑥 ) , ( 2 . 1 ) or in the alternative form [3] 𝑓 ( 𝑥 ) = ∞ 𝑘 = − ∞ 𝑎 𝜙 , 𝑘 𝜙 ( 𝑥 − 𝑘 ) + ∞ ∞ 𝑗 = 0 𝑘 = − ∞ 𝑎 𝑗 , 𝑘 𝜓 𝑗 𝑘 ( 𝑥 ) . ( 2 . 2 )

The first sum is a smooth approximation of 𝑓 ( 𝑥 ) , where the wavelets for 𝑗 ≤ 0 have been rolled together into scaling functions. The second sum is an addition of the details of 𝑓 ( 𝑥 ) at a specific level of resolution.

For complex wavelet coefficients, we have to define two amplitude coefficients 𝑎 𝑗 , 𝑘 = 2 𝑗 ∞ − ∞ 𝑓 𝑥 𝜓 ∗ 2 𝑗 𝑥 − 𝑘 𝑑 𝑥 , 𝑎 𝑗 , 𝑘 = 2 𝑗 ∞ − ∞ 𝑓 𝑥 𝜓 2 𝑗 𝑥 − 𝑘 𝑑 𝑥 , ( 2 . 3 ) and the corresponding pair of complex coefficients for the terms of scaling function, 𝑎 𝜑 , 𝑘 = ∞ − ∞ 𝑓 𝑥 𝜑 ∗ 𝑥 − 𝑘 𝑑 𝑥 , 𝑎 𝜑 , 𝑘 = ∞ − ∞ 𝑓 𝑥 𝜑 𝑥 − 𝑘 𝑑 𝑥 . ( 2 . 4 ) If 𝑓 ( 𝑥 ) is real, then 𝑎 𝑗 , 𝑘 is the complex conjugate of 𝑎 𝑗 , 𝑘 , that is, 𝑎 𝑗 , 𝑘 = 𝑎 ∗ 𝑗 , 𝑘 , but to allow the general case, when 𝑓 ( 𝑥 ) is complex, we will consider 𝑎 𝑗 , 𝑘 and 𝑎 ∗ 𝑗 , 𝑘 as two different amplitudes. Then the expansion formulas (2.1) and (2.2) become [2] 𝑓 𝑥 = ∞ ∞ 𝑗 = − ∞ 𝑘 = − ∞ 𝑎 𝑗 , 𝑘 𝜓 2 𝑗 𝑥 − 𝑘 + 𝑎 𝑗 , 𝑘 𝜓 ∗ 2 𝑗 , 𝑓 𝑥 = 𝑥 − 𝑘 ∞ 𝑘 = − ∞ 𝑎 𝜑 , 𝑘 𝜑 𝑥 − 𝑘 + 𝑎 𝜑 , 𝑘 𝜑 ∗ + 𝑥 − 𝑘 ∞ ∞ 𝑗 = 0 𝑘 = − ∞ 𝑎 𝑗 , 𝑘 𝜓 2 𝑗 𝑥 − 𝑘 + 𝑎 𝑗 , 𝑘 𝜓 ∗ 2 𝑗 . 𝑥 − 𝑘 ( 2 . 5 )

Our primary purpose is to compute the coefficients 𝑎 𝜑 , 𝑘 , 𝑎 𝜑 , 𝑘 , 𝑎 𝑗 , 𝑘 and 𝑎 𝑗 , 𝑘 of this expansion.

An important condition for the function is that ∞ − ∞ | | 𝑓 𝑥 | | 2 𝑑 𝑥 < ∞ . ( 2 . 6 )

Let us consider a real-valued function 𝑓 ( 𝑥 ) , represented by its discrete sequence 𝑓 𝑟 , 𝑟 = 0 , 1 , … , 𝑁 − 1 , ( 2 . 7 ) where 𝑁 = 2 𝑗 . Recalling the definition of the discrete Fourier transform, the corresponding Fourier coefficients are 𝑓 𝑚 = 1 𝑁 𝑁 − 1 𝑟 = 0 𝑓 𝑟 𝑒 − 2 𝜋 𝑖 𝑚 𝑟 / 𝑁 , 𝑚 = 0 , 1 , … , 𝑁 − 1 . ( 2 . 8 ) Note that 𝑓 𝑁 − 𝑚 = 1 𝑁 𝑁 − 1 𝑟 = 0 𝑓 𝑟 𝑒 − 2 𝜋 𝑖 ( 𝑁 − 𝑚 ) 𝑟 / 𝑁 = 1 𝑁 𝑁 − 1 𝑟 = 0 𝑓 𝑟 𝑒 − 2 𝜋 𝑖 𝑟 𝑒 2 𝜋 𝑖 𝑚 𝑟 / 𝑁 = 𝑓 ∗ 𝑚 , ( 2 . 9 ) where the asterisk stands for the complex conjugate; 𝑓 0 and 𝑓 𝑁 / 2 are always real numbers.

Furthermore, we will consider the coefficient 𝑎 𝑗 , 𝑘 , defined by the first formula in (2.3). Firstly, we will substitute 𝜓 ∗ 𝑗 , 𝑘 ( 𝑥 ) in terms of its Fourier transform (1.1) 𝜓 ∗ 𝑗 , 𝑘 1 ( 𝑥 ) = 2 𝑗 4 𝜋 2 𝑗 2 𝜋 2 𝑗 1 𝑒 2 𝜋 𝑖 𝜔 𝑘 / 2 𝑗 𝑒 − 𝑖 𝜔 𝑥 𝑑 𝜔 ( 2 . 1 0 ) into the first formula of (2.3), and we obtain the following integral 𝑎 𝑗 , 𝑘 = 1 2 𝜋 4 𝜋 2 𝑗 2 𝜋 2 𝑗 𝑒 𝑖 𝜔 𝑘 / 2 𝑗 𝑑 𝜔 ∞ − ∞ 𝑓 ( 𝑥 ) 𝑒 − 𝑖 𝜔 𝑥 𝑑 𝑥 , ( 2 . 1 1 ) where we have reversed the order of integration. The second integral over 𝑥 represents the Fourier transform of 𝑓 ( 𝑥 ) multiplied by 2 𝜋 , and (2.11) becomes 𝑎 𝑗 , 𝑘 = 4 𝜋 2 𝑗 2 𝜋 2 𝑗 𝑓 ( 𝜔 ) 𝑒 − 𝑖 𝜔 𝑘 / 2 𝑗 𝑑 𝜔 . ( 2 . 1 2 )

To derive a discrete algorithm of decomposition of the function, we must replace the operation of integration by summation, and (2.12) becomes 𝑎 2 𝑗 + 𝑘 = 2 𝑗 − 1 𝑠 = 0 𝑓 2 𝑗 + 𝑠 𝑒 2 𝜋 𝑖 𝑠 𝑘 / 2 𝑗 , 𝑘 = 0 , … , 2 𝑗 − 1 . ( 2 . 1 3 ) This identity represents the inverse discrete Fourier transform for the sequence of frequency coefficients 𝑓 2 𝑗 + 𝑠 .

Analogous transformation towards the computation of 𝑎 2 𝑗 + 𝑘 will lead us to the following [2]: 𝑎 2 𝑗 + 𝑘 = 2 𝑗 − 1 𝑠 = 0 𝑓 𝑁 − ( 2 𝑗 + 𝑠 ) 𝑒 2 𝜋 𝑖 𝑠 𝑘 / 2 𝑗 , 𝑘 = 0 , … , 2 𝑗 − 1 . ( 2 . 1 4 )

Computation of the amplitudes 𝑎 0 and 𝑎 𝑁 / 2 in the reviewed algorithm involves special approach, and 𝑎 0 = 𝑓 0 and 𝑎 𝑁 / 2 = 𝑓 𝑁 / 2 [2].

Also, it is easy to show from (2.13) that if 𝑗 = 0 , then 𝑘 = 0 and 𝑎 1 = 𝑓 1 . ( 2 . 1 5 )

Summarizing the stated above, the sequence of operations for computation of wavelet amplitude coefficients is as follows:

(i) represent the given function 𝑓 ( 𝑥 ) by a discrete sequence 𝑓 𝑟 , where 𝑟 = 0 , 1 , … , 𝑁 − 1 ;(ii) compute the set of frequency coefficients by fast Fourier transform 𝑓 𝑚 , where 𝑚 = 0 , 1 , … , 𝑁 − 1 ;(iii) the inverse fast Fourier transform of the octave blocks 𝑓 𝑚 generates the amplitudes of the harmonic wavelet expansion of the function 𝑓 𝑟 .

It is important to mention that this algorithm works for only the functions which satisfy the following conditions.

(i) The discrete transform covers the unit internal of 𝑥 .(ii) The analysed function is periodic in 𝑥 with period 1 .The algorithm was applied to the given functions which satisfy the mentioned conditions.

3. Implementation of Newland's Algorithm towards a Given Function

Let us review functions which satisfy the stated conditions. For example, it is 𝑓 ( 𝑥 ) = 2 s i n 2 𝜋 𝑥 and 𝑓 ( 𝑥 ) = 2 c o s 2 𝜋 𝑥 . Following the algorithm, we discretise the interval [ 0 ; 1 ] into 𝑁 = 2 𝑗 equally spaced nods, and obtain discrete set of values of functions 𝑓 𝑟 = 2 s i n 2 𝜋 𝑟 𝑁 , 𝑓 𝑟 = 2 c o s 2 𝜋 𝑟 𝑁 , 𝑟 = 0 , … , 𝑁 − 1 . ( 3 . 1 )

The fast Fourier transform (2.8) of the obtained discrete sequence gives us the set Fourier coefficients 𝑓 𝑚 . Recalling that 𝑎 0 = 𝑓 0 , 𝑎 1 = 𝑓 1 , and 𝑎 𝑁 / 2 = 𝑓 𝑁 / 2 , we can easily find these three coefficients. Another part of coefficients from 𝑎 2 𝑗 to 𝑎 2 𝑗 + 1 − 1 is obtained by computation of the inverse fast s Fourier transform (2.13) of coefficients from 𝑓 2 𝑗 to 𝑓 2 𝑗 + 1 − 1 .

To reconstruct the function from its wavelet coefficients, we followed the reverse algorithm of decomposition, that is: the fast Fourier transform of the wavelet coefficients 𝑎 2 𝑗 + 𝑘 represents the discrete Fourier transform of the reconstructed function 𝑓 𝑟 . Then, taking into account the shifting property (2.9), we can find 𝑓 as inverse fast Fourier transform of 𝑓 .

The results of decomposition and reconstruction of functions 𝑓 ( 𝑥 ) = 2 s i n 2 𝜋 𝑥 and 𝑓 ( 𝑥 ) = 2 c o s 2 𝜋 𝑥 are presented in Figures 2 and 3.

fig2
Figure 2: Arbitrary given function: (a) s i n 2 𝜋 𝑥 , (b) c o s 2 𝜋 𝑥 (dashed line), and its reconstructed clone (solid line) from wavelet coefficients for 𝑁 = 8 .
fig3
Figure 3: Arbitrary given function: (a) s i n 2 𝜋 𝑥 , (b) c o s 2 𝜋 𝑥 (dashed line), and its reconstructed clone (solid line) from wavelet coefficients for 𝑁 = 1 6 .

One can notice that the plots of the reconstructed functions are defined within the interval from 𝑟 = 1 to 𝑟 = 𝑁 . The difference between the algorithm and its corresponding computer code consists in that we put 𝑎 1 in the code instead of 𝑎 0 , and so forth . Therefore, the reconstruction of the function begins from point 1 / 𝑁 to 1, and not from 0 to 𝑁 − 1 .

To show the efficiency of the algorithm, it is worth to estimate the absolute error of the reconstructed function in the discrete nods. It is well known that the absolute error is given by 𝜖 𝑁 = | | 𝑓 ( 𝑥 𝑟 ) − 𝑓 r e c ( 𝑥 𝑟 ) | | , 𝑟 = 0 , … , 𝑁 − 1 , ( 3 . 2 ) where 𝑓 r e c ( 𝑥 𝑟 ) is the value of the reconstructed point. The dependence of absolute error of the reconstruction of the function from l n 𝑁 is represented in Figure 5 and for two partial cases, when 𝑁 = 8 and 𝑁 = 1 6 can be found in Figure 4. As we can see, small numbers of the level of decomposition 𝑗 give a very good approximation, when we reconstruct the function.

687318.fig.004
Figure 4: Absolute error of the reconstruction of 𝑓 ( 𝑥 ) = 2 s i n 2 𝜋 𝑥 for 𝑁 = 8 (solid line) and 𝑁 = 1 6 (dashed line).
687318.fig.005
Figure 5: Absolute error of the reconstruction of 2 s i n 2 𝜋 𝑥 after regression analysis.

4. Discussion of Results and Conclusion

Wavelets are considered as a new powerful tool for time-frequency analysis of nonlinear phenomena. In our paper, we discussed the harmonic wavelet transform and applied its algorithm towards decomposition and reconstruction of functions with a unit period. This algorithm might be useful for the wavelet solution of partial differential equations, when it is reduced to a system of ordinary differential equations [4, 5]. The algorithm of the decomposition consists of fast Fourier transform of the given discredited vector function, in which approximation error is proportional to l n 𝑁 and the corresponding approximation was obtained in our simulations (see Figure 5). It means that the increase of the length of 𝑁 leads us to a slow, but steady increase of the approximation error. The line of the dependence of the error from 𝑁 was obtained by implementing the method of least squares [6]. Note that the line of the plot takes discrete values due to the fact that 𝑁 takes only integer values of 2 𝑗 .

The only disadvantage of harmonic wavelets is that its decay rate is relatively low (proportional to 𝑥 − 1 ), therefore, its localisation is not precise. However, we have this disadvantage for the restricted Fourier transform of a harmonic wavelet of a specific level.

The application of harmonic wavelets towards particular problems is still new. The subject is developing very fast, however, there are still many questions remain unanswered. For example, what is the best choice of wavelet to use for a particular problem? How far does the harmonic wavelet transform computational simplicity compensate its slow decay rate in the 𝑥 -domain? How it can be used for the solution of integrodifferential equations, and many others. This work is in progress.

Acknowledgment

The work of Aleksey Kudreyko has been supported by Istituto Nazionale di Alta Matematica Francesco Severi (Rome-IT) under Scholarship U 2007/000394, 02/07/2007.

References

  1. D. E. Newland, “Harmonic wavelet analysis,” Proceedings of the Royal Society of London. Series A, vol. 443, no. 1917, pp. 203–225, 1993. View at MathSciNet
  2. D. E. Newland, An Introduction to Random Vibrations, Spectral & Wavelet Analysis, John Wiley & Sons, New York, NY, USA, 3rd edition, 1993. View at MathSciNet
  3. C. Cattani and J. Rushchitsky, Wavelet and Wave Analysis as Applied to Materials with Micro or Nanostructure, vol. 74 of Series on Advances in Mathematics for Applied Sciences, World Scientific, Hackensack, NJ, USA, 2007. View at MathSciNet
  4. C. Cattani and A. Kudreyko, “Mutiscale analysis of the Fisher equation,” in Proceedings of International Conference on Computational Science and Its Applications (ICCSA '08), vol. 5072 of Lecture Notes in Computer Science, pp. 1171–1180, Springer, Perugia, Italy, June-July 2008. View at Publisher · View at Google Scholar
  5. S. V. Muniandy and I. M. Moroz, “Galerkin modelling of the Burgers equation using harmonic wavelets,” Physics Letters A, vol. 235, no. 4, pp. 352–356, 1997. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  6. B. P. Demidovich, I. A. Maron, and E. Z. Shuvalova, Chislennie Metody Analiza, Lan', Moscow, Russia, 2008. View at MathSciNet