# On the Discrete Harmonic Wavelet Transform

Authors: Carlo Cattani; Aleksey Kudreyko
Journal: Mathematical Problems in Engineering
 (2008)

---

## 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.


---

## Body

## 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)ψ^(ω)={(12π)(12j)for2π2j≤ω<4π2j,0elsewhere.Calculating its inverse Fourier
transform, we obtain(1.2)ψkj(x)=e4πi(2jx−k)−e2πi(2jx−k)2πi(2jx−k),where j=0,…,∞ and k=−∞,…,∞.
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 functionf(x) into an infinite summation of wavelets at
different scales according to the expansion(2.1)f(x)=∑j=−∞∞∑k=−∞∞aj,kψkj(x),or in  the alternative form
[3](2.2)f(x)=∑k=−∞∞aϕ,kϕ(x−k)+∑j=0∞∑k=−∞∞aj,kψkj(x).The first sum is a smooth approximation off(x),
where the wavelets for j≤0 have been rolled together into scaling
functions. The second sum is an addition of the details of f(x) at a specific level of resolution.For complex wavelet coefficients, we have to define
two amplitude coefficients(2.3)aj,k=2j∫−∞∞f(x)ψ∗(2jx−k)dx,a˜j,k=2j∫−∞∞f(x)ψ(2jx−k)dx,and the corresponding pair of
complex coefficients for the terms of scaling function,(2.4)aφ,k=∫−∞∞f(x)φ∗(x−k)dx,a˜φ,k=∫−∞∞f(x)φ(x−k)dx.If f(x) is real, then a˜j,k is the complex conjugate of aj,k,
that is, a˜j,k=aj,k∗,
but to allow the general case, when f(x) is complex, we will consider a˜j,k and aj,k∗ as two different amplitudes. Then the
expansion formulas (2.1) and
(2.2) become [2](2.5)f(x)=∑j=−∞∞∑k=−∞∞{aj,kψ(2jx−k)+a˜j,kψ∗(2jx−k)},f(x)=∑k=−∞∞{aφ,kφ(x−k)+a˜φ,kφ∗(x−k)}+∑j=0∞∑k=−∞∞{aj,kψ(2jx−k)+a˜j,kψ∗(2jx−k)}.Our primary purpose is to compute the coefficientsaφ,k,a˜φ,k,aj,k and a˜j,k of this expansion.An important condition for the function is
that(2.6)∫−∞∞|f(x)|2dx<∞.Let us consider a real-valued functionf(x),
represented by its discrete sequence(2.7)fr,r=0,1,…,N−1,where N=2j.
Recalling the definition of the discrete Fourier transform, the corresponding
Fourier coefficients are(2.8)f^m=1N∑r=0N−1fre−2πimr/N,m=0,1,…,N−1.Note that(2.9)f^N−m=1N∑r=0N−1fre−2πi(N−m)r/N=1N∑r=0N−1fre−2πire2πimr/N=f^m∗,where the asterisk stands for
the complex conjugate; f^0 and f^N/2 are always real numbers.Furthermore, we will consider the coefficientaj,k,
defined by the first formula in (2.3). Firstly, we will substitute ψj,k∗(x) in terms of its Fourier transform
(1.1)(2.10)ψj,k∗(x)=12j∫2π2j4π2j12πeiωk/2je−iωxdωinto the first formula of (2.3),
and we obtain the following integral(2.11)aj,k=12π∫2π2j4π2jeiωk/2jdω∫−∞∞f(x)e−iωxdx,where we have reversed the order
of integration. The second integral over x represents the Fourier transform of f(x) multiplied by 2π,
and (2.11) becomes(2.12)aj,k=∫2π2j4π2jf^(ω)e−iωk/2jdω.To derive a discrete algorithm of decomposition of the
function, we must replace the operation of integration by summation, and (2.12)
becomes(2.13)a2j+k=∑s=02j−1f^2j+se2πisk/2j,k=0,…,2j−1.This identity represents the
inverse discrete Fourier transform for the sequence of frequency coefficients f^2j+s.Analogous transformation towards the computation ofa˜2j+k will lead us to the following [2]:(2.14)a˜2j+k=∑s=02j−1f^N−(2j+s)e2πisk/2j,k=0,…,2j−1.Computation of the amplitudesa0 and aN/2 in the reviewed algorithm involves special
approach, and a0=f^0 and aN/2=f^N/2 [2].Also, it is easy to show from (2.13) that if j=0,
then k=0 and(2.15)a1=f^1.Summarizing the stated above, the sequence of
operations for computation of wavelet amplitude coefficients is as
follows:(i)
represent the given functionf(x) by a discrete sequence fr,
where r=0,1,…,N−1;(ii)
compute the set of frequency coefficients by
fast Fourier transformf^m,
where m=0,1,…,N−1;(iii)
the inverse fast Fourier transform of the
octave blocksf^m generates the amplitudes of the harmonic
wavelet expansion of the function fr.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
ofx.(ii)
The analysed function is periodic inx 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 isf(x)=2sin2πx and f(x)=2cos2πx.
Following the algorithm, we discretise the interval [0;1] into N=2j equally spaced nods, and obtain discrete set
of values of functions(3.1)fr=2sin2πrN,fr=2cos2πrN,r=0,…,N−1.The fast Fourier transform (2.8) of the obtained
discrete sequence gives us the set Fourier coefficients f^m.
Recalling that a0=f^0,a1=f^1,
and aN/2=f^N/2,
we can easily find these three coefficients. Another part of coefficients from a2j to a2j+1−1 is obtained by computation of the inverse fast
s Fourier transform
(2.13) of coefficients from f^2j to f^2j+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 coefficientsa2j+k represents the discrete Fourier transform of
the reconstructed function fr.
Then, taking into account the shifting property (2.9), we can find f as inverse fast Fourier transform of f^.The results of decomposition and reconstruction of
functionsf(x)=2sin2πx and f(x)=2cos2πx are presented in Figures 2 and
3.Arbitrary given function: (a)sin2πx,
(b) cos2πx (dashed line), and its reconstructed clone
(solid line) from wavelet coefficients for N=8.
(a)(b)Arbitrary given function: (a)sin2πx,
(b) cos2πx (dashed line), and its reconstructed clone
(solid line) from wavelet coefficients for N=16.
(a)(b)One can notice that the plots of the reconstructed
functions are defined within the interval fromr=1 to r=N.
The difference between the algorithm and its corresponding computer code consists
in that we put
a1
in the code instead of a0,
and so forth . Therefore, the reconstruction of the function begins from point 1/N to 1, and not from 0 to N−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(3.2)ϵN=|f(xr)−frec(xr)|,r=0,…,N−1,where frec(xr) is the value of the reconstructed point. The
dependence of  absolute error of the  reconstruction of the function from lnN is represented in Figure 5 and for two partial
cases, when N=8 and N=16 can be found in Figure 4. As we can see, small
numbers of the level of decomposition j give a very good approximation, when we
reconstruct the function.Figure 4
Absolute error of the reconstruction off(x)=2sin2πx for N=8 (solid line) and N=16 (dashed line).Figure 5
Absolute error of the reconstruction
of2sin2πx 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 lnN and the corresponding approximation was
obtained in our simulations (see Figure 5). It means that the
increase of the length of N leads us to a slow, but steady
increase of the approximation error. The line of
the dependence of the error from N 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 N takes only integer values of 2j.The only disadvantage of harmonic wavelets is that its
decay rate is relatively low (proportional tox−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
thex-domain? How it can be used for the solution
of integrodifferential equations, and many others. This work is in progress.



---

*Source: 687318-2008-10-07.xml*