function sy = synth2d(an,fopt,levs)
%Routine to reconstruct the 2-D image from its orthonormal DWT. Circular
%extension used to extend the subbands.
%
%This routine is complement to the analysis routine 'anal2d'.
%
%sy = synth2d(an,lpf,L) takes the L level DWT of the 2-D image 
%and reconstructs it using the orthonormal synthesis filter bank
%defined by the filter lpf. The reconstructed image is stored in 
%the matrix 'sy'. If lpf corresponds to the analysis filter used 
%for the DWT, then perfect reconstruction will be obtained.
%
%sy = synth2d(an,fopt,L) takes the L level DWT of the 2-D image
%and reconstructs it using the synthesis filter bank defined by a 
%filter stored in the routine. This is chosen using the variable fopt.
%The reconstructed image is stored in the matrix 'sy'. If the 
%synthesis filter bank thus chosen corresponds to the analysis 
%filter bank used then perfect recontruction will be obtained.
%
%The input variables are:
%an   : 2-D array containing the DWT.
%lpf  : Lowpass filter used to construct the synthesis filter bank.
%OR 
%fopt : which can take a value from 1 to 5 and corresponds to Daubechies
%       2*fopt tap filters respectively
%L    : The number of levels down to which the DWT has been performed.
%
%The reconstructed image stored in 'sy' can be displayed using the 
%image or the imagesc command along with the colormap command. Please
%refer to MATLAB help for more information on these commands.
%
%The reconstructed image is displayed at the end of the routine.
%
%Refer to Chapter 4 for more information on 2-D separable DWT. 
%
%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
       lpf = [1/sqrt(2) 1/sqrt(2)];
    elseif (fopt == 2) %Daubechies 4
      lpf =[0.48296291314453 0.83651630373781 0.22414386804201 -0.12940952255126];
    elseif (fopt == 3) %Daubechies 6
      lpf =[0.33267055295000 0.80689150931100 0.45877502118000 -0.13501102001000 -0.08544127388200 0.03522629188200];
    elseif (fopt == 4) %Daubechies 8
      lpf =[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; 
      lpf =[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
    lpf = fopt;
  end %end if
  lf = length(lpf);
 
  lo = lpf(lf:-1:1);

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

   [m,n] = size(an); %dimensions of the original image

   m = m/(2^levs); %the dimensions of 
   n = n/(2^levs); %lowest LL subband

   an1 = an;

   for i=1:levs
     sy = synth21(an1(1:2*m,1:2*n),lo,hi);
     an1(1:2*m,1:2*n) = sy;
     m = 2*m;
     n = 2*n;
   end; %endfor
   sy=sy+128;
   %figure;colormap(gray);imagesc(sy);title('IDWT of input DWT array'); 
   %imshow(mat2gray(sy));title('IDWT of input DWT array');