1
#close all;
clear
#1
##################################################################
#Φ(0:N)の値を求めるプログラム。 Nは整数
##################################################################
close all;

D=input("\n  ==calc Daubechies's phi[integers] values==\n　please input D. D is DaubechiesD??\n  D=4,8,20.  D:");


#Daubechies coefficients size D=N*2
p=zeros(1,D);

switch D
 case {4}
	 p=[0.6830,1.1830,0.3170,-0.1830];
 case {8}
	p=[0.32580343,1.01094572,0.8922014,-0.03957503,-0.26450717,0.0436163,0.0465036,-0.01498699];
 case {20}
	p=[0.03771716,0.26612218,0.74557507,0.97362811,\
	   0.39763774,-0.35333620,-0.27710988,0.18012745,\
	   0.13160299,-0.10096657,-0.04165925,0.04696981,5.10043697e-3,\
	   -0.01517900,1.97332536e-3,2.81768659e-3,-9.69947840e-4,-1.64709006e-4,1.32354367e-4,-1.875841e-5];
 otherwise
	error ("invalid value");
endswitch

a=zeros(D-2,D-2);
sub=ones(1,D-2);

#先ず、連立方程式を行列で表現する。(D-2)*(D-2)の行列を作る
for i=1:D-2
	for k=1:D-2

		if (i==k) retval=p(i+1)-1;
		else 
		   if(i*2-k+1<=D && (i*2-k+1>0) )
			retval=p(i*2-k+1);
		   else retval=0;
		   endif
		endif
		a(i,k)=retval;
	endfor
endfor   
# 最上段に１行の追加
 b=[sub;a];

#連立方程式を解く
 x=zeros(D-2,1);
 c=zeros(D-1,1);
 c(1)=1;
# b*x=c, solve
 x= b\c;
#転置行列
 x=x';
#配列の最初と最後に、Φ(0)=Φ(D-1)=０を加える
phi_value=[0,x,0];
#printf("\n  daubechies D%g phi[integers] values\n\n",D);
#phi_value

q=zeros(1,D);

for i=-D:-1
	if( ((1-i)>=1) &&  ((1-i)<21) )
	   q(i+D+1)=(-1)^(i+1)*p(-i);
	 endif
endfor

#q(1)=p(4);
#q(2)=-p(3);
#q(3)=p(2);
#q(4)=-p(1);

#p
#q

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

#initialize phi array
#phi(n)　n:integers　の値を収納したものは　phi_value(n)
 for i=1:D-1
  phi(reso*i+1)=phi_value(i+1);
 endfor
 
#sub function###############################################
function retval=calc_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=0;
		for m=1:D
                 x=x+p(m)*calc_phi_value(phi,index,my_size);
		   index=index-reso;
		endfor
      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=0;
		for m=1:D
                 x=x+q(m)*calc_phi_value(phi,index,my_size);
		   index=index-reso;
		endfor

      psi(i)=x; 
    endfor
endfor

# plot Daubechies functions & their spectrum

figure(2);
subplot(2,1,1);

t=linspace(0,D-1,my_size);
plot(t,phi,"b-;scaling function;");
s=sprintf("DaubechiesD%g  functions. resolution =1/%d",D,reso);
title(s);

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



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

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

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

plot(x(1:my_size/4),"r-;mother wavelet;");
