NAG FL Interface
c09baf (dim1_​cont)

1 Purpose

c09baf computes the real, continuous wavelet transform in one dimension.

2 Specification

Fortran Interface
Subroutine c09baf ( wavnam, wparam, n, x, nscal, scales, c, ifail)
Integer, Intent (In) :: wparam, n, nscal, scales(nscal)
Integer, Intent (Inout) :: ifail
Real (Kind=nag_wp), Intent (In) :: x(n)
Real (Kind=nag_wp), Intent (Out) :: c(nscal,n)
Character (*), Intent (In) :: wavnam
C Header Interface
#include <nag.h>
void  c09baf_ (const char *wavnam, const Integer *wparam, const Integer *n, const double x[], const Integer *nscal, const Integer scales[], double c[], Integer *ifail, const Charlen length_wavnam)
The routine may be called by the names c09baf or nagf_wav_dim1_cont.

3 Description

c09baf computes the real part of the one-dimensional, continuous wavelet transform
Cs,k = xt 1 s ψ* t-k s dt ,  
of a signal xt at scale s and position k, where the signal is sampled discretely at n equidistant points xi, for i=1,2,,n. ψ is the wavelet function, which can be chosen to be the Morlet wavelet, the derivatives of a Gaussian or the Mexican hat wavelet (* denotes the complex conjugate). The integrals of the scaled, shifted wavelet function are approximated and the convolution is then computed.
The mother wavelets supplied for use with this routine are defined as follows.
  1. 1.The Morlet wavelet (real part) with nondimensional wave number κ is
    ψx = 1 π 1/4 cosκx - e -κ2/2 e -x2/2 ,  
    where the correction term, e-κ2/2 (required to satisfy the admissibility condition) is included.
  2. 2.The derivatives of a Gaussian are obtained from
    ψ^m x = dm e -x2 d xm ,  
    taking m=1,,8. These are the Hermite polynomials multiplied by the Gaussian. The sign is then adjusted to give ψ^m0>0 when m is even while the sign of the succeeding odd derivative, ψ^m+1, is made consistent with the preceding even numbered derivative. They are normalized by the L2-norm,
    pm = - ψ^ m x 2 d x 1/2  
    The resulting normalized derivatives can be written in terms of the Hermite polynomials, Hmx, as
    ψ m x = α Hmx e -x2 pm ,  
    where
    α = 1, when ​m=0,3 mod 4; -1, when ​m=1,2 mod 4.  
    Thus, the derivatives of a Gaussian provided here are,
    ψ1 x = - 2π 1/4 2 x e -x2 ,  
    ψ2 x = - 2π 1/4 1 3 4x2 - 2 e -x2 ,  
    ψ3 x = 2π 1/4 115 8x3 - 12x e -x2 ,  
    ψ4 x = 2π 1/4 1105 16x4 - 48x2 + 12 e -x2 ,  
    ψ5 x = - 2π 1/4 1 3105 32x5 - 160x3 + 120x e -x2 ,  
    ψ6 x = - 2π 1/4 1 31155 64x6 - 480x4 + 720x2 - 120 e -x2 ,  
    ψ7 x = 2π 1/4 1 315015 128x7 - 1344x5 + 3360x3 - 1680x e -x2 ,  
    ψ8 x = 2π 1/4 1 451001 256x8 - 3584x6 + 13440x4 - 13440x2 + 1680 e -x2 .  
  3. 3.The second derivative of a Gaussian is known as the Mexican hat wavelet and is supplied as an additional function in the form
    ψx = 2 3 π1/4 1-x2 e -x2/2 .  
    The remaining normalized derivatives of a Gaussian can be expressed as multiples of the exponential e - t2 / 2 by applying the substitution x = t / 2 followed by multiplication with the scaling factor, 1 / 24 .

4 References

Daubechies I (1992) Ten Lectures on Wavelets SIAM, Philadelphia

5 Arguments

1: wavnam Character(*) Input
On entry: the name of the mother wavelet. See the C09 Chapter Introduction for details.
wavnam='MORLET'
Morlet wavelet.
wavnam='DGAUSS'
Derivative of a Gaussian wavelet.
wavnam='MEXHAT'
Mexican hat wavelet.
Constraint: wavnam='MORLET', 'DGAUSS' or 'MEXHAT'.
2: wparam Integer Input
On entry: the nondimensional wave number for the Morlet wavelet or the order of the derivative for the Gaussian wavelet. It is not referenced when wavnam='MEXHAT'.
Constraints:
  • if wavnam='MORLET', 5wparam20;
  • if wavnam='DGAUSS', 1wparam8.
3: n Integer Input
On entry: the size, n, of the input dataset x.
Constraint: n2.
4: xn Real (Kind=nag_wp) array Input
On entry: x contains the input dataset xj=xj, for j=1,2,,n.
5: nscal Integer Input
On entry: the dimension of the array scales and the first dimension of the array c as declared in the (sub)program from which c09baf is called. The number of scales to be computed.
Constraint: nscal1.
6: scalesnscal Integer array Input
On entry: the scales at which the transform is to be computed.
Constraint: scalesi1, for i=1,2,,nscal.
7: cnscaln Real (Kind=nag_wp) array Output
On exit: the transform coefficients at the requested scales, where cij is the transform coefficient Ci,j at scale i and position j.
8: ifail Integer Input/Output
On entry: ifail must be set to 0, -1 or 1. If you are unfamiliar with this argument you should refer to Section 4 in the Introduction to the NAG Library FL Interface for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value -1 or 1 is recommended. If the output of error messages is undesirable, then the value 1 is recommended. Otherwise, if you are not familiar with this argument, the recommended value is 0. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry ifail=0 or -1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
ifail=1
On entry, wavnam=value was an illegal value.
ifail=2
On entry, wavnam='MORLET' and wparam=value.
Constraint: if wavnam='MORLET', 5wparam20.
On entry, wavnam='DGAUSS' and wparam=value.
Constraint: if wavnam='DGAUSS', 1wparam8.
ifail=3
On entry, n=value.
Constraint: n2.
ifail=5
On entry, nscal=value.
Constraint: nscal1.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

The accuracy of c09baf is determined by the fact that the convolution must be computed as a discrete approximation to the continuous form. The input signal, x, is taken to be piecewise constant using the supplied discrete values.

8 Parallelism and Performance

c09baf is not threaded in any implementation.

9 Further Comments

Workspace is internally allocated by c09baf. The total size of these arrays is 213 + n + nk - 1 real elements and nk integer elements, where nk = k × maxscalesi and k=17 when wavnam='MORLET' or 'DGAUSS' and k=11 when wavnam='MEXHAT'.

10 Example

This example computes the continuous wavelet transform of a dataset containing a single nonzero value representing an impulse. The Morlet wavelet is used with wave number κ=5 and scales 1, 2, 3, 4.

10.1 Program Text

Program Text (c09bafe.f90)

10.2 Program Data

Program Data (c09bafe.d)

10.3 Program Results

Program Results (c09bafe.r)