On the Discrete Harmonic Wavelet Transform

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𝜓(𝜔)=12𝜋)(2^𝑗)for2𝜋2^𝑗≤𝜔<4𝜋2^𝑗,0elsewhere. (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.

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^𝑗𝑥−𝑘)𝑑𝑥,𝑎^{0𝑥002𝑑𝑐}_𝑗,𝑘=2^𝑗∞_−∞𝑓(𝑥)𝜓(2^𝑗𝑥−𝑘)𝑑𝑥, (2.3)and the corresponding pair of complex coefficients for the terms of scaling function,𝑎_𝜑,𝑘=^∞_−∞𝑓(𝑥)𝜑^∗(𝑥−𝑘)𝑑𝑥,𝑎^{0𝑥002𝑑𝑐}_𝜑,𝑘=^∞_−∞𝑓(𝑥)𝜑(𝑥−𝑘)𝑑𝑥. (2.4)If 𝑓(𝑥) is real, then ^0𝑥0𝑎^02𝑑𝑐_𝑗,𝑘 is the complex conjugate of 𝑎_𝑗,𝑘, that is, ^0𝑥0𝑎^02𝑑𝑐_𝑗,𝑘=𝑎^∗_𝑗,𝑘, but to allow the general case, when 𝑓(𝑥) is complex, we will consider ^0𝑥0𝑎^02𝑑𝑐_𝑗,𝑘 and 𝑎^∗_𝑗,𝑘 as two different amplitudes. Then the expansion formulas (2.1) and (2.2) become [2]𝑓(𝑥)=^∞∞_{𝑗=−∞𝑘=−∞}{𝑎_𝑗,𝑘𝜓(2^𝑗𝑥−𝑘)+𝑎^{0𝑥002𝑑𝑐}_𝑗,𝑘𝜓^∗(2^𝑗𝑥−𝑘)},𝑓(𝑥)=^∞_𝑘=−∞{𝑎_𝜑,𝑘𝜑(𝑥−𝑘)+𝑎^{0𝑥002𝑑𝑐}_𝜑,𝑘𝜑^∗+(𝑥−𝑘)}_{∞∞𝑗=0𝑘=−∞}{𝑎_𝑗,𝑘𝜓(2^𝑗𝑥−𝑘)+𝑎^{0𝑥002𝑑𝑐}_𝑗,𝑘𝜓^∗(2^𝑗𝑥−𝑘)}. (2.5)

Our primary purpose is to compute the coefficients 𝑎_𝜑,𝑘,𝑎^{0𝑥002𝑑𝑐}_𝜑,𝑘,𝑎_𝑗,𝑘 and ^0𝑥0𝑎^02𝑑𝑐_𝑗,𝑘 of this expansion.

An important condition for the function is that^∞_−∞|𝑓(𝑥)|²𝑑𝑥<∞. (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; 𝑓₀ 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.10)into the first formula of (2.3), and we obtain the following integral𝑎_𝑗,𝑘=12𝜋^4𝜋2𝑗_2𝜋2^𝑗𝑒^{𝑖𝜔𝑘/2𝑗}𝑑𝜔^∞_−∞𝑓(𝑥)𝑒^{−𝑖𝜔𝑥}𝑑𝑥, (2.11)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.12)

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

Analogous transformation towards the computation of ^0𝑥0𝑎^02𝑑𝑐_2𝑗+𝑘 will lead us to the following [2]:^0𝑥0𝑎^02𝑑𝑐_2𝑗+𝑘=^2𝑗−1_𝑠=0𝑓_𝑁−(2^𝑗_+𝑠)𝑒^{2𝜋𝑖𝑠𝑘/2𝑗},𝑘=0,…,2^𝑗−1. (2.14)

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

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

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 𝑓(𝑥)=2sin2𝜋𝑥 and 𝑓(𝑥)=2cos2𝜋𝑥. Following the algorithm, we discretise the interval [0;1] into 𝑁=2^𝑗 equally spaced nods, and obtain discrete set of values of functions𝑓_𝑟=2sin2𝜋𝑟𝑁,𝑓_𝑟=2cos2𝜋𝑟𝑁,𝑟=0,…,𝑁−1. (3.1)

The fast Fourier transform (2.8) of the obtained discrete sequence gives us the set Fourier coefficients 𝑓_𝑚. Recalling that 𝑎₀=𝑓₀,𝑎₁=𝑓₁, and 𝑎_𝑁/2=𝑓_𝑁/2, we can easily find these three coefficients. Another part of coefficients from 𝑎₂^𝑗 to 𝑎₂^𝑗+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 𝑎₂^𝑗_+𝑘 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 𝑓(𝑥)=2sin2𝜋𝑥 and 𝑓(𝑥)=2cos2𝜋𝑥 are presented in Figures 2 and 3.

(a)

(b)

Figure 2 

Arbitrary given function: (a) sin2𝜋𝑥, (b) cos2𝜋𝑥 (dashed line), and its reconstructed clone (solid line) from wavelet coefficients for 𝑁=8.

(a)

(b)

Figure 3 

Arbitrary given function: (a) sin2𝜋𝑥, (b) cos2𝜋𝑥 (dashed line), and its reconstructed clone (solid line) from wavelet coefficients for 𝑁=16.

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 𝑎₁ in the code instead of 𝑎₀, 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𝜖_𝑁=|𝑓(𝑥_𝑟)−𝑓_rec(𝑥_𝑟)|,𝑟=0,…,𝑁−1, (3.2)where 𝑓_rec(𝑥_𝑟) is the value of the reconstructed point. The dependence of absolute error of the reconstruction of the function from ln𝑁 is represented in Figure 5 and for two partial cases, when 𝑁=8 and 𝑁=16 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.

Figure 4 

Absolute error of the reconstruction of 𝑓(𝑥)=2sin2𝜋𝑥 for 𝑁=8 (solid line) and 𝑁=16 (dashed line).

Figure 5 

Absolute error of the reconstruction of 2sin2𝜋𝑥 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 ln𝑁 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 𝑥⁻¹), 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.