analw.m

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

function [an,pl] = analw(dat,al,ah,levs)
%Routine to compute  DWT using odd length or whole sample symmetric 
%(WSS) biorthogonal filter bank. Symmetric extensions are used to 
%extend the data to perform the DWT.
%
%[an,p] = analw(array,lpf,hpf,L) effects an L level DWT of the input array. 
%The analysis filter bank is defined by lpf and hpf. The output is 
%returned in array 'an', while 'p' contains the end index of each 
%sequence.
%
%[an,p] = analw(array, fopt,L) effects an L level  DWT of the input array. 
%The analysis filter bank is chosen using fopt to be one of the 
%prestored filter banks. The output is returned in array 'an', while 'p'
%contains the end index of each sequence.
%
%
%Input variables are
%array   : 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 list 'p' that contains the end
%index of each subband, starting with the low pass subband.
%
%Individual subband coefficients are plotted after the decomposition 
%starting with the low pass subband positioned in the top left corner 
%of the subplot and ending with the detail function at the first level 
%positioned at the bottom right of the subplot.
%
%It is required that at each level of decomposition, the input to the
%filter bank at that level have length 'len' such that '2*len-2 >=lm',
%where, 'lm' is the length of the larger of the filters in the filter 
%bank. 
%
%The routine checks to see if this condition is met at each level. If 
%not, then the routine calculates the maximum number of levels for 
%which this criteria is met and returns this value. The input array is
%then decomposed down to this new number of levels.  
%
%See also the complementary synthesis function 'synthw'.
%
%Refer to Chapter 4 for information on biorthogonal wavelet 
%decompositions.
%
%Author: Ajit S. Bopardikar
%Copyright (c) 1998 by Addison Wesley Longman, Inc.
%

num = nargin; %number of input arguments

  if (num ==3)
    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

ll = length(dat); %initialize
an = dat;         %initialize
pl = [ll];        %initialize the subband length array

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

  lev1=0;          %initialize
  ll1 = ll;        %initialize
  lo = length(al); %initialize
  hi = length(ah); %initialize  
  lm = max(lo,hi); %determine the larger of the filter lengths
  
  while ((2*ll1-2) > lm)
    if (ll1/2 == round(ll1/2)) %even length sequence
      ll1=ll1/2;
    else
      if (lo>hi) %if low pass filter has larger length...
        ll1=ceil(ll1/2);
      else       %if low pass filter has smaller length...
        ll1=floor(ll1/2);
      end %endif
    end %endif
   lev1=lev1+1;
   pl=[pl ll1];
  end %endwhile

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

%begin DWT
 for i=1:levs
   if (pl(i)/2 == round(pl(i)/2))       %even length
    an1 = bianal1(an(1:pl(i)),al,ah); %1 level DWT for the even length case
   else                           %odd length 
    an1 = bianal3(an(1:pl(i)),al,ah); %1 level DWT for the odd length case
   end                            %endif
    %one level of decomposition
    
    an(1:pl(i)) = an1(1:pl(i)); %update the DWT array
 end; %after all the levels of wavelet decomposition

 if (levs/2 == round(levs/2)) %even number of levels
   sr = levs/2 + 1;
 else
   sr = round(levs/2); %number of rows in the subplot
 end
 
 pl =pl((levs+1):-1:1);

 subplot(sr,2,1); %top left corner
 plot(an(1:pl(1)));title('Coarse Approximation Coefficients');
%plot the low pass part....

%plot the detail functions

 for i= 1:levs 
   subplot(sr,2,i+1);
   t = ['Level' ' ' num2str(levs-i+1) ' ' 'Detail Coefficients'];
   plot(an(pl(i)+1:pl(i+1)));title(t);
 end %end for