Derivation of the Daubechies D4 wavelet coefficients

The basic idea:
we want to find a function phi(x); that satisfies the recursive "scaling" relationship
        phi(t) = sum(h[k]*phi(2*t-k),k = 0 .. 3);
where the h[k]; are constants that satisfy a set of constraints that we choose.
(These constraints can be chosen so that the function phi(x); is "nice".)

Ingrid Daubechies imposed the following constraints:
1.  The even coefficients should total 1
2.  The odd coefficients should total 1
3.  The inner product of phi(t); with itself should total 2.
     Using the scaling relationship above, we findsum(h[k]^2,k = 0 .. 3) = 2;  
4.  sum(h[k] * h[k+2] , k) = 0
5. An extra "moment condition":   sum(k*(-1)^k*h[k],k) = 0;
> restart;
> solve({
>    	h0 + h2 = 1,
>    	h1 + h3 = 1,
>    	h0^2 + h1^2 + h2^2 + h3^2 = 2,
>    	h0*h2 + h1*h3 = 0,
>    	0*h0 - 1*h1 + 2*h2 - 3*h3 = 0
> });
> 

  {h1 = - 1/2 %1 + 1, h3 = 1/2 %1, h2 = 1/2 %1 + 1/2,

        h0 = - 1/2 %1 + 1/2}

                   2
  %1 := RootOf(2 _Z  - 2 _Z - 1)

> # convert the ``RootOf'' values above to explicit solutions
> IngridSolutions := allvalues([ % ]);
> 

  IngridSolutions := [{h1 = 3/4 - 1/4 sqrt(3), h3 = 1/4 + 1/4 sqrt(3),

        h2 = 3/4 + 1/4 sqrt(3), h0 = 1/4 - 1/4 sqrt(3)}], [{

        h1 = 3/4 + 1/4 sqrt(3), h3 = 1/4 - 1/4 sqrt(3),

        h2 = 3/4 - 1/4 sqrt(3), h0 = 1/4 + 1/4 sqrt(3)}]

> 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;
> 


           hvalue := i -> subs(op(IngridSolutions[2]), h.i)


                         [1/4 + 1/4 sqrt(3)]
                         [                 ]
                         [3/4 + 1/4 sqrt(3)]
                         [                 ]
                         [3/4 - 1/4 sqrt(3)]
                         [                 ]
                         [1/4 - 1/4 sqrt(3)]


                        h := array(0 .. 3, [])


                         h[0] := .6830127020


                         h[1] := 1.183012702


                         h[2] := .3169872980


                         h[3] := -.1830127020

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

                                M := 4


                             dilate := 3


                            points := 128


                               T := 383


                           iterations := 5


                           printlevel := 0


  phi := proc(i, t)
option operator, arrow;
    if 0 <= t and t < T then RETURN(tphi[i][t])
    else RETURN(0)
    fi
end


                      tphi := array(0 .. 5, [])


                                i := 0


                                i := 5


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

The Wavelet function psi(t); (also sometimes called W(t)) is then defined in terms of phi(t); as follows:

psi(t) = sum((-1)^k*h[M-1-k]*phi(2*t-k),k = 0 .. M-1); 
> 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)` );
> 

                       W := array(0 .. 383, [])


                                i := 5


There are other wavelets possible!!
For example:
... a rational wavelet:  {h[0]=3/5, h[1]=6/5, h[2]=2/5, h[3]=-1/5}
... another wavelet:  {h[0] = 5/sqrt(578), h[1]=20/sqrt(578), h[2]=12/sqrt(578), h[3]=-3/sqrt(578)}

> h := map(evalf,
>      array(0..3,[ 3/5, 6/5, 2/5, -1/5 ])
>      );
> 

  h := array(0 .. 3, [

        (0) = .6000000000

        (1) = 1.200000000

        (2) = .4000000000

        (3) = -.2000000000

        ])

> h := map(evalf,
>      array(0..3,[ (5)/sqrt(578), (20)/sqrt(578),
>                  (12)/sqrt(578), (-3)/sqrt(578) ])
>      );
> 

  h := array(0 .. 3, [

        (0) = .2079725826

        (1) = .8318903306

        (2) = .4991341984

        (3) = -.1247835496

        ])

>