synth.m

来自「这是伯克里wavelet transforms一书中的例子的代码」· M 代码 · 共 78 行

function sy = synth(a,fopt,lev)
%Routine to reconstruct the signal from its DWT
%using orthonormal filter bank. Circular extensions are 
%used to extend the subbands. This routine is the complement to
%the DTWT analysis routine 'anal'.
%
%sy =synth(an,lpf,L) takes the L level DWT, 'an', which could be 
%the output of the 'anal' routine, and reconstructs the original 
%sequence using the orthonormal synthesis filter bank defined by lpf. 
%If the filter is chosen to correspond to the analysis filter used 
%for DWT,then we get perfect reconstruction. The reconstructed signal 
%is returned in the array 'sy'.
% 
%sy = synth(an,fopt,L) takes the L level DWT, 'an', which could be
%the output of the 'anal' routine, and reconstructs the original 
%signal using the orthonormal synthesis filter bank provided in the 
%routine. This filter bank is chosen using the variable fopt. If the 
%filter bank used here corresponds to the analysis filter bank used 
%for DWT, then we get perfect reconstruction. The reconstructed signal 
%is returned in the array 'sy'.
%
%The final result of the synthesis procedure is plotted in the end.
%
%The input variables are:
%an   : array containing DWT/subband coefficients.
%lpf  : Lowpass filter of the ANALYSIS filter bank.
%OR
%fopt : which can take a value from 1 to 5 and corresponds to
%       Daubechies 2*fopt tap filters respectively. 
%levs : Specifies the number of levels to which the analysis signal
%       has been decomposed.
%
%Refer to Chapter 3 for more information on Orthonormal 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
       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.

  la = length(a);
  ll = la/(2^lev); %initialize
  an =a;           %initialize
  for i=1:lev %reconstruct each level
    sy = an(1:2*ll);
    sy = syn1(sy,lo,hi); %one level reconstruction routine
    an(1:2*ll) = sy(1:2*ll);
    ll = ll*2;
  end; %endfor
  figure;plot(sy);title('IDWT of Input DWT Array');
%plot the reconstructed signal