analh2d.m

这是伯克里wavelet transforms一书中的例子的代码

function [an,pr,pc] = analh2d(img,al,ah,levs)
%Routine effects a 2-D separable multilevel DWT using biorthogonal
%even length or half sample symmetric (HSS) filters. Symmetric 
%extensions are used to extend the data to perform the DWT.
%
%[an,pr,pc] = analh2d(array,lpf,hpf,L) effects an L level separable DWT
%of the 2-D input array. The analysis filter bank is defined by 
%lpf and hpf. The DWT is returned as the array, 'an'. The arrays 
%'pr' and 'pc' contain the end indices of the rows and columns
%of the subbands.
%
%[an,pr,pc] = analh2d(array,fopt,L) effects an L level separable DWT of the 
%2-D input array. The analysis filter is chosen using the variable 
%fopt from the ones prestored in the routine. The DWT is returned as 
%the array, 'an'. The arrays 'pr' and 'pc' contain the end indices 
%of the rows and columns of the subbands.
%
%After the completion of the routine, the DWT of the array is displayed.
%
%Input variable are:
%array   : 2-D input array
%lpf, hpf: Lowpass and highpass filters constituting the biorthogonal
%          WSS filter bank.
%OR
%fopt    : Choose fopt to be 1, 2 or 3. They correspond to 2/6, 6/2 and 
%          4/4 tap biorthogonal HSS filter banks.
%L       : Number of levels of decomposition
%
%One can separate the subbands using the lists 'pr' and 'pc'that 
%contains the end coordinates of the row and column decomposition.
%
%It is required that the input image at each level of decomposition 
%satisfy the following criteria:
%The number of rows, r. of the input image at that level be such that
%2*r>=lm, where 'lm' is the length of the larger of the filters in the 
%filter bank.
%
%The number of columns, c, of the input image at that level be such that
%2*c>=lm.
%
%The routine checks to see if these conditions 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 image is then 
%decomposed down to this new number of levels.
%
%The DWT is stored in the 2-D array, 'an' and can be displayed using
%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.
%
%Also see the corresponding synthesis routine 'synthh2d'.
% 
%Refer Chapter 4 for information on biorthogonal wavelet 
%decompositions and 2-D separable wavelet transform.
%
%Author: Ajit S. Bopardikar
%Copyright (c) 1998 by Addison Wesley Longman, Inc.
%

%if a filter option is chosen
  num = nargin; %number of input arguments
  if (num ==3)
    levs = ah;
   if(al == 1)
     al =[0.70710678118655   0.70710678118655];
     ah =[-0.17677669529664 -0.17677669529664 0.70710678118655 -0.70710678118655 0.17677669529664 0.17677669529664];
   elseif (al ==2)
     al =[-0.17677669529664 0.17677669529664 0.70710678118655 0.70710678118655 0.17677669529664  -0.17677669529664];
     ah =[-0.70710678118655   0.70710678118655];
   elseif (al >=3)
      if (al > 3)
      fprintf('fopt chosen to be greater than 3. Using fopt=3 instead\n');
      end; 
     al =[0.17677669529664 0.53033008588991 0.53033008588991 0.17677669529664];
     ah =[-0.35355339059327 -1.06066017177982 1.06066017177982 0.35355339059327];
   end %end inner if
  end %enf if
  
  [m,n] = size(img);
  an = img;
  pr =[m];
  pc =[n];

  levr = 0; %initialize
  levc = 0; %initialize
  lr   = m; %initialize
  lc   = n; %initialize
  lo   = length(al); %initialize
  hi   = length(ah); %initialize
  lm   = max(lo,hi); %determine the larger of the filter lengths  

  pr  = [lr]; %initialize the row coordinates
  pc  = [lc]; %initialize the column coordinates

  while (2*lr >=lm)  %for rows
    if (lo/4 == round(lo/4) ) %center odd
      if (lr/2 == round(lr/2)) %even number of rows
        lr = lr/2; 
      else %odd number of rows
        lr = (lr+1)/2;  
      end %end inner if
    else %center even
      if (lr/2 == round(lr/2))  %even number of rows
        lr = lr/2+1; 
      else  %odd number of rows
        lr = (lr+1)/2;  
      end %end inner if
    end %end if
   levr = levr+1; 
   pr = [lr pr];
   end %endwhile

    %columns
   while(2*lc >lm)
    if (lo/4 == round(lo/4) ) %center odd
      if (lc/2 == round(lc/2)) %even number of columns
        lc =lc/2;
      else  %odd number of columns
        lc = (lc+1)/2; 
      end %end inner if
    else  %center even
      if (lc/2 == round(lc/2)) %even number of columns
        lc = lc/2+1; 
      else %odd number of columns
      lc = (lc+1)/2; 
      end %end inner if
    end %end if
   levc = levc+1;  
   pc = [lc pc];
   end %endwhile

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

  lpr = length(pr);
  lpc = length(pc);
  prr = pr(lpr:-1:(lpr-levs)); %equalize the two arrays...
  pcr = pc(lpc:-1:(lpc-levs));
       
  for i = 1:levs
    an1 = hssa1(an(1:prr(i),1:pcr(i)),al,ah); %one level of decomposition
    an(1:prr(i),1:pcr(i)) = an1;  
  end %end for

  pr = prr((levs+1):-1:1); %reverse the lists for returning them...
  pc = pcr((levs+1):-1:1);
  figure;colormap(gray);imagesc(an);title('DWT of Input Array');