%% iwts.m --- for Daubechies' order m inverse FWT transforms (sparse) vs iwt.m
%% Usage: [1] y =iwts(x)           FWT for vector x
%%        [2] y =iwts(x,m)         FWT for vector x (m order)
%%        [3] B =iwts(A)           FWT for matrix A
%%        [4] B =iwts(A,m)         FWT for matrix A (m order)
%%        [5] B =iwts(A,m,d1,d2)   FWT for matrix A with 1-sided options
%%        [6] B =iwts(A,m,d1,d2,t) for t-levels of FWT for matrix A
%%                                [dj=1 on,=0 off] with 1-sided options 
%%     [7] B =iwts(A,m,d1,d2,t,s,iout) for t-levels and s=0 (normal) 1 (NS)
%%e.g. A=cz(128,2);                    % iout = 0 (default) or =1 (printout)
%      B = iwts(A,4, 1, 1, 3); % for 3-levels of Daub_4 (both sides) A=A*W^T 
%      B2= iwts(A,6, 1, 0, 2); % for 2-levels of Daub_6 (left only)  A=W*A
function M=iwts(A,m_in,d1,d2,t,formy,iout) % formy 0/normal 1/NS
if nargin <1, eval(['help ' mfilename]), 
              disp('! At least one parameter needed');  return, end
   [n1 n2]=size(A); n=max(n1,n2); M=A;
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
if nargin <6,formy=0; end
if nargin <7, iout=0; end
   t_max=round( log(n)/log(2) );          % n = 2^t_max
if nargin <5  | length(t)<1 %-------------------- Set t
   t=t_max+1 - ceil( log(m_in)/log(2) );  % Max active levels up to m
end
if iout==1
   Pat='A=W''AW'; if d2==0; Pat='A=W''A'; end, if d1==0; Pat='A=AW'; end
fprintf(' iwts for %d levels of Daub%d to A of %d x %d [%s,%d]\n\t\tDone levs:',...
               t,m_in,n1,n2, Pat,formy);
end
N = 2*floor(n/2); 

for iter = t:-1:1
k=floor( n/2^iter ) * 2; 
  T=speye(n,n); p=1:n;
% NS form ----------------------------------------
if formy==1
  T(N-k+1:N,N-k+1:N)=w_matrix(m_in,k);
  p(N-k+1:2:N)=N+1-k/2:N;     % Odd ...... | A B |
  p(N-k+2:2:N)= N-k+1:N-k/2;  % Even ..... | C T |
else
% Normal form ------------------------------------
  T(1:k,1:k)=w_matrix(m_in,k);
  p(1:2:k)= 1:k/2;     % Odd ...... | T B |
  p(2:2:k)= k/2+1:k;   % Even ..... | C A |
end

  if (d1==1 &  d2==1),M=M(p,p);end
  if (d1==1 &  d2==0),M=M(p,:);end
  if (d1==0 &  d2==1),M=M(:,p);end
  
  if (d1==1),M=T'*M; end
  if (d2==1),M=M*T;end

  k1=k/2; if k1*2 ~= k, iter=1; end
  if iout==1, fprintf(' %d',iter); end
end % iter - t levels;
  if iout==1, disp(' '), end
%% ________________________________________________________
%% function W=w_matrix(m, n) - matrix size n with m filters
function W = w_matrix(m,n)
   p = find([2 4 6 8 16] == m);
if nargin<1, help w_matrix; return; end
if m<2, W=speye(n); return; end
 %  Note these DIAG repeating patterns:   (0,c_{m-1}), (c_{m-2} d_{m-1})
 %   (c_{m-4} d_{m-3}) (c_{m-6} d_{m-5}) ... (c_{0} d_1) and (c_{m-1} 0)
if m==2
c1=sqrt(2)/2; c0=sqrt(2)/2; 
d1=zeros(n,1); d1(1:2:n)=-c0*ones(n/2,1);
d3=c0*ones(n,1);
d5=zeros(n,1); d5(2:2:n)= c0*ones(n/2,1);
W=spdiags([d1 d3 d5],-1:1,n,n);
elseif m==4      %====================================================
c3=sqrt(2)*(1+sqrt(3))/8; c2=sqrt(2)*(3+sqrt(3))/8;
c1=sqrt(2)*(3-sqrt(3))/8; c0=sqrt(2)*(1-sqrt(3))/8;

d1=zeros(n,1); d1(1:2:n)=c0*ones(n/2,1);
d2=zeros(n,1); d2(1:2:n)=c3*ones(n/2,1); d2(2:2:n)=-c1*ones(n/2,1);
d3=c2*ones(n,1);
d4=zeros(n,1); d4(1:2:n)=c1*ones(n/2,1); d4(2:2:n)=-c3*ones(n/2,1);
d5=zeros(n,1); d5(2:2:n)=c0*ones(n/2,1);

W=spdiags([d1 d2 d3 d4 d5],-1:3,n,n);
W(n-1:n,1:2)=[c1 c0; c2 -c3];
elseif m==6      %====================================================
c0= 3.326705529500826d-01; c1= 8.068915093110927d-01;
c2= 4.598775021184915d-01; c3=-1.350110200102546d-01;
c4=-8.544127388202666d-02; c5= 3.522629188570955d-02;

d1=zeros(n,1); d1(1:2:n)=c5*ones(n/2,1);
d2=zeros(n,1); d2(1:2:n)=c0*ones(n/2,1); d2(2:2:n)=-c4*ones(n/2,1);
d3=zeros(n,1); d3(1:2:n)=c3*ones(n/2,1); d3(2:2:n)= c1*ones(n/2,1);
d4=zeros(n,1); d4(1:2:n)=c2*ones(n/2,1); d4(2:2:n)=-c2*ones(n/2,1);
d5=zeros(n,1); d5(1:2:n)=c1*ones(n/2,1); d5(2:2:n)= c3*ones(n/2,1);
d6=zeros(n,1); d6(1:2:n)=c4*ones(n/2,1); d6(2:2:n)=-c0*ones(n/2,1);
d7=zeros(n,1);                           d7(2:2:n)= c5*ones(n/2,1);

W=spdiags([d1 d2 d3 d4 d5 d6 d7],-1:5,n,n);
W(n-3:n-2,1:2)=[c4 c5; c1 -c0];
W(n-1:n,3:4)=[c4 c5; c1 -c0];
W(n-1:n,1:2)=[c2 c3; c3 -c2];
elseif m==8      %==================================================== 
c0=-7.5765714789357d-2; c1=-2.9635527645960d-2;
c2= 4.9761866763256d-1; c3= 8.0373875180539d-1;
c4= 2.9785779560561d-1; c5=-9.9219543576956d-2;
c6=-1.2603967262264d-2; c7= 3.2223100604078d-2;

d1=zeros(n,1); d1(1:2:n)=c7*ones(n/2,1);
d2=zeros(n,1); d2(1:2:n)=c0*ones(n/2,1); d2(2:2:n)=-c6*ones(n/2,1);
d3=zeros(n,1); d3(1:2:n)=c5*ones(n/2,1); d3(2:2:n)= c1*ones(n/2,1);
d4=zeros(n,1); d4(1:2:n)=c2*ones(n/2,1); d4(2:2:n)=-c4*ones(n/2,1);
d5=c3*ones(n,1); 
d6=zeros(n,1); d6(1:2:n)=c4*ones(n/2,1); d6(2:2:n)=-c2*ones(n/2,1);
d7=zeros(n,1); d7(1:2:n)=c1*ones(n/2,1); d7(2:2:n)= c5*ones(n/2,1);
d8=zeros(n,1); d8(1:2:n)=c6*ones(n/2,1); d8(2:2:n)=-c0*ones(n/2,1);
d9=zeros(n,1);                           d9(2:2:n)= c7*ones(n/2,1);

W=spdiags([d1 d2 d3 d4 d5 d6 d7 d8 d9],-1:7,n,n);
W(n-5:n-4,1:2)=[c6 c7; c1 -c0];
W(n-3:n-2,3:4)=[c6 c7; c1 -c0];
W(n-1:n  ,5:6)=[c6 c7; c1 -c0];
W(n-3:n-2,1:2)=[c4 c5; c3 -c2];
W(n-1:n,3:4)=[c4 c5; c3 -c2];
W(n-1:n,1:2)=[c2 c3; c5 -c4];
elseif m==8      %==================================================== 
    c0  = 1.8899503329009d-03; c1  =-3.0292051455166d-04;
    c2  =-1.4952258336794d-02; c3  = 3.8087520140604d-03;
    c4  = 4.9137179673481d-02; c5  =-2.7219029916816d-02;
    c6  =-5.1945838107879d-02; c7  = 3.6444189483599d-01;
    c8  = 7.7718575169981d-01; c9  = 4.8135965125924d-01;
    c10 =-6.1273359067914d-02; c11 =-1.4329423835107d-01;
    c12 = 7.6074873252854d-03; c13 = 3.1695087810348d-02;
    c14 =-5.4213233163559d-04; c15 =-3.3824159513597d-03;
d1=zeros(n,1); d1(1:2:n)=c15*ones(n/2,1);
d2=zeros(n,1); d2(1:2:n)=c0 *ones(n/2,1); d2(2:2:n)=-c14*ones(n/2,1);
d3=zeros(n,1); d3(1:2:n)=c13*ones(n/2,1); d3(2:2:n)= c1 *ones(n/2,1);
d4=zeros(n,1); d4(1:2:n)=c2 *ones(n/2,1); d4(2:2:n)=-c12*ones(n/2,1);
d5=zeros(n,1); d5(1:2:n)=c11*ones(n/2,1); d5(2:2:n)=c3  *ones(n/2,1);
d6=zeros(n,1); d6(1:2:n)=c4 *ones(n/2,1); d6(2:2:n)=-c10*ones(n/2,1);
d7=zeros(n,1); d7(1:2:n)=c9 *ones(n/2,1); d7(2:2:n)= c5 *ones(n/2,1);
d8=zeros(n,1); d8(1:2:n)=c6 *ones(n/2,1); d8(2:2:n)=-c8 *ones(n/2,1);
d9=c7*ones(n,1); 
d10=zeros(n,1); d10(1:2:n)=c8 *ones(n/2,1); d10(2:2:n)=-c6*ones(n/2,1);
d11=zeros(n,1); d11(1:2:n)=c5 *ones(n/2,1); d11(2:2:n)=c9 *ones(n/2,1);
d12=zeros(n,1); d12(1:2:n)=c10*ones(n/2,1); d12(2:2:n)=-c4*ones(n/2,1);
d13=zeros(n,1); d13(1:2:n)=c3 *ones(n/2,1); d13(2:2:n)=c11*ones(n/2,1);
d14=zeros(n,1); d14(1:2:n)=c12*ones(n/2,1); d14(2:2:n)=-c2*ones(n/2,1);
d15=zeros(n,1); d15(1:2:n)=c1 *ones(n/2,1); d15(2:2:n)=c13*ones(n/2,1);
d16=zeros(n,1); d16(1:2:n)=c14*ones(n/2,1); d16(2:2:n)=-c0*ones(n/2,1);
d17=zeros(n,1);                             d17(2:2:n)=c15*ones(n/2,1);

W=spdiags([d1 d2 d3 d4 d5 d6 d7 d8 d9 d10 d11 d12 d13 d14 d15 d16 d17],-1:15,n,n);
W=W+spdiags([d3 d4 d5 d6 d7 d8 d9 d10 d11 d12 d13 d14 d15 d16 d17],-n+1:-n+15,n,n);
end
