1
#to show daubechies N=2 scaling function & mother wavelet

close all;
 
#coefficients of Daubechies N=2
p=[0.6830,1.1830,0.3170,-0.1830];
q=zeros(1,4);
q(1)=p(4);
q(2)=-p(3);
q(3)=p(2);
q(4)=-p(1);

phi_origin=[0,1.366,-0.366,0];
psi_origin=[0,-0.36603,-1.366,0];
 #resolution must be 2^n
n=input("  please input N, N is 2^N.\n  N=2 to 8 is useful.  N:");
 reso=2**n;
 my_size=reso*3+1;
printf("\n  daubechies[%g]\n\n",my_size);
phi=zeros(1,my_size); 
psi=zeros(1,my_size);

#initialize phi array
 phi(reso+1)=phi_origin(2);
 phi(reso*2+1)=phi_origin(3);
 
#sub function###############################################
function retval=phi_value(array,index,my_size)
                if( (index<=1) | (index>my_size) ) retval=0;
                else
                  retval=array(index);
                endif
 endfunction
############################################################
# calc phi function value
for k=1:n
    for i=1:my_size
      index=i*2-1;
     
      x=p(1)*phi_value(phi,index,my_size);
      x=x+p(2)*phi_value(phi,index-reso,my_size);
      x=x+p(3)*phi_value(phi,index-reso*2,my_size);
      x=x+p(4)*phi_value(phi,index-reso*3,my_size);

      phi(i)=x; 
      x=0;
      index=0;
    endfor
endfor

# calc psi value
for k=1:n
    for i=1:my_size
      index=i*2-1;
     
      x=q(1)*phi_value(phi,index,my_size);
      x=x+q(2)*phi_value(phi,index-reso,my_size);
      x=x+q(3)*phi_value(phi,index-reso*2,my_size);
      x=x+q(4)*phi_value(phi,index-reso*3,my_size);

      psi(i)=x; 
      x=0;
      index=0;
    endfor
endfor

# plot Daubechies functions & their spectrum

subplot(2,1,1);

t=linspace(0,3,my_size);
plot(t,phi,"b-;sacaling function;");
title("Daubechies N=2 functions");

hold on;
subplot(2,1,1);
t=linspace(-1,2,my_size);
plot(t,psi,"r-;mother wavelet;");



y=fft(phi,512);
x=abs(y);
max_value=max(x);
x=x/max_value;

subplot(2,1,2);
plot(x(1:128),"b-;scaling function;");
title("spectrum");
hold on;

y=fft(psi,512);
x=abs(y);
max_value=max(x);
x=x/max_value;
subplot(2,1,2);

plot(x(1:128),"r-;mother wavelet;");

