%% New Wavelet Permutation Ordering rversed:  iperm2.m 
% Usage:  [Ao,x] = Mr_iperm(Ai,d1,d2,lev)
%    or      Ao  = Mr_iperm(Ai,d1,d2,lev) for x ordering
%
%e.g. A=cz(64,3); close all  % 64 x 64 CZ matrix %             from FWT2.m 
%     A1=fwt(A,4, 1,1,3);  subplot(221);spyc(A1)
%     A2=perm2(A1, 1,1,3); subplot(222);spyc(A2),title(''), hold on
%    Ab=iperm2(A2, 1,1,3); xlabel('');subplot(212);spyc(Ab);hold off
%    check=norm(A1-Ab,1),  title('Undo the perm2')
 
function [Ao,x] = Mr_perm(Ai,d1,d2,lev)
      %% n=input('Size N = ? '); lev=input('DWT levels = ? ');
if nargin ==1,
   d1 = 1; d2 = 1;
end
if nargin <1, help iperm2; return; end
if nargin <4, lev=floor(log(n)/log(2)); end
     [m1 m2] = size(Ai); n=max([m1 m2]);
       if m1 ==1, d1=0; end
         if m2 ==1, d2=0; end
    p=speye(n);           x = 1:n; n1=n;
  Pat='A=PAP'''; if d2==0; Pat='A=PA'; end, if d1==0; Pat='A=AP'''; end 
fprintf('  iperm2.m for %d levels of A permutation %d x %d [%s]\n\t\tDone levs:',...
               lev,m1,m2, Pat);

for k=lev:-1:1     
    n1 = floor(n/2^k)*2;         % Active length- allow odd n
    n2 = floor(n/2^k);           % Active's Half
         x0 = x(1:n1);   % Part 1 of x
         x2 = x(n1+1:n); % Part 2 of x - fixed
    x1 = [];
    for j=1:n2; x1=[x1 j n2+j]; end
    x1 = x0(x1);   %% Re-order this Part 1
    x = [x1 x2];   %% Re-set x
    fprintf(' %d',k)
end,disp(' ')
%  P=p(:,x);  % for A_new = P*A_dwt*P'
%  Ao = Ai;
%     if d1==1, Ao = P*Ao; end
%     if d2==1, Ao = Ao*P'; end
  Ao = Ai; [x,j,s]=find(p(x,:)); %% NEW x found!!!
     if d1==1, Ao = Ao(x,:); end
     if d2==1, Ao = Ao(:,x); end
%   check = norm(An-Ao,1)
%--------------------------------------------------------------
