getwave.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

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 frames. 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).

%----------------------------------------------------------------------
% Copyright (c) 2000.  Karl Skretting.  All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail:  [e-mail]   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
%----------------------------------------------------------------------


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')
   % 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.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];
   if strcmp(wname,'db79a') 
      % the inverse filters (analysis)
      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];
   end
else
    if exist('wfilters.m')==2
        [h0,h1]=wfilters(wname,'r');   % get the wavelet reconstruction filters
   else
      disp([Mfile,': can not find wfilters.m, is wavelet toolbox installed?']);
      h0=[1,1]/sqrt(2);h1=[1,-1]/sqrt(2);    % Haar wavelet
   end
end
Level=abs(Level);

f=[];        % f stores the different filters 
gh=h1;       % the high pass filter
hh=1;
j=J;
for l=1:Level
   L(j)=length(gh);
   f=[gh(:);f];
   U(j)=2^l;
   j=j-1;
   if 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(length(gh))];   % gh was a row vector
      gh=gh(:)';gh=gh(1:(end-1));  % and is still a row vector
      hh=[hh;zeros(length(hh))];   % hh was a row vector
      hh=hh(:)';hh=hh(1:(end-1));  % and is still a row vector
   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),'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