# Haar Filtering and Condensing
# Filter out high frequency noise
# 
# The goal here is to use Haar Wavelets to 
# filter out high frequency noise. First the signal 
# is defined by defining a function f. Then the 
# signal is discretized over the interval [0,1]
# with step size 2^N (N=8 or 9 is typical).
# Then the signal is decomposed in a wavlet
# decomposition. The various levels of decomposition
# represent the various frequency components of the signal
# from 0 to 2^N.
# 
# First define the signal.

> g:=x->x*exp(-x^2);

                                           2
                         g := x -> x exp(-x )

> f:=x->sin(2*Pi*x)+cos(4*Pi*x)+sin(8*Pi*x)+4*g((x-1/3)*64)+4*g((x-2/3)*
> 128);

f := x -> sin(2 Pi x) + cos(4 Pi x) + sin(8 Pi x) + 4 g(64 x - 64/3)

     + 4 g(128 x - 256/3)

> plot(f,0..1,numpoints=400,thickness=0);
# Bascially, the signal f  involves a smoothly varying term (given
# first) together with
# some high frequency  noise (the terms involving g) at x=1/3 and 2/3.
> 
> j:='j';

                                j := j

# Now define the value of N which represents the number 
# of levels - - the highest level has resolution 2^N
# 
> N:=8;

                                N := 8

# Now store the number of data points 2^j  that will
# be stored at the jth level where j ranges from 0 to N
# 
> for j from 0 to N do
> L[j]:=2^j;
> od:
> j:='j';

                                j := j

> l:='l';

                                l := l

# Now discretize the function on [0,1] into the highest
# Nth level.
# 
> for l from 0 to L[N] do
> a[N,l]:=evalf(f(l/L[N]));
> od:
> 
> a[8,3];

                             1.353025751

> evalf(f(3/2^8));

                             1.353025751

> j:='j';

                                j := j

> l:='l';

                                l := l

# 
# Now execute the decomposition algorithm.
# The projection of the signal onto the jth
# level (V_j) is given in the array a[j,*].
# The jth wavelet level, i.e. the projection
# of the signal onto W_j is represented
# by the array b[j,*]. Care must be taken
# to use the periodicity of the data.
# 
> for j from N by -1 to 1 do
> for l from 0 to L[j-1]-1 do
> a[j-1,l]:=(a[j,2*l]+a[j,1+2*l])/2;
> b[j-1,l]:=(a[j,2*l]-a[j,1+2*l])/2;
> od;
> a[j-1,L[j-1]]:=(a[j,L[j]]+a[j,1])/2;
> b[j-1,L[j-1]]:=(a[j,L[j]]-a[j,1])/2;
> od:
> b[7,5];

                             -.0249957890

# If the goal is data compression, then skip to that section.
# 
> j:='j'; l:='l';

                                j := j


                                l := l

# Now we plot the projection of the function onto
# the jth level (the array a[j,*]) by storing
# the x and y values into the array level[j,*]
# whose elements are of the form [x,y]
# with x=k/L[j] and y=a[j,k]
# 
> for j from 0 to N do
> for l from 1 to L[j] do
> level[j,2*l]:=[l/L[j],a[j,l-1]];
> level[j,2*l-1]:=[(l-1)/L[j],a[j,l-1]];
> od;
> 
> 
> od;
> j:='j'; l:='l';

                                j := j


                                l := l

# Now plot level j where you can choose which 
# level you wish from 0 to N
> for j from 0 to N do
> plotlevel[j]:=seq(level[j,l],l=1..2*L[j]);
> od:
> plot([plotlevel[8]],style=line, thickness=0);
> 
> j:='j'; l:='l';

                                j := j


                                l := l

# Reconstruction and Data Compression
# 
# Here we compress the signal.
# We define a tolerance tol and a variable
# count which will keep track of how
# many coefficients that are set equal to 
# zero.
# 
# 
# 
> j:='j'; l:='l'; tol:=.05; count:=0;

                                j := j


                                l := l


                              tol := .05


                              count := 0

> for j from 0 to N-1 do
> for l from 1 to L[j] do
> if abs(b[j,l])< tol then b[j,l]:=0; count:=count+1; fi;
> od:
> od:
> j:='j'; l:='l'; count;

                                j := j


                                l := l


                                 123

# Note that 2^N-count represents the number of NON-zero
# wavelet coefficients.
# 
# Now we reconstruct the signal using the (now modified)
# wavelet coefficients.
# 
# 
> 
> for j from 1 to N do
> for l from 1 to L[j-1] do
> a[j,2*l]:=a[j-1,l]+b[j-1,l];
> a[j,2*l-1]:=a[j-1,l-1]-b[j-1,l-1];
> od:
> a[j,0]:=a[j-1,0]+b[j-1,0];
> od:
> 
> 
> 
> 
> j:='j'; l:='l';

                                j := j


                                l := l

# Now we plot the compressed signal by forming
# an array of x,y values in rlevel[j,*]. This array
# is of the form [x,y] where x is the x-coordinate,
# i.e. k/L[j]) and whose y coordinate is a[j,k]
# where k ranges from 1 to L[j].
# Note that the same y value a[j,k]
# is used for 2 x-values - - one at either
# end of the interval from k/L[j] to (k+1)/L[j].
# 
> 
> for j from 0 to N do
> for l from 1 to L[j] do
> rlevel[j,2*l]:=[l/L[j],a[j,l-1]];
> rlevel[j,2*l-1]:=[(l-1)/L[j],a[j,l-1]];
> od:
> od:
> 
> j:='j'; l:='l';

                                j := j


                                l := l

> for j from 0 to N do
> plotrlevel[j]:=seq(rlevel[j,l],l=1..2*L[j]);
> od:
> plot([plotrlevel[N]],style=line, thickness=0);
> plot(f,0..1, thickness =0);
> j:='j'; l:='l';

                                j := j


                                l := l

> 
> 
> 
>