function yy = Psijk(z, j, k, filt, n)
%--------------------------------------------------------------
%  yy=Psijk(z, j, k, filter, n)
%  Evaluation of the wavelet function corresponding to an Orthogonal 
%           MRA by Daubechies-Lagarias Algorithm.
%  inputs:  z -- the argument
%           j -- scale
%           k -- shift
%           filter -- ON finite wavelet filter, might be an
%                     output of WaveLab's: MakeONFilter
%           n -- precision of approximation maesured by the number
%                of Daubechies-Lagarias steps (default n=20)
%--------------------------------------------------------
%  output:  yy -- value of mother wavelet (j,k) coresponding to
%                  'filter' at z.
%--------------------------------------------------------------
% Example of use: 
% > xx = -1:0.01:5;  yy=[]; 
% > for i=1:length(xx)
% > yy =[yy Psijk(x(i), 0, 1, MakeONFilter('Daubechies',4), 25)];
% > end
% > plot(x,yy)
%---------------------------------------------------------------- 
    if (nargin == 4) n=20;
    end
    N=length(filt)-1;
    daun = (N+1)/2;
    x=(2^j)*z-k;
if(x<=daun-N|x>=daun) yy=0;
     %if(x<=0|x>=N) yy=0;
else
    twox = 2 * x;
    inti=floor(twox);
    dec=twox-inti;
    dy=dec2bin(dec,n);
    t0=t0(filt);
    t1=t1(filt);
    prod=eye(N);
    for i=1:n
            if dy(i)==1
                prod=prod*t1;
            else prod=prod*t0;
            end
        end
    uu=[];
    for i=1:N
        index = i + 1 - inti;
        if ( index > 0 & index < N + 2 ) 
             fi= (-1)^(-index)*filt(index);
        else fi=0;
        end   
        uu =[uu   fi];
    end
   %----------------------------------------
   v = 1/N * ones(1,N) * prod' ;
    yy=2^(j/2)* uu * v';
   end
        
  %--------------------functions needed----------------------------
  %---------
  function a = dec2bin(t,n)
            a=[];
             for i = 1:n
             if(t <= 0.5) a=[a 0]; t=2*t;
              else a=[a 1]; t=2*t-1;
             end
            end
 %-----------
 function t0 = t0(filt)
%
n = length(filt);
nn = n - 1;
%
t0 = zeros(nn);
for i = 1:nn
    for j= 1:nn
        if (2*i - j > 0 & 2*i - j <= n) 
            t0(i,j) = sqrt(2) * filt( 2*i - j ); 
        end
    end
end
%------------------
function t1 = t1(filt)
%
n = length(filt);
nn = n - 1;
%
t1 = zeros(nn);
for i = 1:nn
    for j= 1:nn
        if (2*i -j+1 > 0 & 2*i - j+1 <= n) 
            t1(i,j) = sqrt(2) * filt( 2*i - j+1 ); 
        end
    end
end
%---------------------  B. Vidakovic, 2002    --------------------
%-----------------------------------------------------------------