%% New Wavelet Permutation Ordering  perm2.m
% Usage:  [Ao,x] = Mr_perm(Ai,d1,d2,lev)
%    or      Ao  = Mr_perm(Ai,d1,d2,lev) for x ordering
%
%e.g. A=cz(64,3);  % 64 x 64 CZ matrix %                        from FWT2.m 
%     A1=fwt(A,4, 1,1,3); spyc(A1), figure
%    A2=perm2(A1, 1,1,3); spyc(A2)
 
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 perm2; 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=P''AP'; if d2==0; Pat='A=P''A'; end, if d1==0; Pat='A=AP'; end 
fprintf('  perm2.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;
   if d1==1, Ao = Ao(x,:); end
   if d2==1, Ao = Ao(:,x); end
   %% check = norm(An-Ao,1)
%--------------------------------------------------------------
