%% New Wavelet Permutation Ordering   perm0.m
% Usage:  [Ao,x] = perm0(Ai,d1,d2,lev,lev_border,id)
%    or      Ao  = perm0(Ai,d1,d2,lev)      % Exclude the lev_border
%                  for  x = ordering  and lev_border in [1,lev] (id=1 extra)
% If lev-border=lev, perm0 is the same as DWTPer (perm2) !
% If lev-border=0,   perm0 is the same as DWTPerMod (omit J1)
%
%e.g.  n=128,A=cz(n,1);  % n x n CZ matrix %                      from perm2.m 
%        A1=fwt(A,4, 1,1,4); spy(abs(A1)>1/n), figure
%       A0=perm0(A1, 1,1,4); spy(abs(A0)>1/n), figure,
%       A0=perm0(A1, 1,1,4,0); spy(abs(A0)>1/n), figure
%       A0=perm0(A1, 1,1,4,1); spy(abs(A0)>1/n), figure
%       A0=perm0(A1, 1,1,4,2); spy(abs(A0)>1/n) 
 
function [Ao,x] = perm0(Ai,d1,d2,lev,levuse,L2)
      
if nargin ==1, d1 = 1; d2 = 1; end
if nargin <1, help perm0; return; end
     [m1 m2] = size(Ai); n=max([m1 m2]);
     levmax=floor(log(n)/log(2))-1;
if nargin <4, lev=ceil(log(n)/log(2)); end
if nargin <5, levuse=lev; end
       if m1 ==1, d1=0; end
         if m2 ==1, d2=0; end
n_sav = n; 
exclud = n/2^lev; %n=n-exclud;
    p=speye(n_sav);           x = 1:n; n1=n;
  Pat='A=WAW'''; if d2==0; Pat='A=PA'; end, if d1==0; Pat='A=AP'''; end 
fprintf('  perm0.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_sav/2^k)*2;         % Active length- allow odd n
    n2 = floor(n_sav/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
Ao = Ai;
P=p(:,x);  % for A_new = P'*A_dwt*P
   if d1==1, Ao = P'*Ao; end
   if d2==1, Ao = Ao*P; end
%%% ========================= Possibly A Dense Block at top-left corner :
if nargin<6, L2=0; end
if nargin>4
fprintf('  [used border levels = %d]\n',levuse)
   kb=1:2^(lev-levuse):n; %%%%%%%%%%%%

if L2>0 % Once more recursion (could carry on)
   lb=length(kb); 
   [A1 x]=perm0(Ai(kb,kb),d1,d2,lev,levuse); kb=kb(x);
end
   k1=setdiff(1:n,kb); 
   x=[kb k1]; %=====================================
   P=p(:,x);  % for A_new = P'*A_dwt*P
   if d1==1, Ao = P'*Ao; end,  if d2==1, Ao = Ao*P; end
else,disp(' '),end
