analh.m

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

function [an,pl] = analh(data,al,ah,levs)
%Routine to compute  DWT using an even length half sample symmetric 
%(HSS) biorthogonal filter bank. Symmetric extensions are used to 
%extend the data to perform the DWT.
%
%[an,p] = analh(array,lpf,hpf,L) effects an L level DWT of the input array. 
%The analysis filter bank is defined by lpf and hpf. The decomposed 
%sequence is returned in 'an', while 'p' contains the end index of 
%each subband.
%
%[an,p] = analh(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 decomposed sequence is returned in 'an', 
%while 'p' contains the end index of each subband.
%
%Individual subbands are plotted after the decomposition starting with
%the low pass band and ending with the detail function at the first level.
%
%Input variables are
%array   : Input array
%lpf, hpf: Lowpass and highpass filters constituting the biorthogonal
%          HSS filter bank.
%OR
%fopt    : Choose fopt to be 1,2 or 3. They correspond to the 2/6 and 6/2
%           4/4 tap filter banks respectively.
%L       : Number of levels of decomposition.
%
%One can separate the subbands using the list 'p' that contains the end
%coordinates of each subband, starting with the lowpass 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 >=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 data is
%then decomposed down to this new number of levels.  
%
%See also the complementary synthesis routine 'synthh'.
%
%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 the number of input arguments is 3 then use the second argument to
%pick up the filter bank and the third argument gives the number of levels
  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
  
  l =length(data);
  an = data;
  pl = [l]; %initialize an array for plotting data....

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

  lev1 = 0; %initialize
  l1   = l; %initialize
  lo   = length(al); %initialize
  hi   = length(ah); %initialize
  lm   = max(lo,hi); %determine the larger of the filter lengths
  
  while(2*l1>= lm) 
    if (lo/4 == round(lo/4)) %center odd
      if(l1/2 == round(l1/2)) %length of data even
        l1=l1/2; 
      else %length odd
        l1=(l1+1)/2; 
      end %end inner if
    else   %center even
      if(l1/2 == round(l1/2))
        l1 = l1/2 +1; 
      else
        l1 = (l1+1)/2; 
      end %end inner if
    end %end outer if
  lev1 =lev1+1; 
  pl = [l1 pl]; %update the plot array
  end %end while
 
 if (lev1<levs)
   fprintf('Cant decompose to %d levels. Decomposing to %d levels instead.\n',levs,lev1);
   levs=lev1;
 end  %endif

 lp = length(pl);
 plr= pl(lp:-1:lp-levs);
  for i=1:levs    
     an1 = bianal2(an(1:plr(i)),al,ah); %one level decomposition
     an(1:plr(i)) = an1;
 end %endfor 

 if(levs/2 == round(levs/2)) %determine the number of rows in subplot
   sr = levs/2 +1;
 else
   sr = round(levs/2);
 end
  
  pl = plr((levs+1):-1:1);
  subplot(sr,2,1)
  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