anal.m

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

function an = anal(x,fopt,lev)
%Routine for multilevel DWT decomposition using an orthonormal filter bank.
%
%an = anal(data,lpf,L) computes an L-level DWT of the input vector, 
%data, using the filter bank generated by the lowpass orthonormal 
%filter lpf. The result is returned as an array, an. Circular extension 
%is used.
% 
%an = anal(data,fltopt,L) computes an L-level DWT of the input 
%vector, data, using the filter bank generated by a lowpass filter 
%chosen from the ones provided in the routine, using fltopt. The 
%result is returned as an array, an. Circular extension is used.
%
%One level of decomposition of a length L sequence results in 
%a lowpass sequence of length L/2 and a highpass (detail) sequence 
%of length L/2. 
%
%The DWT coefficients computed by the routine are arranged in array an 
%such that the lowpass coefficents (Lth level lowpass coefficients) 
%occupy the first locations, followed by the Lth level detail 
%coefficients and so on till the 1st level detail coefficients which 
%occupy the last locations in the array. 
%    For example, if we have a 256 sample sequence decomposed down to 
%3 levels then the three levels are arranged in array an as follows:
%Locations 1-32 in an contain the lowpass coefficients. Locations 32-64
%contain the level 3 detail coefficients. Locations 65-128 contain the 
%second level detail coefficients and locations 129-256 contain the 
%level 1 detail coefficients.
%    
%It is required that the input sequences at each level satisfy the 
%following criteria:
%1. It always have EVEN number of samples.
%2. The number of samples be greater than 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 sequence is decomposed down to this new number of levels.
%
%The subband coefficients of the decompositon are plotted starting 
%with the level 1 detail subband at top left of the subplot and 
%ending with with the final level lowpass subband at bottom right.
%
%The input variables are:
%data: The input sequence.
%lpf : Lowpass filter of the orthonormal filter bank
%OR
%fltopt : which can take a value of 1 to 5 and corresponds to 
%         Daubechies 2*fltopt tap filters respectively.
%levs   : Specifies the number of levels of decomposition.
%
%See also the complementary DWT synthesis function 'synth'.
%
%Refer to Chapter 3 for more information on DWT using orthonormal
%filters.

  if(prod(size(fopt))==1) %you want to use one of the filter options...
    if (fopt==1) %Daubechies 2 or Haar case
       f = [1/sqrt(2) 1/sqrt(2)];
    elseif (fopt == 2) %Daubechies 4
      f =[0.48296291314453 0.83651630373781 0.22414386804201 -0.12940952255126];
    elseif (fopt == 3) %Daubechies 6
      f =[0.33267055295000 0.80689150931100 0.45877502118000 -0.13501102001000 -0.08544127388200 0.03522629188200];
    elseif (fopt == 4) %Daubechies 8
      f =[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; 
      f =[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
    f = fopt;
  end %end if

 %generate the high pass filter for the given low pass filter
 lf = length(f);
  
 for i=0:(lf-1)
   h(i+1) = (-1)^i*f(lf-i);
 end; %endfor
%here we have the high pass filter here.

 ll = length(x); %initialize
 an = x;        %initialize

%verify if you can really go down to the specified number of levels...
 lev1 =0;   %initilaize 
 ll1 = ll;  %initialize

 while ((ll1/2==round(ll1/2)) & (ll1>= lf))
   ll1 = ll1/2;
   lev1=lev1+1;
 end

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

%initialize the number of rows in the subplot
 if(lev/2 == round(lev/2))
   sr = lev/2 +1;
 else
   sr = round(lev/2);
 end %end if

 for i=1:lev
    an1 = anal1(an(1:ll),f,h); 
    %one level of decomposition
    
    an(1:ll) = an1(1:ll); 
    ll = ll/2;
    subplot(sr,2,i);
    t = ['Level' ' ' num2str(i) ' ' 'detail coefficients']; 
    plot(an(ll+1:2*ll));title(t);
    %plot high pass subbands

 end; %after all the levels of wavelet decomposition, plot the lowpass
%subband
   subplot(sr,2,lev+1);
   plot(an(1:ll));title('Lowpass Coefficients');