vsblock.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

function W=VSblock(X,Fin,S,VSalg,B,el,Bf)
% VSblock   Vector Selection for block-oriented frame, 
%           several (=B) blocks may be done at each step (each call to VSalg)
% At each step we select vectors for B signal blocks, the actual frame used is 
% a block-diagonal matrix formed by repeating input frame B times.
% The program works for all block-oriented frames, the signal may be 1D, 2D or
% multi-dimensional, but it must be organized as a matrix of size NxL.
% For overlapping frames (1D signal) use the function VSolap1 or BlockVS
%
% W=VSblock(X,F,S,VSalg,B);       
% W=VSblock(X,F,S,VSalg,B,el,Bf);       
%--------------------------------------------------------------------------------
% arguments:
%  W     - The weight matrix, W is a sparse matrix of size KxL
%  X     - The signal reshaped into blocks of length N, size NxL
%  F     - The frame of synthesis vectors (dictionary), size NxKx1  (P==1)
%          the vectors of F must be normalized, F(:,k)'*F(:,k)==1 for all k
%  S     - the third input argument may have different meaning depending on its size.
%    1x1   S is scalar, then it is assumed to be sparsness factor, 0<S<1,
%          a total of S*N*L weights will be selected, evenly distributed or by GMP.
%    1xL   number of vectors to select for each block of length N 
%          S(l) is used for vector X(:,l) to find W(:,l)
%          Note: only when B==1, if B>1 the distribution of weights may be changed
%    KxL   and the third argument should be the previous weights, W.
%          Now S=full(sum(W~=0)); and used as previous case (size 1xL).
%          If VSalg='VSab2' previous weights are used as input when VSalg is called
%  VSalg - Which vector selection algorithm to use, it must be a function
%          called like 'w=VSfomp2(x,S);' and have F and FF as as global variables. 
%          Recommended (=tested) are 'VSfomp2', 'VSmp2' or 'VSab2'.
%  B     - number of blocks to use at each call to VSalg, the actual frame will
%          be a block-diagonal matrix of size NBxKB, where the input frame (size NxK)
%          is used in the blocks on the diagonal.
%  el    - extra loops to do, vector selection will be done (1+el) times 
%          for each block,   default is el=0
%  Bf    - The increment of l when going from one block to the next
%          We should have (1 <= Bf <= B), and Bf a factor of L, default is Bf=B
%--------------------------------------------------------------------------------

%--------------------------------------------------------------------------------
% Copyright (c) 2000.  Karl Skretting.  All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail:  [e-mail]   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  04.10.2001  KS: made function based on BlockVS (based on FindW)
% Ver. 1.1  26.10.2001  KS: some minor changes, el and Bf arguments added
% Ver. 1.2  25.11.2002  KS: moved from ..\Frames to ..\FrameTools
%--------------------------------------------------------------------------------

global F FF 
Mfile='VSblock';
Display=0;        % decide if function display information, 2 also waitbar
UseGMP=0;         % GMP: Global Matching Pursuit

if (nargin<5)
   error([Mfile,' should have at least 5 arguments, see help.']);
end
[k,Ls]=size(S);
[n,L]=size(X);
[N,K,P]=size(Fin);           % note: input frame is called Fin, not F
if (k==K) & (Ls==L)
   W=S;             % the old (previous) weights 
   S=full(sum(W~=0));
   [k,Ls]=size(S);
else
   W=sparse(K,L);    % the weight matrix is sparse
end
if P>1
   error([Mfile,': P>1, VSolap1 or VSolap2 should be used for overlapping frame.']);
end
if k~=1
   error([Mfile,': size of third argument (S or W) is wrong, see help.']);
end
if n~=N
   error([Mfile,': size of X and F do not correspond, see help.']);
end   
if Ls==1   % S is sparseness factor, 0<S<1
   Savg=S;
   S=zeros(1,L);
   if UseGMP
      R=X;
      Q=Fin'*R;                                               % Q is KxL
      SNL=floor(Savg*N*L);     % number of non-zeros in W left to select
      L1=(L-max([floor(Savg*L),1])+1);
      count=0;
      while SNL>0
         count=count+1;
         [temp,k]=max(abs(Q));   
         [temp,l]=sort(temp);    
         k=k(l);
         for i=L:(-1):L1
            ki=k(i);li=l(i);  %  k(i) and l(i) are the indexes in Q (and W)
            if ~W(ki,li); SNL=SNL-1; end;
            if (SNL<0); break; end;
            W(ki,li)=W(ki,li)+Q(ki,li);
            R(:,li)=R(:,li)-Fin(:,ki)*Q(ki,li);
            Q(:,li)=Fin'*R(:,li);
         end
      end   
      clear R Q SNL L1 count temp k l ki li i
      S=full(sum(W~=0));
      % return     % this may be removed (or included)
   else
      % distrubute S*N*L evenly
      t1=0;
      for l=1:L
         t1=t1+N*Savg;
         S(l)=floor(t1);
         t1=t1-S(l);
      end
   end
   [k,Ls]=size(S);
end
if Ls~=L
   error([Mfile,': size of X and S do not correspond, see help.']);
end   

if (nargin<6); el=0; end;
if (nargin<7); Bf=B; end;
if length(el)==0; el=0; end;
if length(Bf)==0; Bf=B; end;
if (nargout < 1); 
   error([Mfile,': function must have one output arguments, see help.']); 
end
if (el>20); el=20; end;   % not too many extra loops
if (B==0); B=1; end;
if (Bf==0); Bf=B; end;

% now build the frame to use here from input frame
% F and FF are global variables from now on
F=zeros(N*B,K*B);
for b=0:(B-1)
   F((1:N)+b*N,(1:K)+b*K)=Fin;
end
FF=F'*F;

Stot=sum(S);  % S is now 1xL

if Display
   disp([Mfile,': size of frame is ',int2str(N*B),'x',int2str(K*B),...
         ' and average number of weights is ',int2str(floor(Stot*B/L+0.5))]);
   disp([Mfile,': start selecting ',int2str(Stot),' weights. (B=',int2str(B),...
         ') VSalg=',VSalg]);
end

if B==1
   if ~strcmp(VSalg,'VSab2'); el=0; end;
   if Display
      hwbL=(el+1)*L;
      t1=['Please wait while ',int2str(hwbL),' calls to ',VSalg,' is done.'];
      if Display>1
         hwbi = 0;
         hwb = waitbar(0,t1);
      else
         disp([Mfile,': ',t1]);
      end
   end
   for i_el=0:el
      % this is the easy case
      for l=1:L
         s=S(l);                         % number of vectors for current block
         if s==1
            % should only find one weight, then find the best
            W(:,l)=zeros(K,1);
            c=(X(:,l)'*F);                % the inner products
            [temp,i]=max(abs(c));i=i(1);
            W(i,l)=c(i);
         elseif s>1
            if strcmp(VSalg,'VSab2')
               temp=full(sum(W(:,l)~=0));
               if (temp==s)
                  W(:,l)=VSab2(X(:,l),W(:,l));   % ! Vector Selection
               elseif (temp==0)
                  W(:,l)=VSfomp2(X(:,l),s);  %  Vector Selection selecting s vectors
                  % W(:,l)=VSab2(X(:,l),s);  
               else
                  disp([Mfile,': full(sum(W(:,',int2str(l),')~=0))=',int2str(temp),...
                        ' not equal S(',int2str(l),')=',int2str(s),'.']);
                  W(:,l)=VSfomp2(X(:,l),s);  %  Vector Selection selecting s vectors
                  % W(:,l)=VSab2(X(:,l),s);  
               end
            elseif strcmp(VSalg,'VSfs')
               W(:,l)=VSfs(F,X(:,l),s);            % here F is not global!!
            elseif strcmp(VSalg,'VSps2')
                % we have not given M as argument, so we pick a random value
                W(:,l)=VSps2(X(:,l),s,floor(5+10*rand(1)+1000*power(rand(1)-0.25,4)));   
               % W(:,l)=VSps2(X(:,l),s,20);   
            else
               W(:,l)=feval(VSalg,X(:,l),s);       % ! Vector Selection 
            end
         else
            W(:,l)=zeros(K,1);
         end
         S(l)=full(sum(W(:,l)~=0));
         if S(l)<s
            % if fewer than available was used, we may select more in next block
            i=l+1; if (i>L); i=1; end;
            S(i)=S(i)+s-S(l);
         end
         if Display>1
            hwbi = hwbi+1;
            waitbar(hwbi/hwbL,hwb);
         end
      end
   end
   if Display>1
      close(hwb);
   end
else    % B>1
   % we want Bf<=B to be a factor of L
   while rem(L,Bf); Bf=Bf-1; end;
   if Display
      hwbL=ceil((el+1)*L/Bf);
      t1=['Please wait while ',int2str(hwbL),' calls to ',VSalg,' is done.'];
      if Display>1
         hwbi = 0;
         hwb = waitbar(0,t1);
      else
         disp([Mfile,': ',t1]);
      end
   end
   for i_el=0:el
      for l=1:Bf:L
         ll=l+(0:(Bf-1));
         ll=mod(ll-1+i_el,L)+1;      % this line may be skipped
         % and add some more
         while length(ll)<B
            temp=floor(rand(1,1)*L+1);
            if ~ismember(temp,ll); ll=[ll,temp]; end
         end
         x=X(:,ll);x=x(:);     % BNx1
         w=W(:,ll);w=w(:);     % BKx1
         s=sum(S(ll));
         %
         if s==1
            % should only find one weight, then find the best
            w=zeros(size(w));
            c=(x'*F);                % the inner products
            [temp,i]=max(abs(c));i=i(1);
            w(i)=c(i);
         elseif s>1
            if strcmp(VSalg,'VSab2')
               if (full(sum(w~=0))==s)
                  w=VSab2(x,w);    %  Vector Selection using previous w
               else
                  w=VSfomp2(x,s);  %  Vector Selection selecting s vectors
                  % w=VSab2(r,w);    %  Vector Selection using previous w
               end
            else
               w=feval(VSalg,x,s);       % ! Vector Selection
            end
         else
            w=zeros(size(w));
         end
         %
         W(:,ll)=reshape(w,K,B);
         S(ll)=full(sum(W(:,ll)~=0));
         % it may be that fewer are selected, then 
         
         if sum(S(ll))<s
            % if fewer than available was used, we may select more in next block
            i=l+Bf+i_el; if (i>L); i=i-L; end;
            S(i)=S(i)+s-sum(S(ll));
         end
         if Display>1
            hwbi = hwbi+1;
            waitbar(hwbi/hwbL,hwb);
         end
      end
   end
   if Display>1
      close(hwb);
   end
end

% check that S is correct, this test is not relevant??
% temp=full(sum(W~=0));
% if sum(temp==S)~=L
%    disp([Mfile,': Logical error?, program did not track changes in S.']);
%    S=temp;
% end

Ssum=sum(S);
if Display
   disp([Mfile,': Number of weights selected is ',int2str(Ssum),'.']);
else
   if Ssum~=Stot
      disp([Mfile,': Number of weights used has changed from ',int2str(Stot),...
            ' to ',int2str(Ssum),'.']);
   end
end
if Ssum<Stot
   % we should try to select some more weights
   F=Fin;       % use F as if B==1
   FF=F'*F;
   I=find(S<(Stot/L));
   length(I);
   temp=floor(length(I)/(Stot-Ssum));
   t=ceil(rand(1,1)*temp);
   i=t:temp:length(I);
   ll=I(i);
   for i=1:(Stot-Ssum)
      l=ll(i);
      S(l)=S(l)+1;
      if S(l)==1
         % should only find one weight, then find the best
         W(:,l)=zeros(K,1);
         c=(X(:,l)'*F);                % the inner products
         [temp,k]=max(abs(c));k=k(1);
         W(k,l)=c(k);
      else
         W(:,l)=feval(VSalg,X(:,l),S(l));       % ! Vector Selection
      end
   end
   Ssum=sum(S);
   disp([Mfile,': Number of weights is increased to ',int2str(Ssum),'.']);
end

return