prob1.html

Derivation of the Daubechies D4 wavelet coefficients

>

>    restart;
solve({
    h0 + h2 = 2,
    h1 + h3 = 1,
    h0^2 + h1^2 + h2^2 + h3^2 = 5,
    h0*h2 + h1*h3 = 0,
    0*h0 - 1*h1 + 2*h2 - 3*h3 = 0
});

${h3 = 1/2*RootOf(2*_Z^2-4*_Z-3,label = _L1), h2 = 1/2*RootOf(2*_Z^2-4*_Z-3,label = _L1)+1/2, h1 = -1/2*RootOf(2*_Z^2-4*_Z-3,label = _L1)+1, h0 = -1/2*RootOf(2*_Z^2-4*_Z-3,label = _L1)+3/2}$
${h3 = 1/2*RootOf(2*_Z^2-4*_Z-3,label = _L1), h2 = 1/2*RootOf(2*_Z^2-4*_Z-3,label = _L1)+1/2, h1 = -1/2*RootOf(2*_Z^2-4*_Z-3,label = _L1)+1, h0 = -1/2*RootOf(2*_Z^2-4*_Z-3,label = _L1)+3/2}$
${h3 = 1/2*RootOf(2*_Z^2-4*_Z-3,label = _L1), h2 = 1/2*RootOf(2*_Z^2-4*_Z-3,label = _L1)+1/2, h1 = -1/2*RootOf(2*_Z^2-4*_Z-3,label = _L1)+1, h0 = -1/2*RootOf(2*_Z^2-4*_Z-3,label = _L1)+3/2}$

> # convert the ``RootOf'' values above to explicit solutions
### WARNING: allvalues now returns a list of symbolic values instead of a sequence of lists of numeric values
IngridSolutions := allvalues([ % ]);

$IngridSolutions := [{h3 = 1/2+1/4*10^(1/2), h2 = 1+1/4*10^(1/2), h1 = 1/2-1/4*10^(1/2), h0 = 1-1/4*10^(1/2)}], [{h3 = 1/2-1/4*10^(1/2), h2 = 1-1/4*10^(1/2), h1 = 1/2+1/4*10^(1/2), h0 = 1+1/4*10^(1/2)}]...$
$IngridSolutions := [{h3 = 1/2+1/4*10^(1/2), h2 = 1+1/4*10^(1/2), h1 = 1/2-1/4*10^(1/2), h0 = 1-1/4*10^(1/2)}], [{h3 = 1/2-1/4*10^(1/2), h2 = 1-1/4*10^(1/2), h1 = 1/2+1/4*10^(1/2), h0 = 1+1/4*10^(1/2)}]...$

> hvalue := i -> subs( op(IngridSolutions[2]), h||i);
matrix(4,1, i -> hvalue(i-1) );
h := array(0..3);
for i_ from 0 to 3 do h[i_] := evalf(hvalue(i_)) od;

>    # Compute wavelets by iteration
# of the "Mother function" phi's dilation equation.

# Define wavelet of order M = 4:

M := 4;
dilate := (M-1);
points := 128;
T := dilate*points - 1;

#
# The dilation equation should hold for all values of t:
#
#      phi(t) = sum_{k=0}^{M-1} h_k phi(2t - k)
#
# We force this equation to hold with a kind of `Newton's Method'':
# starting with phi[0](t) = 1, we iteratively compute the dilation
#      phi[i](t) = sum_{k=0}^{M-1} h_k phi[i-1](2t - k).
# Below we iterate five times, but better results come from doing more.
#

iterations := 5;
printlevel := 0;

phi :=
(i,t) -> if (0<=t and t<T) then RETURN(tphi[i][t]) else RETURN(0) fi;

tphi := array(0..iterations);
for i from 0 to iterations do tphi[i] := array(0..T) od;

i := 0; # current iteration
for ij from 0 to T do tphi[i][ij] := 1.0 od; # phi = 1 initially

for i from 1 to iterations do
    for t from 0 to T do
        tphi[i][t] := sum( 'h[k] * phi(i-1,(2*t - k*points))',
                           'k'=0..(M-1) )
    od
od;

i := iterations;
plot( [[j, tphi[i][j]] $ j=0..T], title=`mother function phi(t)` );
# plot the final value of phi(t)

> plots[animate]( (p,q) -> tphi[q][p],
0..T, 0..iterations, frames=iterations, color=blue );

The Wavelet function (also sometimes called ) is then defined in terms of as follows:

>    W := array(0..T);

i := iterations;       # compute psi using the final value of phi
for t from 0 to T do
    W[t] := sum( '(-1)^(M-k) * h[(M-1)-k] * phi(i,(2*t - k*points))',
                 'k'=0..(M-1) )
od;

plot( [[j, W[j]] $ j=0..T], title=`wavelet function psi(t)` );

>