function F=getwave(wname,Level)
% getwave         Get a synthesis transform matrix based on a dyadic wavelet filter bank
%                 Note the NormLim parameter given in this file, which may be set to
% trucate the filters, i.e. remove ends that are close to zero.
% This function returns the synthesis filters, which are derived from the
% wavelet (reconstruction) filters found by wfilters, as columns in a matrix F. 
% This is like the matrices used for overlapping transforms. For orthogonal filter banks
% (orthogonal wavelets) we have F'*F=I. The normalization of the filters in
% the last part of the function cause problems for biorthogonal wavelets, but
% biorthogonal wavelets are in general problematic.
%
% F=getwave(wname,Level);  
% F=getwave('db3',1);  % returns the filters in 'db3' as columns in F, 6x2 
% F=getwave('db3',3);  % F is size 40x8, N=8 and P=5. 
% F=getwave('db79s',3);G=getwave('db79a',3);  % F'*F~=I, but G'*F=I (and F'*G=I)
% ---------------------------------------------------------------------------
% arguments:
%  F        The matrix with the synthesis filters, size of F (without truncation)
%           is NPxN, where N=2^Level and P depends on the length of the wavelet
%           filters and the number of Levels in the filter bank.
%  wname    name of the wavelet filter, as used in wfilters.m,
%           the synthesis (or analysis) part is returned.
%           also Daubechies 7-9 biorthogonal filter is possible
%  Level    number of levels to use, 
%           use Level=1 to return just the synthesis filters, 
% ---------------------------------------------------------------------------

% The Daubechies 7-9 biorthogonal wavelet is like in table 7.2 page 119
% in C.S.Burrows et.al "Introduction to Wavelets and Wavelet Transforms"   
% It is normalized. Use 'db79a' for analysis and 'db79s' for synthesis
% or only 'db79' (which use the synthesis filters).
% 'db79s' er (nesten) lik 'bior4.4' og 'db79a' er (nesten) lik 'rbio4.4',
% forskjellen er fortegn for det ene filteret og justering i matrisene.

%----------------------------------------------------------------------
% Copyright (c) 2000.  Karl Skretting.  All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail:  karl.skretting@tn.his.no   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  20.10.2000  KS: m-file made
% Ver. 1.1  09.01.2001  KS: added Daubechies 7-9 biorthogonal wavelet
% Ver. 1.2  11.01.2001  KS: added the NormLim parameter
% Ver. 1.3  29.05.2002  KS: Names 'db79a' and 'db79s' are used
% Ver. 1.4  28.11.2002  KS: moved from ..\Frames to ..\FrameTools
% Ver. 1.5  02.10.2008  KS: Exchange synthesis and analysis for db79
% Ver. 1.6  12.03.2010  KS: copied from ..\FrameTools\GetWave.m to
%                       ..\ICTools\getwave.m
% Ver. 1.7  12.03.2012  KS: Added numbers so 'db2', 'db3', 'db4' work even if wavelet
%                       toolbox (wfilters.m) is not available
%----------------------------------------------------------------------
% det er noen problemer med dette her en plass, linje 124 kan få feil
% indekser, TestFrames(x,x,x,[28,29])  30.6.2005 KS


Mfile='getwave';
NormLim=0;            % keeps also small ends, this may lead to F(:,:,p)
                      % is only zeros
%NormLim=1e-4;         % remove small ends of the filters 

J=abs(Level)+1;
L=zeros(J,1);         % length of the different synthesis filters
U=zeros(J,1);         % upsampling factor for the different synthesis filters

if (strcmp(wname,'db79') || strcmp(wname,'db79a') || strcmp(wname,'db79s') || ...
    strcmp(wname,'db97') || strcmp(wname,'db97a') || strcmp(wname,'db97s'))     
   % makes the filters of length 10 (this makes them fit into the F matrix
   % of size 10x2, or 2x2x5. F=[h0',h1']
   h0=[0,0,-0.06453888262876,-0.04068941760920,0.41809227322204,...
         0.78848561640637,0.41809227322204,-0.04068941760920,-0.06453888262876,0];
   h1=[0.03782845550726,0.02384946501956,-0.11062440441844,-0.37740285561283,...
         0.85269867900889,-0.37740285561283,-0.11062440441844,0.02384946501956,...
         0.03782845550726,0];
   if strcmp(wname,'db79a') || strcmp(wname,'db97a') 
      % the inverse filters (analysis)
       h0=[0,0.03782845550726,-0.02384946501956,-0.11062440441844,0.37740285561283,...
             0.85269867900889,0.37740285561283,-0.11062440441844,-0.02384946501956,...
             0.03782845550726];
       h1=[0,0.06453888262876,-0.04068941760920,-0.41809227322204,0.78848561640637,...
             -0.41809227322204,-0.04068941760920,0.06453888262876,0,0];
   end
else
    if exist('wfilters.m','file')
        [h0,h1]=wfilters(wname,'r');   % get the wavelet reconstruction filters
   else
      disp([Mfile,': can not find wfilters.m, is wavelet toolbox installed?']);
      if strcmpi(wname,'db2')
           h0 = [0.482962913144690   0.836516303737469   0.224143868041857  -0.129409522550921];
           h1 = [-0.129409522550921  -0.224143868041857   0.836516303737469  -0.482962913144690];
      elseif strcmpi(wname,'db3')
          h0 = [    0.332670552950957
                    0.806891509313339
                    0.459877502119331
                   -0.135011020010391
                   -0.085441273882241
                    0.035226291882101 ]';
          h1 = [    0.035226291882101
                    0.085441273882241
                   -0.135011020010391
                   -0.459877502119331
                    0.806891509313339
                   -0.332670552950957 ]';
      elseif strcmpi(wname,'db4')
          h0 = [      0.230377813308855
   0.714846570552542
   0.630880767929590
  -0.027983769416984
  -0.187034811718881
   0.030841381835987
   0.032883011666983
  -0.010597401784997 ]';
          h1 = [    -0.010597401784997
  -0.032883011666983
   0.030841381835987
   0.187034811718881
  -0.027983769416984
  -0.630880767929590
   0.714846570552542
  -0.230377813308855 ]';
      else
          disp([Mfile,': use Haar-wavelet = db1 ']);
          h0=[1,1]/sqrt(2);h1=[1,-1]/sqrt(2);    % Haar wavelet
      end
   end
end
Level=abs(Level);

f=[];        % f stores the different filters 
gh=h1;       % the high pass filter
hh=1;        % bare 1 (en) se 10-15 linjer nedenfor 
j=J;
for l=1:Level
    L(j)=length(gh);
    f=[gh(:);f];
    U(j)=2^l;
    j=j-1;
    if 0 % exist('wfilters.m')==2      % wavelet toolbox
        gh=conv(dyadup(gh,0),h0);
        hh=conv(dyadup(hh,0),h0);   % note: l=1 (first time) this is hh=h0
    else
        gh=[gh;zeros(1,length(gh))];   % gh was a row vector
        gh=gh(:)';gh=gh(1:(end-1));  % and is still a row vector
        gh=conv(gh,h0);
        hh=[hh;zeros(1,length(hh))];   % hh was a row vector
        hh=hh(:)';hh=hh(1:(end-1));  % and is still a row vector
        hh=conv(hh,h0);
    end
end
L(1)=length(hh);
f=[hh(:);f];
U(1)=U(2);

N=U(1);       % N is least common multiple of upsampling factors
K=sum(N./U);
P=max(ceil((L-U)./N))+1;  % overlap factor

F=zeros(N*P,K);
Q=length(f);

% top align or center the vectors in F
nj=zeros(J,1);
for j=2:J; nj(j)=floor((L(1)-L(j)-U(1)+U(j))/(2*U(j)))*U(j); end;
k=0;q=1;
for j=1:J
    qi=q:(q+L(j)-1);     % indexes if f
    q=q+L(j);
    g=qi;                % values to write into F
    n=nj(j);
    for k=(k+1):(k+N/U(j))
        F((n+1):(n+L(j)),k)=f(g');
        n=n+U(j);
    end
end

% check if beginning or end of F is small
if NormLim
   NP=size(F,1);
   while norm(F(1:N,:))<NormLim
      F=F((N+1):NP,:);
      NP=size(F,1);
   end
   while norm(F((NP-N+1):NP,:))<NormLim
      F=F(1:(NP-N),:);
      NP=size(F,1);
   end
end   

% now normalize F, but not for the biorthogonal wavelets!
wname=[wname,'    ']; % makes name longer
if ~(strcmp(wname(1:4),'db79') || strcmp(wname(1:4),'db97') || ...
        strcmp(wname(1:4),'bior') || strcmp(wname(1:4),'rbio') )
   for k=1:K;    
      temp=F(:,k)'*F(:,k);
      if (temp > 0); F(:,k)=F(:,k)/sqrt(temp); end;
      if sum(F(:,k))<0; F(:,k)=-F(:,k); end;
   end
end

return

