function an = anal2d(img,fopt,levs)
%Routine for computing a separable, multilevel 2-D DWT 
%using orthonormal filters. Circular extension used to
%extend the input 2-D data.
%
%an = anal2d(image,lpf,L) computes a separable L-level 2-D DWT of
%the input image using the filter bank generated by lpf, the 
%low pass filter of an orthonormal filter bank. The result is 
%returned as a 2-D array, 'an'. 
%
%an = anal2d(image,fopt,L) computes a separable L-level 2-D DWT 
%of the 2-D image array using a filter bank generated by a lowpass 
%filter chosen from the ones provided in the routine, using fopt. 
%The result is returned as a 2-D array, 'an'. 
%
%The one-level separable DWT of an M by N images consists of 4 subbands. 
%Each of these subbands are M/2 by N/2 in size. Thus the matrix an 
%consists of the four subbands arraged as follows:
%(Note: Locations are represented as (row,column)). 
%
%The LL subband occupies the locations from (0,0) to (M/2,N/2).  
%The LH subband occupies the locations from (0,M/2+1) to (M,N/2).
%The HL subband occupies the locations from (N/2+1,0) to (M/2,N).
%The HH subband occupies the locations from (M/2,N/2) to (M,N).
%
%For further levels of decomposition the LL subband is decomposed into
%four subbands that are arranged in a similar manner. 
%
%It is required that the input image at each level satisfy the following
%criteria:
%1. The number of rows and columns be even.
%2. The number of rows and columns be greater than or equal to 
%   the filter length.
%The routine checks to see if the above criteria are met at each level.
%If not, then the routine calculates the maximum number of levels 
%for which these criteria are met and returns this value.
%The input image is decomposed down to this new number of levels. 
%
%
%Input variables are:
%array : 2-D array
%lpf   : Lowpass filter of the orthonormal filter bank
%OR
%fopt  : which can take a value of 1 to 5 and corresponds to 
%        Daubechies 2*fopt tap filter respectively.
%L     : Specifies the number of levels of decomposition.
%
%The DWT stored in array 'an' can be displayed using the 
%'image' or the 'imagesc' command. Use the 'colormap(gray)' command 
%to specify the colormap used. Please refer to online help on these
%commands for more information.
% 
%The DWT of the image is displayed at the end of the computations
%anyway. 
%
%See also the complementary synthesis function 'synth2d'.
%
%Refer to Chapter 4 for more information on 2-D separable DWT using 
%orthonormal filters.
%
%Author: Ajit S. Bopardikar
%Copyright (c) 1998 by Addison Wesley Longman, Inc.
%
img=img-128;
  if(prod(size(fopt))==1) %you want to use one of the filter options...
    if (fopt==1) %Daubechies 2 or Haar case
      lpf = [1/sqrt(2) 1/sqrt(2)];
    elseif (fopt == 2) %Daubechies 4
      lpf =[0.48296291314453 0.83651630373781 0.22414386804201 -0.12940952255126];
    elseif (fopt == 3) %Daubechies 6
      lpf =[0.33267055295000 0.80689150931100 0.45877502118000 -0.13501102001000 -0.08544127388200 0.03522629188200];
    elseif (fopt == 4) %Daubechies 8
      lpf =[0.23037781330900 0.71484657055300 0.63088076793000 -0.02798376941700 -0.18703481171900 0.03084138183600 0.03288301166700 -0.01059740178500];
    elseif (fopt>= 5) %Daubechies 10
      if (fopt > 5)
      fprintf('fopt chosen to be greater than 5. Using fopt=5 instead\n');
      end; 
      lpf =[0.16010239797400 0.60382926979700 0.72430852843800 0.13842814590100 -0.24229488706600 -0.03224486958500 0.07757149384000 -0.00624149021300 -0.01258075199900 0.00333572528500];
    end %end inner if
  else %input filter
    lpf = fopt;
  end %end if
  lf  = length(lpf); 

   for i=0:(lf-1)
    hpf(i+1) = (-1)^i*lpf(lf-i);
   end; %endfor
%so we have the high pass filter here.

   [m,n] = size(img);

  levr =0;   %initilaize 
  lr = m;  %initialize
  levc =0; %initilaize
  lc = n;  %initilaize

 while ((lr/2==round(lr/2)) & (lr>= lf))
   lr = lr/2;
   levr=levr+1;
 end
 
 while ((lc/2==round(lc/2)) & (lc>= lf))
   lc = lc/2;
   levc=levc+1;
 end

 lev1 = min(levr,levc);
 if (lev1<levs)
   fprintf('Cant decompose to %d levels. Decomposing to %d levels instead\n',levs,lev1);
   levs = lev1;
 end %endif

 an = img; %initialize
 for i=1:levs
    an1 = anal21(an(1:m,1:n),lpf,hpf); 
    %one level of decomposition
    
    an(1:m,1:n) = an1; 
    m = m/2;
    n = n/2;
 end; %after all the levels of wavelet decomposition

 %figure;colormap(gray);imagesc(an);title('DWT of Input Matrix');
 imshow(mat2gray(an));title('DWT of Input Matrix');