analw2d.m

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

function [an,pr,pc] = analw2d(img,al,ah,levs)
%Routine to compute 2-D separable biorthogonal DWT using an odd length or
%whole sample symmetric (WSS) biorthogonal filter bank. Symmetric 
%extensions are used to extend the data to perform the DWT.
%
%an = analw2d(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 result is returned as a 2-D array, 'an'. The arrays
%'pr' and 'pc' contain the end indices of the rows and columns
%of the subbands.
%
%an = analw2d(array, fopt,L) effects an L level separable DWT 
%of the 2-D input array. The analysis filter bank is chosen using fopt
%to be one of the prestored filter banks. The result is returned as a 
%2-D array, 'an'. The arrays 'pr' and 'pc' contain the end indices
%of the rows and columns of the subbands.
%
%Input variables 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 or 2. They correspond to the 9/7 and 7/9
%          tap filters used by FBI for their fingerprint image 
%          compression scheme.
%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-2 >= 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-2 >= 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 the above criteria are met and returns this value. The input 
%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.
%
%The DWT of the image is displayed at the end of the computations anyway.
%
%See also the complementary synthesis routine 'syntw2d'.
%
%Refer to 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.
%

num = nargin; %number of input arguments

  if (num ==3) %filter option specified instead of filters
    levs = ah;
   if(al == 1)
     al =[0.03782845550700 -0.02384946501938 -0.11062440441842 0.37740285561265 0.85269867900940 0.37740285561265 -0.11062440441842 -0.02384946501938 0.03782845550700];
     ah =[0.06453888262894 -0.04068941760956 -0.41809227322221 0.78848561640566 -0.41809227322221 -0.04068941760956 0.06453888262894];
   elseif (al >=2)
      if (al > 2)
      fprintf('fopt chosen to be greater than 2. Using fopt=2 instead\n');
      end; 
     al =[-0.06453888262894 -0.04068941760956 0.41809227322221 0.78848561640566 0.41809227322221 -0.04068941760956 -0.06453888262894];
     ah =[0.03782845550700 0.02384946501938 -0.11062440441840 -0.37740285561265 0.85269867900940 -0.37740285561265 -0.11062440441840 0.02384946501938 0.03782845550700];
   end %end inner if
  end %enf if

%Algorithm to check if the number of levels input is valid.

  [m,n]=size(img);
  
  levr=0;          %initialize
  lr = m;          %initialize
  levc=0;          %initialize
  lc = n;          %initialize
  lo = length(al); %initialize
  hi = length(ah); %initialize  
  lm = max(lo,hi); %determine the larger of the filter lengths
  pr = [m];        %initialize the row coordinate vector
  pc = [n];        %initialize the column coordinate vector 

  while ((2*lr-2) >= lm)
    if (lr/2 == round(lr/2)) %even length sequence
      lr=lr/2;
    else
      if (lo>hi) %if low pass filter has larger length...
        lr=ceil(lr/2);
      else       %if low pass filter has smaller length...
        lr=floor(lr/2);
      end %endif
    end %endif
   pr = [pr lr];
   levr=levr+1;
  end %endwhile

  while ((2*lc-2) >= lm)
    if (lc/2 == round(lc/2)) %even length sequence
      lc=lc/2;
    else
      if (lo>hi) %if low pass filter has larger length...
        lc=ceil(lc/2);
      else       %if low pass filter has smaller length...
        lc=floor(lc/2);
      end %endif
    end %endif
   pc = [pc lc];
   levc=levc+1;
  end %endwhile
   
  lev1 = min(levr,levc); %find the smaller of the calculated levels

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

  pc = pc(1:(levs+1));
  pr = pr(1:(levs+1)); 

   an =img;

  for i=1:levs
    an1 = wssa1(an(1:pr(i),1:pc(i)),al,ah);
    an(1:pr(i),1:pc(i))=an1; 
  end %end for 
  
  pr = pr((levs+1):-1:1);
  pc = pc((levs+1):-1:1);

%display the DWT of the input 2-D array
  figure;colormap(gray);imagesc(an);title('DWT of Input Array');