wpanal.m

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

function [an,tree1] = wpanal(seq,fopt,tree)
%Routine for 1-D wavelet packet decomposition.
%
%[an,tree1]= wpanal(seq, lpf, tree) computes a wavelet packet decomposition
%of the input sequence, seq. The wavelet packet decomposition is 
%specified by the array, tree, and the orthogonal filter bank used for 
%decomposition is defined by the lowpass filter, lpf. The wavelet packet
%decomposed sequence is returned in the array, an, while the decomposition
%tree is returned in the array, tree1.
%
%[an,tree1]= wpanal(seq, fopt, tree) computes a wavelet packet decomposition
%of the input sequence, seq. The wavelet packet decomposition is specified
%by the array, tree, and the orthogonal filter bank is chosen using the 
%the variable, fopt, from the filter banks provided in the routines. The
%wavelet packet decomposed sequence is returned in the array, an, while 
%the decomposition tree is returned in the array, tree1.
%
%The subbands are displayed starting with the lowpass subband (or 
%subband 0) at the top left.
%
%The tree array specifies the wavelet packet decomposition as follows:
%a. A 1 represents that a subband is to be decomposed and 0 represents 
%   that the subband is not to be decomposed. 
%b. Subbands that are decomposed at the present level are the only ones 
%   that can be decomposed at the next level. The tree array should
%   contain entries only for these subbands. 
%c. The first element of the tree array should be 1. This represents
%   the fact that the input array, seq, is decomposed into two subbands
%   to begin with.
%
% An example of a tree array is given below. Here, the lower subband 
%is lowpass  and the upper subband is highpass.
%
%                       __
%                      |      __
%                    __|   __|
%                   |  |  |  |__
%                   |  |__|
%                  -|     |   __
%                   |     |__|
%                   |__      |__
%                    
%tree = [1,0,1,1,0,1,1]
%
%The sequence, seq, is first split into two subbands. This is represented 
%by the first entry of the tree array being 1. The next two entries are
%for the resultant subbands. A 0 indicates that the lowpass subband is 
%not to be decomposed and the 1 indicates that the highpass subband is 
%further decomposed into two subbands. For the next level only these two
%subband will be considered. The next two entries of the array specifies 
%that the lowpass subband needs to be decomposed into two subbands while
%the highpass subband is not decomposed. Thus, for the next level, only 
%these two entries are considered. The last two entries of the tree
%array are both 1s. Hence both the subbands are decomposed.  
%
%It is required that the input sequences at each level satisfy the follwing
%criteria:
%1. The sequences should have an EVEN number of samples.
%2. The number of samples in the input sequences should greater than
%   the filter length.
%
%The routine check to see if these criteria are satisfied. If not then
%the routine computes the wavelet packet decomposition down to the 
%level that they are met and returns the resultant decompostion along 
%with the pruned input tree (in the array, tree1).
%
%The input variables are:
%seq           : The input sequence
%lpf           : Lowpass filter of the orthonormal filter bank.
%OR
%fopt          : which can take a value of 1 to 5 and corresponds to 
%                Daubechies 2*fopt tap filters. 
%tree          : Specifies the tree for wavelet packet decomposition.
%
%See also the complementary wavelet packet synthesis routine 'wpsynth'
%
%Refer to Chapter 4 for more information on wavelet packets. 
%
%Author: Ajit S. Bopardikar
%Copyright (c) 1998 by Addison Wesley Longman, Inc.
%

  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
      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.

  ls = length(seq);  %length of the sequence
  lt = length(tree); %length of the tree
  ls1 =ls/2;         %initialize the length counter
  levt =0;           %initialize the level counter
%Determine the number of levels using the tree structure
  an = [];           %initialze the wavelet packet array

 if ((tree(1) == 1) & (ls1 == round(ls1))) %there is an initial decomposition
   levt =1;       %initialize the levels
   count = 2;     %initialize the count variable
   k = 2;         %initialze the other counter
   tree1 = [1];   %initialize the tree to be output
   
%All the following four conditions need to be satisfied...
   while((k<lt) & (count~=0) & (ls1/2==round(ls1/2)) & (ls1>=lf))
     count1 = 0;
     for(i = k:(k+count-1))
       count1 = count1 +2*tree(i);
       tree1 = [tree1 tree(i)];  %Append from the previous tree
     end %endfor

     k = k+count; %update variables
     count = count1; %update count variable
     levt = levt + 1; %increment the levels
     ls1 =ls1/2; %reduce length by half for next level...
   end %end while       
 
%begin the decomposition...

%First only the first level... INITIALIZATION HERE
   pls = [];       %initialize the plot list.
   stl = [0,ls/2]; %Starting point of the subbands...
   ls1 = ls/2;    %initiate the length counter...
   an1 = anal1(seq,f,h); %perform a 1-level initial decomposition
   dc1 = [an1(1:ls1);an1(ls1+1:2*ls1)]; %create the initial decompostion
                                        %list
   count = 2;   %Initialize the initial tree list position variable.
   k = 2;       %Initialize the tree list position counter.
   lt1 = length(tree1); %length of the new tree array.
   an = []; %initialize the analysis array
%here the initialization ends 


   while(k<lt1)  %for the entire tree 
     count1 = 0; %a temporary count variable to be intialized...
    
     j=1; %this the counter for the stl and the dc1 arrays...
     dc1temp =[]; %initialize the temporary decomposition list

     for(i=k:(k+count-1))
       count1 = count1+2*tree1(i); %update the temporary variable...

       if(tree1(i)==0) %then append this to the analysis list
      
         [m,n]=size(dc1);                    %determine the dimensions of dc1
         an(stl(2*j-1)+1:stl(2*j-1)+n) = dc1(j,:); %append to the analysis list
         dc1=[dc1(1:j-1,:);dc1(j+1:m,:)];    %drop the appended subband
         pls = [pls,stl(2*j-1)];                  %Append to the plot list
         stl = drop(stl,2*j-1);              %Remove that entry from the list

       else   %if tree(i)==1
         
         an1 = anal1(dc1(j,:),f,h);   %decompose and separate the subbands
         dc1temp = [dc1temp; an1(1:ls1/2); an1(ls1/2 +1:ls1)];
         stl = insert(stl,stl(2*j-1)+ls1/2,2*j); %Insert the starting point
         j=j+1; %increment j    
    
       end  %endif
      
     end %endfor

    k=k+count; %update k to count for the next level
    dc1 = dc1temp; %update the subband array
    count = count1;
    ls1 = ls1/2;
   end %endwhile  - This finishes the decompostion part. 

  %Now put in other subbands
   [m,n] = size(dc1);
   for(i=1:m)
     an(stl(i)+1:stl(i)+n)=dc1(i,:);
   end %endfor 
    pls = [stl pls ls];
    pls = sort(pls);

%Now display the suubands
   lp = length(pls)-1;

   for(i=1:lp)
     subplot(round(lp/2),2,i);
     plot(an(pls(i)+1:pls(i+1)));
     t =['Subband' ' ' num2str(i-1)];
     title(t);
   end %endfor
   
     
  else  %No decomposition...
    an = seq;    %return the input sequence as no decompostion is possible
    tree1 = [0]; %tree is just this-no decompostion...
    plot(an); title('No decomposition possible.');
    fprintf('No decomposition possible. returning the original sequence');
  end %endif