src/examples/wavelet.c

    Wavelet transforms and filtering

    This simple example illustrates how to compute wavelet transforms and perform filtering and mesh adaptation. We do this in one dimension, using bi-trees.

    #include "grid/bitree.h"

    A simple function to output fields on each level.

    void write_level (scalar * list, const char * tag, FILE * fp)
    {
      for (int l = 0; l <= depth(); l++)
        foreach_level (l, serial) {
          fprintf (fp, "%s%d %g ", tag, l, x);
          for (scalar s in list)
    	fprintf (fp, "%g ", s[]);
          fputc ('\n', fp);
        }
    }
    
    int main()
    {

    We consider a periodic 1D domain.

      init_grid (128);
      periodic (right);
      size (2);
      
      scalar s[], w[], s2[];

    We can optionally change the prolongation operator. We need to make sure that all fields use the same prolongation.

    #if 0
      for (scalar i in {s,w,s2})
        i.refine = i.prolongation = refine_biquadratic;
    #endif

    We initialise the field with a function containing low-frequency and localised high-frequency components.

      foreach()
        s[] = sin(2.*pi*x) + 0.4*sin(15*pi*x)*max(sin(2.*pi*x), 0);
      boundary ({s});

    The w field contains the wavelet transform of s.

      wavelet (s, w);

    We check that the inverse wavelet transform recovers the initial signal (almost) exactly.

      inverse_wavelet (s2, w);
      foreach()
        assert (fabs (s2[] - s[]) < 1e-12);
      write_level ({s,w}, "", stderr);

    The original signal and its wavelet transform.

    maxlevel = 6
    spacing = 1.5
    set ytics -maxlevel,1,0
    set xlabel 'x'
    set ylabel 'Level'
    set grid ytics ls 0
    max(a,b) = a > b ? a : b
    set samples 500
    f(x)=sin(2.*pi*x)+0.4*sin(15.*pi*x)*max(sin(2.*pi*x), 0);
    plot for [i=0:maxlevel] '< grep ^'.i.' log' u 2:($3/spacing - i) \
                               w lp t '' lt 7,			   \
         for [i=0:maxlevel] f(x)/spacing - i t '' lt 1
    Original signal representation on each level (script)

    Original signal representation on each level (script)

    plot for [i=0:maxlevel] '< grep ^'.i.' log' u 2:($4/spacing - i) \
                               w lp t '' lt 7,			   \
         for [i=0:maxlevel] f(x)/spacing - i t '' lt 1
    Corresponding wavelet decomposition (script)

    Corresponding wavelet decomposition (script)

    Localised low-pass filtering

    We want to filter the high-frequency components of the signal, but only locally i.e. for x > 1. To do this we cancel the high-frequency (i.e. high-level) wavelet coefficients and recover the corresponding filtered signal s2 with the inverse wavelet transform.

      for (int l = 5; l <= 6; l++) {
        foreach_level (l)
          w[] *= (x <= 1);
        boundary_level ({w}, l);
      }
      inverse_wavelet (s2, w);
      write_level ({s2,w}, "a", stdout);
    plot for [i=0:maxlevel] '< grep ^a'.i.' out' u 2:($3/spacing - i) \
                               w lp t '' lt 7,			   \
         for [i=0:maxlevel] f(x)/spacing - i t '' lt 1
    Low-pass filtered reconstruction (for x > 1) (script)

    Low-pass filtered reconstruction (for x > 1) (script)

    High-pass filtering

    Here we filter out the low-frequency components everywhere.

      wavelet (s, w);
      for (int l = 0; l < 5; l++) {
        foreach_level (l)
          w[] = 0.;
        boundary_level ({w}, l);
      }
      inverse_wavelet (s2, w);
      write_level ({s2,w}, "b", stdout);
    plot for [i=0:maxlevel] '< grep ^b'.i.' out' u 2:($3/spacing - i) \
                               w lp t '' lt 7,			   \
         for [i=0:maxlevel] f(x)/spacing - i t '' lt 1
    High-pass filtered reconstruction (script)

    High-pass filtered reconstruction (script)

    Mesh adaptation

    We can use the wavelet coefficients to decide where to coarsen the mesh.

      wavelet (s, w);
      unrefine (fabs(w[]) < 0.1);
      write_level ({s,w}, "c", stdout);
    }
    plot for [i=0:maxlevel-1] '< grep ^c'.i.' out' u 2:($3/spacing - i)	\
            w lp t '' lt 7,							\
          '< grep ^c6 out' u 2:($3/spacing - 6) w p t '' lt 7,		\
          for [i=0:maxlevel] f(x)/spacing - i t '' lt 1
    Adapted mesh (script)

    Adapted mesh (script)

    Discussion

    • Note that the restriction operators assume that the discrete quantity is the average of the function over the cell, not its point value. While this makes no difference for second-order (i.e. linear) prolongation operators, it matters for higher-order prolongation but is not taken into account in this example.
    • Boundary conditions for the wavelet transform are important and non-trivial. This is not an issue here since we use a periodic domain.
    • The link between the transform implemented here and formal wavelets is not entirely clear. It is probably a kind of second-generation wavelet.
    • The filters applied in the example could be formulated in a more generic way using a filtering function \phi(x,l) and a code looking like:
      foreach_cell()
        w[] *= phi(x,level);

    combined with proper application of boundary conditions.

    References

    [sweldens1998]

    Wim Sweldens. The lifting scheme: A construction of second generation wavelets. SIAM journal on mathematical analysis, 29(2):511–546, 1998.