www.pudn.com > Upload_EZW128.zip > ANAL2D.M

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