wsss1.m

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

function sy = wsss1(an,sl,sh)
%sy = wsss1(an,lpf,hpf) performs a one level reconstruction of a 
%1 level DWT, an, using whole sample symmetric (WSS) synthesis filters 
%lpf and hpf. The resultant synthesized 2-D function is returned to sy.
%
%This routine is used by the routine syntw2d.
%
%Author: Ajit S. Bopardikar
%Copyright (c) 1998 by Addison Wesley Longman, Inc.
%

  [m,n] = size(an);

  ll = length(sl);
  lh = length(sh);

  le = (ll-1)/2;
  he = (lh-1)/2;

%first reconstruct the columns

  if (m/2 == round(m/2) ) %if even number of rows
    low = an(1:m/2,:)'; %isloate the low pass columns
    hi  = an((m/2+1):m,:)'; %isolate the high pass columns 

    %upsample the subbands
    lu = zeros([n,m]);
    lu(:,1:2:m) = low;
    hu = zeros([n,m]);
    hu(:,1:2:m) = hi;

   if ll<lh %for analysis filters ll>lh
     extl = [lu(:,le+1:-1:2) lu lu(:,m-1:-1:m-le-1)]; %LPF extension
     exth = [hu(:,he:-1:1) zeros([n 1]) hu hu(:,m-3:-1:m-he-4)]; 
     %HPF extension

   else %for analysis filters ll<lh
     extl = [lu(:,le:-1:1) zeros([n 1]) lu lu(:,m-3:-1:m-le-4)]; 
     %LPF extension
     exth = [hu(:,he+1:-1:2) hu hu(:,m-1:-1:m-he-1)]; %HPF extension
   end %end inner if for extensions

   r0 = conv2(extl,sl); %reconstruct the low pass part
   r1 = conv2(exth,sh); %reconstruct the high pass part

   syc = r0(:,ll:1:ll+m-1)'+r1(:,lh:1:lh+m-1)';

 else %if the number of rows is odd...
   
   if (ll>lh)
     %isolate the subbands 
     lo = an(1:floor(m/2),:)';
     hi = an((floor(m/2)+1):m,:)';
     %upsample
     lu = zeros([n m]);
     hu = zeros([n m]);
     lu(:,2:2:m) = lo;
     hu(:,1:2:m) = hi;

   else
     %isolate the subbands
     lo = an(1:ceil(m/2),:)';
     hi = an(ceil(m/2)+1:m,:)';
     %upsample
     lu = zeros([n m]);
     hu = zeros([n m]);
     lu(:,1:2:m) = lo;
     hu(:,2:2:m) = hi;
   end %end inner if 

  extl = [lu(:,le+1:-1:2) lu lu(:,m-1:-1:m-le-1)];
  exth = [hu(:,he+1:-1:2) hu hu(:,m-1:-1:m-he-1)];

  r0 = conv2(extl,sl);
  r1 = conv2(exth,sh);

  syc = r0(:,ll:m+ll-1)'+ r1(:,lh:m+lh-1)';

 end %end if 

%this finishes the columns now for the rows

  if (n/2 == round(n/2) ) %if even number of columns
    low = syc(:,1:n/2); %isloate the low pass columns
    hi  = syc(:,(n/2+1):n); %isolate the high pass columns 

    %upsample the subbands
    lu = zeros([m,n]);
    lu(:,1:2:n) = low;
    
    hu = zeros([m,n]);
    hu(:,1:2:n) = hi;

   if ll<lh %for analysis filters ll>lh
     extl = [lu(:,le+1:-1:2) lu lu(:,n-1:-1:n-le-1)]; %LPF extension
     exth = [hu(:,he:-1:1) zeros([m 1]) hu hu(:,n-3:-1:n-he-4)];
     %HPF extension

   else %for analysis filters ll<lh
     extl = [lu(:,le:-1:1) zeros([m 1]) lu lu(:,n-3:-1:n-le-4)];
     %LPF extension
     exth = [hu(:,he+1:-1:2) hu hu(:,n-1:-1:n-he-1)]; %HPF extension
   end %end inner if for extensions

   r0 = conv2(extl,sl); %reconstruct the low pass part
   r1 = conv2(exth,sh); %reconstruct the high pass part

   sy = r0(:,ll:1:ll+n-1)+r1(:,lh:1:lh+n-1);

   else %if the number of columns is odd...
   
   if (ll>lh)
     %isolate the subbands 
     lo = syc(:,1:floor(n/2));
     hi = syc(:,(floor(n/2)+1):n);
     %upsample
     lu = zeros([m n]);
     hu = zeros([m n]);
     lu(:,2:2:n) = lo;
     hu(:,1:2:n) = hi;

   else
     %isolate the subbands
     lo = syc(:,1:ceil(n/2));
     hi = syc(:,ceil(n/2)+1:n);
     %upsample
     lu = zeros([m n]);
     hu = zeros([m n]);
     lu(:,1:2:n) = lo;
     hu(:,2:2:n) = hi;
   end %end inner if 

  extl = [lu(:,le+1:-1:2) lu lu(:,n-1:-1:n-le-1)];
  exth = [hu(:,he+1:-1:2) hu hu(:,n-1:-1:n-he-1)];

  r0 = conv2(extl,sl);
  r1 = conv2(exth,sh);

  sy = r0(:,ll:n+ll-1)+ r1(:,lh:n+lh-1);

 end %end if