Version 15 of wavelet

Updated 2013-03-30 11:12:01 by AviationPete

The Wikipedia's description [L1 ], while typically accurate, is too abstract for CL's taste. Here's the idea: 'know Laplace or Fourier? 'Seen "sonograms" that resolve, say, a voice signal into its constituent frequencies? Perhaps you have some notion of the ubiquity and potency of such techniques. They rely on the duality of a signal (typically conceived along a time dimension) and its frequency analysis (the signal decomposed into the "harmonics" that make it up).

The signal and its constituent frequencies are ideally "infinite", that is, defined over the entire real (or complex, ...) line. Instead of those pure tones, consider instead signals limited in time--a note that starts and stops. Doesn't that hint at a better match for physical reality? And what if the calculus of signal decomposition into band- and time-limited resolvents were tweaked so that its calculations matched up particularly well with computer capabilities, by, for example, admitting natural expressions in binary? At that point, you'd have powerful, very high-performance techniques for analysis of signals in terms of physically-meaningful fundamental bases.

You'd have wavelets.

As wavelet theory remains young, there's a lot of exploration and interfacing still to do. Tcl, of course, makes a great partner.


aricb: I'm familiar with spectrograms, which are commonly used in phonetics to analyze speech signals (and which the Snack extension can produce). My understanding of spectrograms is that they are derived from Fourier analysis of a signal. Can somebody explain how a wavelet differs from a spectrogram?


Serge Kazantzev Time-scale (wavelet transforms) analysis is often presented as an alternative method to time-frequency (fourier transforms) analysis, and strong links exist between the two. On-line tutorials such as [L2 ] offer an exploration the differences and similarities between the two.

Wavelet spectogram called scalogram, refers to a time-scale energy distribution defined as the square modulus of the wavelet coefficients. Scalogram tastes like spectogram, however there is no trivial transition between scalogram and spectrogram.


TV Well, it's a complicated subject of course, preferably for a background in functional integrals. Theoretical physics already for about a century works with a (complicated form of) wavelett type functions, as a basis for Fock space which can be thought to contain all possible solution of the Schroedinger equation for a certain measurement basis. So the idea is most certainly very not new or a more complicated continuation of the top of mathematics. - Lars H: Are you sure that Fock space thingie isn't just plain old Hilbert space theory? The idea of having, even constructing, an orthogonal basis for the space of all solutions to some PDE is certainly that old (second half 19th century, I'd say, even though back then linear algebra hadn't been recognised yet, so people didn't know that this was what they were doing), but wavelets are more than that. That the entire basis can be generated from one function by scaling and shifting is rather remarkable (first known example was the Haar wavelet), and that it can furthermore be continuous or even differentiable has not been known for that many decades.

The comparison with Fourier analysis is that the point of fourier analysis is a frequency tranform over the whole of the time of the signal, in principle from -inf (or 0) to inf, or everywhere where the signal is non-zero. Practically, for instance with the mathematical often draconic FFT one usually wants to also localize the occurence of certain frequencies in time, so that the question arises how to say where a certain spectral line appears in the signal, and where not.

This is not so easily generally answerable, for two main reasons: it requires a non-trivial analysis of the signal, and the fourier transform doesn't answer this question at all. Mainly, the help of complicated (given.. with functional integrals) statistics can make analysis of every bit of the signal (like in physics) so that a fourier transform could be made for every point of the time axis, and then also with all kinds of widths of the gaussian which then appears. So the outcome of the transform becomes higher dimensional quite a bit. The integrals for standard functions which are computable for well known Fourier basic functions are not so computable anymore either, at least that becomes a lot more complicated.

In the more popular sense of the Fast Fourier Transform one takes a crude approcimation of a time segment based fourier transform, which is very non-trivial to analyse (try any FFT on a single sine wave which isn't of equivalent multiple frequency of the FFT interval, you'll get a well filled spectrum while there's just one spectral component..) , which can of course be practical. To solve the 'hard boundary' and also that the nature of that transform isn't like signal correlation analysis type that fits mathematically formulated problems better, the wavelet is sort of a completely smooth basis instead of the hard cut FFT and the full length fourier transform.

I guess wavelets can be infinite like an exponential multiplied by a sine, or (which has different fundamental mathematical properties) can be time interval limited, but smooth, like some polynomial coming from and going to 1 and then again zero, multiplied with a sine.

When an FFT result is achieved by average over for instance a few low level FFTs with overlapping intervals, all kinds of filtered results are achieved, which in general are not a good localized frequency analysis, but a bit smoother.

Ingrid Daubechies (IIRC) has done well known work on practical wavelet analysis, I think she was at Princeton at some point [L3 ].


Serge Kazantzev: The aim of the example below is to give a beginner some understanding of the wavelet transform using a practical application: image compression. We submitted the following image

(Image here is 404) http://yawtcl.free.fr/wiki/orig-pict.jpg

to a square wavelet transform following Stephane Mallat's decomposition

(Image here is 404) http://yawtcl.free.fr/wiki/mallat.jpg

based on the CRF(13,7) for which the lowpass and the highpass step filter have the following transfer functions:

(Image here is 404) http://yawtcl.free.fr/wiki/crf13-7.jpg

The image being a set of 2D data, we apply the transform to the columns and to the rows and create four smaller image-like sub-regions: the subbands.

(Image here is 404) http://yawtcl.free.fr/wiki/xform-crf.jpg

The upper left subband (subband LL) contains the averages of the pixel values, whilst the 3 other subbands (LH, HL and HH) contains the wavelet coefficients that we need to reconstruct the image. Another way to say is that the upper left quarter of the image is a lower resolution version of the original image and the 3 other quarters contain the details of the original image. To go back to the original, the inverse transform will have to put the details back into the small image. Looking now at the histogram and

http://yawtcl.free.fr/wiki/xcrf-histo.jpg

at the wavelet coefficients for the image

http://yawtcl.free.fr/wiki/xcrf-coef2.jpg

the temptation is high to quantize the coefficients in order to compress the image. After quantization (i.e throwing away coefficients that are 0), and inverse transform, we obtain the following image:

http://yawtcl.free.fr/wiki/xinv-pict.jpg

with the wavelet coefficients becoming:

http://yawtcl.free.fr/wiki/xt10-coef2.jpg

for a compression ratio of 63.7%.

Nb: compression ratio here means the number of zero coefficients divided by the total number of coefficients expressed as a %.


Peter K: I found no sample code, so I decided to try my luck with something based on the excellent work of Ian Kaplan, Bear Products International, www.bearcave.com, 2001. It is meant to filter out high frequencies and comes with no warranty.

#
proc Einrichten {} {
#
# Einrichten des Hauptfensters
#
    global Canv
    global Zaehler
    set Zaehler 1
#
    toplevel .wavelet -class Toplevel -cursor left_ptr
    wm focusmodel .wavelet passive
    wm overrideredirect .wavelet 0
    wm resizable .wavelet 1 1
    wm geometry .wavelet 1280x720+20+40
    wm title .wavelet "Wavelet"
    wm deiconify .wavelet
#
    frame .wavelet.oben
    frame .wavelet.unten -height 60 -background systemDialogBackgroundActive
    set Canv [canvas .wavelet.oben.plotRect -background LightSkyBlue \
        -borderwidth 3 -highlightthickness 3]
    ttk::spinbox .wavelet.unten.zahl -values [list 1 2 4 8 16 32 64]
    .wavelet.unten.zahl set 64
    ttk::button  .wavelet.unten.rechnen   -text "Rechnen"   -width 10
    ttk::button  .wavelet.unten.cancelBut -text "Abbrechen" -width 10 \
        -command exit
    ttk::button  .wavelet.unten.fertigBut -text "Fertig"    -width 10 \
        -default active -command exit
#
    pack .wavelet.oben.plotRect -fill both -expand yes
    grid rowconfigure .wavelet { 0 } -weight 5 -pad 0
    grid rowconfigure .wavelet { 1 } -weight 0 -minsize 60
    grid columnconfigure .wavelet { 0 } -weight 1 -pad 0
    grid .wavelet.oben  -column 0 -row 0 -sticky nsew
    grid .wavelet.unten -column 0 -row 1 -sticky nsew
#
    grid columnconfigure .wavelet.unten  0      -weight 0 -minsize 18
    grid columnconfigure .wavelet.unten {1 2 3} -weight 1 -minsize 18
    grid .wavelet.unten.zahl      -column 0 -row 0 -padx 18 -sticky ew
    grid .wavelet.unten.rechnen   -column 1 -row 0 -padx 18 -sticky ew
    grid .wavelet.unten.cancelBut -column 2 -row 0 -padx 12 -sticky ew
    grid .wavelet.unten.fertigBut -column 3 -row 0 -padx 18 -sticky ew
#
    raise .wavelet
}
#
############################################################
#
proc Plotten { Y n } {
#
    global Farbe
    global Canv
#
    $Canv delete Linie$n
    set Liste [list ]
#    set Grenze [expr {[.wavelet2 * .unten.zahl get]}]
    set Grenze [llength $Y]
    set X  50
    set dX [expr {16 * 128 / $Grenze}]
    set i  1
    while {$i < $Grenze} {
        lappend Liste $X
        lappend Liste [expr {640 - 3 * [lindex $Y $i]}]
        incr i 2
        incr X $dX
    }
    if {$n >= [array size Farbe]} {
        set n [expr {$n - [array size Farbe]}]
    }
    $Canv create line $Liste -fill $Farbe($n) -width 3 -tags Linie$n
}
#
############################################################
#
proc Rechnen { Y } {
#
    global Farbe
    global Canv
    global Zaehler
#
    Haar::forwardTrans $Y
    set Grenze [.wavelet.unten.zahl get]
#
# Hier gezielt Koeffizienten nullen, wenn gefiltert werden soll
#
    Haar::filter $Grenze
#
    Plotten [Haar::inverseTrans $Grenze] $Zaehler
    incr Zaehler
}
#
############################################################
#
proc main {} {
#
    global Farbe
    global Canv
#
    set Y [ list \
     0  77.6875  1  78.1875  2  82.0625  3  85.5625  4  86.7500  5  82.4375 \
     6  82.2500  7  82.7500  8  81.2500  9  79.5625 10  80.2813 11  79.8750 \
    12  77.7500 13  74.7500 14  78.5000 15  79.1875 16  78.8125 17  80.3125 \
    18  80.1250 19  79.3125 20  83.7500 21  89.8125 22  87.7500 23  91.1250 \
    24  94.4375 25  92.7500 26  98.0000 27  97.1875 28  99.4375 29 101.7500 \
    30 108.5000 31 109.0000 32 105.2500 33 105.5000 34 110.0000 35 107.0000 \
    36 107.2500 37 103.3125 38 102.8750 39 102.4375 40 102.0000 41 101.3125 \
    42  97.4375 43 100.5000 44 107.7500 45 110.2500 46 114.3125 47 111.2500 \
    48 114.8125 49 112.6875 50 109.4375 51 108.0625 52 104.5625 53 103.2500 \
    54 110.5625 55 110.7500 56 116.3125 57 123.6250 58 120.9375 59 121.6250 \
    60 127.6875 61 126.0625 62 126.3750 63 124.3750 ]
    array set Farbe [ list 0 black 1 blue 2 green 3 red ]
#
    Einrichten
    .wavelet.unten.rechnen configure -command [list Rechnen $Y]
#
    Plotten $Y 0
}
#
############################################################
#
# Tiefpassfilter mit modischen Wavelets
#
#     This software was written and is copyrighted by Ian Kaplan, Bear
#     Products International, www.bearcave.com, 2001.
#
# allerdings in Java, dies hier ist Tcl. Und die Arbeit von Herrn Kaplan
# ist ein schöner Einstieg, aber mehr auch nicht.
#
# Um alles so einfach wie sinnvoll zu halten, werden nur die sehr simplen
# Haar-Wavelets verwendet. Die sind so simpel, daß die Wavelet-Bezeichnung
# schon fast irreführend ist.
#
package provide Haar 1.0

namespace eval Haar {
    variable forward 1
    variable inverse 2
    variable vec
    variable avg
    variable dif
    variable logN
}
#
############################################################
#
proc Haar::predict { N direction } {
#
# Haar predict step
#
    variable forward
    variable inverse
    variable vec
#
    set half [expr {$N >> 1}]
#
    for {set i 0} {$i < $half} {incr i} {
        set predictVal $vec($i)
        set j [expr {$i + $half}]

        if {$direction == $forward} {
            set vec($j) [expr {$vec($j) - $predictVal}]
        } elseif {$direction == $inverse} {
            set vec($j) [expr {$vec($j) + $predictVal}]
        } else {
            tk_messageBox -icon warning -type ok \
              -title "Fehler in Funktion Haar::predict" \
              -message "Bad direction value $direction"
        }
    }
}
#
############################################################
#
proc Haar::update { N direction } {
#
# Update step of the Haar wavelet transform.
#
    variable forward
    variable inverse
    variable vec
#
    set half [expr {$N >> 1}]
#
    for {set i 0} {$i < $half} {incr i} {
        set j [expr {$i + $half}]
        set updateVal [expr {$vec($j) / 2.0}]

        if {$direction == $forward} {
            set vec($i) [expr {$vec($i) - $updateVal}]
        } elseif {$direction == $inverse} {
            set vec($i) [expr {$vec($i) + $updateVal}]
        } else {
            tk_messageBox -icon warning -type ok \
              -title "Fehler in Funktion Haar::update" \
              -message "Bad direction value $direction"
        }
    }
}
#
############################################################
#
proc Haar::split { N } {
#
#     Split the vec into even and odd elements,
#     where the even elements are in the first half
#     of the vector and the odd elements are in the
#     second half.
#
    variable vec
#
    set start  1
    set end  [expr {$N - 1}]
#
    while {$start < $end} {
        for {set i $start} {$i < $end} {incr i} {
           set tmp $vec($i)
           set vec($i) $vec([incr i])
           set vec($i) $tmp
        }
        incr start
        incr end -1
    }
}
#
############################################################
#
proc Haar::merge { N } {
#
#    Merge the odd elements from the second half of the N element
#    region in the array with the even elements in the first
#    half of the N element region. The result will be the
#    combination of the odd and even elements in a region
#    of length N.
#
    variable vec
#
    set half  [expr {$N >> 1}]
    set start [expr {$half - 1}]
    set end  $half
#
    while {$start > 0} {
        for {set i $start} {$i < $end} {incr i} {
            set tmp $vec($i)
            set vec($i) $vec([incr i])
            set vec($i) $tmp
        }
        incr start -1
        incr end
    }
}
#
############################################################
#
proc Haar::forwardTrans { Y } {
#
    variable forward
    variable inverse
    variable vec
    variable avg
    variable dif
    variable logN
#
    set N [expr {[llength $Y] / 2}]
    array set vec $Y
#
    set i 0
    for {set n 0} {$n < $N} {incr n 2} {
        set m [expr {$n + 1}]
        set avg($i,0) [expr {($vec($n) + $vec($m)) / 2}]
        set dif($i,0) [expr {($vec($n) - $vec($m)) / 2}]
        incr i
    }
    set max [expr {$N >> 1}]
    set j   1
    while {$max > 1} {
        set i 0
        set j1 [expr {$j - 1}]
        for {set n 0} {$n < $max} {incr n 2} {
            set m [expr {$n + 1}]
            set avg($i,$j) [expr {($avg($n,$j1) + $avg($m,$j1)) / 2}]
            set dif($i,$j) [expr {($avg($n,$j1) - $avg($m,$j1)) / 2}]
            incr i
        }
        set max [expr {$max >> 1}]
        incr j
    }
    set logN [expr {$j - 1}]
}
#
############################################################
#
proc Haar::filter { n } {
#
    variable forward
    variable inverse
    variable vec
    variable dif
    variable logN
#
    set j [expr int($logN - [ld $n] + 1)]
    set max $n
    while {$j >= 0} {
        for {set i 0} {$i < $max} {incr i} {
            set dif($i,$j) 0.0
        }
        incr j -1
        set max [expr {$max << 1}]
    }
}
#
############################################################
#
proc Haar::inverseTrans { Grenze } {
#
    variable forward
    variable inverse
    variable vec
    variable avg
    variable dif
    variable logN
#
# Der Trick hier: Resampling nur bis zur gewünschten Auflösung.
# Das Nullen der Koeffizienten und volles Resampling bedeutet, daß man die
# hohe zeitliche Auflösung regeneriert, ohne die hohen Frequenzen wiedergeben
# zu wollen. In Folge gäbe es lange Plateaus mit sehr steilen Gradienten.
#
    set N [array size vec]
    puts stdout "inverseTrans : N = $N,  logN = $logN"
#
    set j1 [expr {$logN - 1}]
    set loc(0,$j1) [expr {$avg(0,$logN) + $dif(0,$logN)}]
    set loc(1,$j1) [expr {$avg(0,$logN) - $dif(0,$logN)}]
    set max 4
    while {$j1 > 0} {
        set  i  0
        set  j  $j1
        incr j1 -1
###        puts stdout "Runde $j1 mit max = $max, Grenze = $Grenze"
        for {set n 0} {$n < $max} {incr n 2} {
            set m [expr {$n + 1}]
            set loc($n,$j1) [expr {$loc($i,$j) + $dif($i,$j)}]
            set loc($m,$j1) [expr {$loc($i,$j) - $dif($i,$j)}]
            incr i
        }
        set max [expr {$max << 1}]
        if {$max == $Grenze} break
    }
    set  i  0
    set  j  $j1
###    puts stdout "Runde $j1 mit max = $Grenze"
    for {set n 0} {$n < $Grenze} {incr n 2} {
        set m [expr {$n + 1}]
        set vec($n) [expr {$loc($i,$j) + $dif($i,$j)}]
        set vec($m) [expr {$loc($i,$j) - $dif($i,$j)}]
###        puts stdout "dif($i,$j) = $dif($i,$j) => vec($n) = $vec($n)   vec($m) = $vec($m)"
###        flush stdout
        incr i
    }
#
    set Y [list ]
    for {set n 0} {$n < $Grenze} {incr n} {
        lappend Y $n $vec($n)
    }
    return $Y
}
#
############################################################
#
# Logarithmus zur Basis 2 (wiki.tcl.tk/819)
#
proc ld x "expr {log(\$x)/[expr log(2)]}"
#
############################################################
#
console show
main