%% fwt.m --- for Daubechies' order m FWT transforms (see also sparse fwts.m)
%% Usage: [1]   y =fwt(x)           FWT for vector x
%%        [2]   y =fwt(x,m)         FWT for vector x (m order)
%%        [3]   B =fwt(A)           FWT for matrix A
%%        [4]   B =fwt(A,m)         FWT for matrix A (m order)
%%        [5]   B =fwt(A,m,d1,d2)   FWT for matrix A with 1-sided options
%%        [6]   B =fwt(A,m,d1,d2,t) for t-levels of FWT for matrix A
%%                                  [dj=1 on,=0 off] with 1-sided options 
%%e.g. A=cz(128,2);
%      B = fwt(A,4, 1, 1, 3); % for 3-levels of Daub_4 (both sides) A=A*W^T 
%      B2= fwt(A,6, 1, 0, 2); % for 2-levels of Daub_6 (left only)  A=W*A
function M=fwt(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=WAW'''; if d2==0; Pat='A=WA'; end, if d1==0; Pat='A=AW'''; end
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Internal Control
fprintf('  fwt.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 = 1:t

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

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

   for row=1:k        %%%%% Sums for c 

     k_bar=row-1;          % Normal range (2*row-1)   to    (2*row+m-2)
     pos = m_in + 2*k_bar; % i.e. start=2*(row-1) +1  to    2*(row-1)+m
     if pos <= K,
        j1=1; j2=m_in;   j3=1; j4=0; Shift=0;
     else
        j1=1; j2=m_in - (pos-K);
        j3=j2+1;  j4=m_in; Shift=K;  % as (row-Shift) will be wrap around
     end 
     for j=j1:j2  % c_j Sums
       pos = j + 2*k_bar; 
       Sum(row) = Sum(row) + c(j)*A(pos,col); 
     end
     for j=j3:j4  % c_j Sums
       pos = j + 2*k_bar - Shift; 
       Sum(row) = Sum(row) + c(j)*A(pos,col); 
     end

   end %row

   for f_2nd=(1+k):K    %%%%% Differences d
       row = f_2nd - k;

     k_bar=row-1;          % Normal range (2*row-1)   to    (2*row+m-2)
     pos = m_in + 2*k_bar; % i.e. start=2*(row-1) +1  to    2*(row-1)+m
     if pos <= K,
        j1=1; j2=m_in;   j3=1; j4=0; Shift=0;
     else
        j1=1; j2=m_in - (pos-K);
        j3=j2+1;  j4=m_in; Shift=K;  % as (row-Shift) will be wrap around
     end 
     for j=j1:j2  % d_j Difs
       pos = j + 2*k_bar; 
       Sum(f_2nd) = Sum(f_2nd) + d(j)*A(pos,col); 
     end
     for j=j3:j4  % d_j Difs
       pos = j + 2*k_bar - Shift; 
       Sum(f_2nd) = Sum(f_2nd) + d(j)*A(pos,col); 
     end

   end %f_2nd

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

   end %col

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

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

   for col=1:k        %%%%% Sums for c 

     k_bar=col-1;          % Normal range (2*col-1)   to    (2*col+m-2)
     pos = m_in + 2*k_bar; % i.e. start=2*(col-1) +1  to    2*(col-1)+m
     if pos <= K,
        j1=1; j2=m_in;   j3=1; j4=0; Shift=0;
     else
        j1=1; j2=m_in - (pos-K);
        j3=j2+1;  j4=m_in; Shift=K;  % as (col-Shift) will be wrap around
     end 
     for j=j1:j2  % c_j Sums
       pos = j + 2*k_bar; 
       Sum(col) = Sum(col) + c(j)*A(row,pos); 
     end
     for j=j3:j4  % c_j Sums
       pos = j + 2*k_bar - Shift; 
       Sum(col) = Sum(col) + c(j)*A(row,pos); 
     end

   end %col

   for f_2nd=(1+k):K    %%%%% Differences d
       col = f_2nd - k;

     k_bar=col-1;          % Normal range (2*col-1)   to    (2*col+m-2)
     pos = m_in + 2*k_bar; % i.e. start=2*(col-1) +1  to    2*(col-1)+m
     if pos <= K,
        j1=1; j2=m_in;   j3=1; j4=0; Shift=0;
     else
        j1=1; j2=m_in - (pos-K);
        j3=j2+1;  j4=m_in; Shift=K;  % as (col-Shift) will be wrap around
     end 
     for j=j1:j2  % d_j Difs
       pos = j + 2*k_bar; 
       Sum(f_2nd) = Sum(f_2nd) + d(j)*A(row,pos); 
     end
     for j=j3:j4  % d_j Difs
       pos = j + 2*k_bar - Shift; 
       Sum(f_2nd) = Sum(f_2nd) + d(j)*A(row,pos); 
     end

   end %f_2nd

   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 implemented');
   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
