%% ifwt.m --- for Daubechies' order m IFWT transforms
%% Usage: [1]   y =ifwt(x)           IFWT for vector x
%%        [2]   y =ifwt(x,m)         IFWT for vector x (m order)
%%        [3]   B =ifwt(A)           IFWT for matrix A
%%        [4]   B =ifwt(A,m)         IFWT for matrix A (m order)
%%        [5]   B =ifwt(A,m,d1,d2)   IFWT for matrix A with 1-sided options
%%        [6]   B =ifwt(A,m,d1,d2,t) for t-levels of IFWT for matrix A
%%                                  [dj=1 on,=0 off] with 1-sided options 
%% e.g  B = ifwt(A,6, 0, 1, 3) for 3-levels of Daub_6 (right only) A=A*W
%%      B = ifwt(A,6, 1, 0, 2) for 2-levels of Daub_6 (left only) A=W'*A
function M=ifwt(A,m_in,d1,d2,t)
%_______________________________________ Internal Control
if nargin <1, eval(['help ' mfilename]), 
              disp('! At least one parameter needed');  return, end
   [n1 n2]=size(A); n=max(n1,n2); 
if nargin <2, m_in=4; end %------- Set Daub4 (as default)
if nargin <3 %-------------------- Set d1 / d2
   d1=1; d2=1; 
   if n1==1, d1=0; d2=1; end
   if n2==1, d2=0; d1=1; end
end
   t_max=round( log(n)/log(2) );          % n = 2^t_max
if nargin <5 %-------------------- Set t
   t=t_max+1 - ceil( log(m_in)/log(2) );  % Max active levels up to m
end
   Pat='A=W''AW'; if d2==0; Pat='A=W''A'; end, if d1==0; Pat='A=AW'; end
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Internal Control
fprintf(' ifwt.m for %d levels of Daub%d to A of %d x %d [%s]\n\t\tDone levs:',...
               t,m_in,n1,n2, Pat);

  [c d]=get_filters(m_in);

for iter = t:-1:1

  k=floor( n/2^iter ); %k=Half way  % Allow A to have odd dimension=n
  K=2*k;  %Whole block  (m_in alyays even)

if d1==1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Left W'T
  for col=1:n2
     Sum=zeros(K,1);  % tmp storage for whole col

   for row=1:2:K        %%%%% even c / d Sums  (k-1)=(row-1)

     k_bar=row-1;          % Normal range (row+1-m)/2  to  (row/2)
     pos = 1 + (k_bar-(m_in-2)) / 2; % far left position
     if pos <=0,           % as (row-Shift) will be wrap around
        j1=1;  j2=(m_in-1)-2*(pos+1);    j3=j2+2; j4=m_in; Shift=k;
     else
        j1=1; j2=m_in;  j3=1; j4=0; Shift=0;
     end 
     for j=j1:2:j2  % c_j/d_j even Sums
       pos = 1 + (k_bar-(j-1))/2;   %\ell=k-1
       Sum(row) = Sum(row) + c(j)*A(pos,col); 
       Sum(row) = Sum(row) + d(j)*A(pos+k,col); 
     end
     for j=j3:2:j4  % c_j/d_j even Sums
       pos = 1 + (k_bar-(j-1))/2 + Shift;   %\ell=k-1
       Sum(row) = Sum(row) + c(j)*A(pos,col); 
       Sum(row) = Sum(row) + d(j)*A(pos+k,col); 
     end

   end %row even
   for row=2:2:K        %%%%%  odd c / d Sums  (k-1)=(row-1)

     k_bar=row-1;          % Normal range (row+1-m)/2  to  (row/2)
     pos = 1 + (k_bar-(m_in-1)) / 2; % far left position
     if pos <=0,           % as (row-Shift) will be wrap around
        j1=2;  j2=m_in-2*(pos+1);    j3=j2+2; j4=m_in; Shift=k;
     else
        j1=2; j2=m_in;  j3=1; j4=0; Shift=0;
     end 
     for j=j1:2:j2  % c_j/d_j  odd Sums
       pos = 1 + (k_bar-(j-1))/2;   %\ell=k-1
       Sum(row) = Sum(row) + c(j)*A(pos,col); 
       Sum(row) = Sum(row) + d(j)*A(pos+k,col); 
     end
     for j=j3:2:j4  % c_j/d_j  odd Sums
       pos = 1 + (k_bar-(j-1))/2 + Shift;   %\ell=k-1
       Sum(row) = Sum(row) + c(j)*A(pos,col); 
       Sum(row) = Sum(row) + d(j)*A(pos+k,col); 
     end

   end %row  odd 

   A(1:K,col) = Sum;  % Updated - Done Col

   end %col

end % d1=1
if d2==1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Right WT i.e. A=AW

  for row=1:n1
     Sum=zeros(1,K);  % tmp storage for whole row

   for col=1:2:K        %%%%% even c / d Sums  (k-1)=(col-1)

     k_bar=col-1;          % Normal range (col+1-m)/2  to  (col/2)
     pos = 1 + (k_bar-(m_in-2)) / 2; % far left position
     if pos <=0,           % as (col-Shift) will be wrap around
        j1=1;  j2=(m_in-1)-2*(pos+1);    j3=j2+2; j4=m_in; Shift=k;
     else
        j1=1; j2=m_in;  j3=1; j4=0; Shift=0;
     end 
     for j=j1:2:j2  % c_j/d_j even Sums
       pos = 1 + (k_bar-(j-1))/2;   %\ell=k-1
       Sum(col) = Sum(col) + c(j)*A(row,pos); 
       Sum(col) = Sum(col) + d(j)*A(row,pos+k); 
     end
     for j=j3:2:j4  % c_j/d_j even Sums
       pos = 1 + (k_bar-(j-1))/2 + Shift;   %\ell=k-1
       Sum(col) = Sum(col) + c(j)*A(row,pos); 
       Sum(col) = Sum(col) + d(j)*A(row,pos+k); 
     end

   end %row even
   for col=2:2:K        %%%%%  odd c / d Sums  (k-1)=(col-1)

     k_bar=col-1;          % Normal range (col+1-m)/2  to  (col/2)
     pos = 1 + (k_bar-(m_in-1)) / 2; % far left position
     if pos <=0,           % as (col-Shift) will be wrap around
        j1=2;  j2=m_in-2*(pos+1);    j3=j2+2; j4=m_in; Shift=k;
     else
        j1=2; j2=m_in;  j3=1; j4=0; Shift=0;
     end 
     for j=j1:2:j2  % c_j/d_j  odd Sums
       pos = 1 + (k_bar-(j-1))/2;   %\ell=k-1
       Sum(col) = Sum(col) + c(j)*A(row,pos); 
       Sum(col) = Sum(col) + d(j)*A(row,pos+k); 
     end
     for j=j3:2:j4  % c_j/d_j  odd Sums
       pos = 1 + (k_bar-(j-1))/2 + Shift;   %\ell=k-1
       Sum(col) = Sum(col) + c(j)*A(row,pos); 
       Sum(col) = Sum(col) + d(j)*A(row,pos+k); 
     end

   end %col  odd 

   A(row,1:K) = Sum;  % Updated - Done Col

   end %row

end % id=2
   fprintf(' %d',iter) 
end % iter - t levels;
   M = A;    disp(' ')
return

%===================================================================
function [c,d] = get_filters(m)
   p = find([2 4 6 8 16] == m);
if isempty(p)
   m,error('Wrong m asked or Requested m not yet implemeneted');
   return
end,  c=[];d=[]; Sign=-1;
if p==1
    c(1)=sqrt(2)/2;
    c(2)=sqrt(2)/2; 
elseif p==2
    c(1)=sqrt(2)*(1+sqrt(3))/8;
    c(2)=sqrt(2)*(3+sqrt(3))/8;
    c(3)=sqrt(2)*(3-sqrt(3))/8;
    c(4)=sqrt(2)*(1-sqrt(3))/8;
elseif p==3
    c(1)= 3.326705529500826d-01;
    c(2)= 8.068915093110927d-01;
    c(3)= 4.598775021184915d-01;
    c(4)=-1.350110200102546d-01;
    c(5)=-8.544127388202666d-02;
    c(6)= 3.522629188570955d-02;
elseif p==4
    c(1)=-7.5765714789357d-2;
    c(2)=-2.9635527645960d-2; 
    c(3)= 4.9761866763256d-1;
    c(4)= 8.0373875180539d-1;
    c(5)= 2.9785779560561d-1;
    c(6)=-9.9219543576956d-2;
    c(7)=-1.2603967262264d-2;
    c(8)= 3.2223100604078d-2;
elseif p==5
    c(1) = 1.8899503329009d-03;
    c(2) =-3.0292051455166d-04;
    c(3) =-1.4952258336794d-02;
    c(4) = 3.8087520140604d-03;
    c(5) = 4.9137179673481d-02;
    c(6) =-2.7219029916816d-02;
    c(7) =-5.1945838107879d-02;
    c(8) = 3.6444189483599d-01;
    c(9) = 7.7718575169981d-01;
    c(10)= 4.8135965125924d-01;
    c(11)=-6.1273359067914d-02;
    c(12)=-1.4329423835107d-01;
    c(13)= 7.6074873252854d-03;
    c(14)= 3.1695087810348d-02;
    c(15)=-5.4213233163559d-04;
    c(16)=-3.3824159513597d-03;
else
    error('More filters to be added')
end
%========= Set d
    for j=1:m
      d(m+1-j)=Sign*c(j); Sign=-Sign;
    end
return
