wpsynth.m

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

function sy = wpsynth(an,fopt,tree)
%Routine for 1-D wavelet packet reconstruction.
%This routine is the complement to the wavelet packet analysis routine
%'wpanal'.
%
%sy = wpsynth(an, lpf, tree) takes the wavelet packet decomposed signal,
%an, and reconstructs the sequence using the filter bank defined by the 
%lowpass filter, lpf. The reconstruction follows the path specified by 
%by the tree structure given in the array, tree.
%
%sy = wpsynth(an, fopt, tree) takes the wavelet packet decomposed signal,
%an, and reconstructs the sequence using the filter bank specified by the
%option, fopt. The filter bank is chosen from those  provided in the 
%routine. The reconstruction follows the path specified by the tree 
%structure given in the array, tree.
%
%If the filter bank and the tree structure correspond to the filter bank
%and the  tree structure used for wavelet packet decomposition then 
%we get perfect reconstruction of the original sequence.
%
%The input variables are:
%an            : wavelet packet decomposed 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 respectively. 
%tree          : Specifies the tree for wavelet packet synthesis.
%
%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
      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 low pass and the high pass synthesis filter
  lf = length(f);
 
  lo = f(lf:-1:1);

   for i=0:(lf-1)
     h(i+1) = (-1)^i*f(lf-i);
   end; %end for
   hi = h(lf:-1:1);
%so we have the high pass filter here.

%Now determine the tree and the number of levels...this is necessary 
%because the tree that has been inputted may be overspecified or may 
%require pruning.
  sy = [];  %remove this later
  ls = length(an);  %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

 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       

%Now begin the reconstruction.
 
   lt1 = length(tree1); %Compute the length of the    

   for lvs = levt-1:-1:1 
      l=1;         %initialize the level counter.
      stl = [0];   %initialize the starting point array.
      ls1 = ls;    %initialize the length.
      k=2;         %initialize counter.
      count =2;    %initialize counter.
  
      %start while loop...this will give you the start points of the subbands
        while (l <= lvs)
          count1 = 0; 
          stlt = []; %initialize temporary starting point array
          n = 1;     %initialize the stlt counter  

          for (i= k:(k+count-1))
            count1 = count1+2*tree1(i); %update count1     
            
            if(tree(i)==0) %No decomposition

              n=n+1; %Simply increment 

            else           %Decomposition

              stlt = [stlt (stl(floor((n+1)/2))+ rem(n+1,2)*ls1/2)];
              n=n+1; %update after appending the starting point to stlt
               
            end %endif
     
          end  %endfor

           stl =stlt;
           k = k+count;
           count = count1;
           ls1 = ls1/2;
           l = l+1;

        end %endwhile            
%here you have the starting coordinates and length of the subbands...
%begin 1 level synthesis

     len = length(stl);

     for(i=1:len)
       %one level synthesis...
       an(stl(i)+1:stl(i)+ls1) = syn1(an(stl(i)+1:stl(i)+ls1),lo,hi);
     end %endfor

 end %endfor
%first level reconstruction
    sy = syn1(an,lo,hi);

 else %if the length of the signal is odd or a tree with no initial 
     %decomposition is specified...
   sy =an;  %simply return the input array.
 end %endif

plot(sy);title('Reconstructed sequence');