vsolap1.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

function W=VSolap1(X,Fin,S,VSalg,B,el)
% VSolap1   Vector Selection for overlapping frame and 1D signal
%           a larger frame is build by repeating the input frame several
% times, i.e actual frame is input frame repeated B times. 
% Then, for each corresponding signal block
% a number of vectors are selected using a block-oriented vector selection
% algorithm, VSxxx, given by VSalg (char array).
% The program works for overlapping frames and 1D dignal.
%
% W=VSolap1(X,F,S,VSalg,B);       
% W=VSolap1(X,F,W,VSalg,B,el);       
%--------------------------------------------------------------------------------
% 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 NxKxP  (P>1)
%          the vectors of F must be normalized, [F,nf]=NormalizeF(F);
%  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.
%    1xL   number of vectors to select for each block of length N 
%          S(l) is number of non-zero weights in 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 N(B+P-1)xKB, where the input frame, 
%          reshaped into size NPxK, is used in the blocks on the diagonal.
%          We should have B larger than or equal to (P-1) to keep overlapping 
%          within the adjacent (large/expanded/actual) frame.
%  el    - extra loops to do, vector selection will be done (1+el) times 
%          for each extended block,   default is el=0
%--------------------------------------------------------------------------------
% Note: Arguments for this function is almost like for VSblock

%--------------------------------------------------------------------------------
% Copyright (c) 2001.  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
% Ver. 1.2  03.12.2002  KS: moved from ..\Frames to ..\FrameTools
%--------------------------------------------------------------------------------
% The main difference from BlockVS is that we do not have extra overlap (Be)
% instead we do vector selection twice for every block, first time when the 
% neighbor weights are zero, second time when they are set.
% This is not done when previous weights are given as input argument.
% Also the blocks may be different each time the function is called

global F FF 
Mfile='VSolap1';
Display=0;        % decide if function display information, 2 also waitbar
UseGMP=1;

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 P<2
   error([Mfile,' should have P>1, for P==1 use VSblock.']);
end
if (k==K) & (Ls==L)
   W=S;             % the old (previous) weights 
   S=full(sum(W~=0));
   [k,Ls]=size(S);
   WeightsGiven=1;
else
   W=sparse(K,L);    % the weight matrix is sparse and initial zeros
   WeightsGiven=0;
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          % the 1D overlapping GMP algorithm is below
      disp([Mfile,': Use GMP to find some weights since only Savg was given.']);
      if Display; ssx=X(:)'*X(:); disp(' '); end;
      Fp=reshape(permute(Fin,[1,3,2]),N*P,K);                 % Fp is NPxK
      R=X; for (p=2:P);  R=[R;[X(:,p:L),X(:,1:(p-1))]];  end; % R is NPxL
      Q=Fp'*R;                                                % Q is KxL
      SNL=floor(Savg*N*L);     % number of non-zeros in W left to select
      L1=(L-max([floor(Savg*L/P),1])+1);
      while SNL>0
         [temp,k]=max(abs(Q));   
         [temp,l]=sort(temp);    
         k=k(l);
         temp=ones(1,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 temp(li) 
               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)-Fp(:,ki)*Q(ki,li);
               ll=(li-P+1):(li+P-1); ll=mod(ll-1,L)+1;
               for p1=1:P
                  for p2=setdiff(1:P,p1)
                     R((1:N)+(p2-1)*N,ll(P+p1-p2))=R((1:N)+(p1-1)*N,ll(P));
                  end
               end
               temp(ll)=0;
               Q(:,ll)=Fp'*R(:,ll);
            end
         end
         if Display
            temp=R(1:N,:);
            temp=temp(:)'*temp(:);
            temp=10*log10(ssx/temp);
            disp([int2str(SNL),' weights left, SNR=',num2str(temp)]);
         end
      end   
      clear Fp R Q SNL L1 temp k l ki li p1 p2
      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 (nargout < 1); 
   error([Mfile,': function must have one output arguments, see help.']); 
end

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

if (B<(P-1))
   disp([Mfile,': (B<(P-1)) is not good.']);
end   

% now build the frame to use here from input frame, Fin
% F and FF are global variables from now on
Fin=reshape(permute(Fin,[1,3,2]),N*P,K);    % Fin is changed from NxKxP to NPxK
% then, build the frame to use here
F=zeros(N*(B+P-1),K*B);
for b=0:(B-1)
   F((1:(N*P))+b*N,(1:K)+b*K)=Fin;
end
FF=F'*F;

Stot=sum(S);  % S is now 1xL
if Display
   disp([Mfile,': size F is ',int2str(size(F,1)),'x',int2str(size(F,2)),...
         ',  Savg=',int2str(floor(Stot*B/L)),...
         ',  L=',int2str(L),',  B=',int2str(B),',  P=',int2str(P),...
         ',  N=',int2str(N),',  K=',int2str(K),',  Stot=',int2str(Stot),'.']);
   hwbL = ceil((el+1)*(L+1)/B);
   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
   disp(' ');
end

for i_el=0:el
   temp=ceil(rand(1,1)*B);         % random offset each time
   Usel=temp:B:(L+B);
   for l=Usel
      xl=l+(0:(B+P-2));  % indexes for this block (x)
      wl=l+(0:(B-1));  % indexes for this block (w)
      wl0=wl-B;        % previous block (w)
      wl1=wl+B;        % next block (w)
      if (l+2*B-1)>L
         xl=rem(xl-1,L)+1;
         wl=rem(wl-1,L)+1;
         wl0=rem(wl0-1,L)+1;
         wl1=rem(wl1-1,L)+1;
      end
      if (l-B)<1
         xl=mod(xl-1,L)+1;
         wl=mod(wl-1,L)+1;
         wl0=mod(wl0-1,L)+1;
         wl1=mod(wl1-1,L)+1;
      end
      %
      x=X(:,xl);x=x(:);     % (B+P-1)Nx1
      w=W(:,wl);w=w(:);     % BKx1
      w0=W(:,wl0);w0=w0(:);     % BKx1
      w1=W(:,wl1);w1=w1(:);     % BKx1
      s=sum(S(wl));         % number of weights available
      %
      xr0=F*w0;
      xr1=F*w1;
      r=x-[xr0((B*N+1):((B+P-1)*N));zeros(B*N,1)]-[zeros(B*N,1);xr1(1:((P-1)*N))];
      %
      if s==1
         % should only find one weight, then find the best
         w=zeros(size(w));
         c=(r'*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(r,w);    %  Vector Selection using previous w
            else
               w=VSfomp2(r,s);  %  Vector Selection selecting s vectors
               % w=VSab2(r,w);    %  Vector Selection using previous w
            end
         else
            w=feval(VSalg,r,s);       % ! Vector Selection
         end
      else
         w=zeros(size(w));
      end
      %
      W(:,wl)=reshape(w,K,B);
      S(wl)=full(sum(W(:,wl)~=0));
      if sum(S(wl))<s
         % if fewer than available was used, we may select more in next block
         i=l+B; if (i>L); i=i-L; end;
         S(i)=S(i)+s-sum(S(wl));
      end
      if Display>1
         hwbi = hwbi+1;
         waitbar(hwbi/hwbL,hwb);
      end
   end
end
if Display>1
   close(hwb);
end

% check that S is correct, is this test 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
   I=find(S<(Stot/L));
   length(I);
   temp=floor(length(I)/(Stot-Ssum));
   t=ceil(rand(1,1)*temp);
   i=t:temp:length(I);
   Usel=I(i);
   Usel=Usel(1:(Stot-Ssum));
   % this for-loop is almost exactly as above  (exceptions marked by **)
   for l=Usel
      xl=l+(0:(B+P-2));  % indexes for this block (x)
      wl=l+(0:(B-1));  % indexes for this block (w)
      wl0=wl-B;        % previous block (w)
      wl1=wl+B;        % next block (w)
      if (l+2*B-1)>L
         xl=rem(xl-1,L)+1;
         wl=rem(wl-1,L)+1;
         wl0=rem(wl0-1,L)+1;
         wl1=rem(wl1-1,L)+1;
      end
      if (l-B)<1
         xl=mod(xl-1,L)+1;
         wl=mod(wl-1,L)+1;
         wl0=mod(wl0-1,L)+1;
         wl1=mod(wl1-1,L)+1;
      end
      %
      x=X(:,xl);x=x(:);     % (B+P-1)Nx1
      w=W(:,wl);w=w(:);     % BKx1
      w0=W(:,wl0);w0=w0(:);     % BKx1
      w1=W(:,wl1);w1=w1(:);     % BKx1
      s=sum(S(wl))+1;           %                                           **
      %
      xr0=F*w0;
      xr1=F*w1;
      r=x-[xr0((B*N+1):((B+P-1)*N));zeros(B*N,1)]-[zeros(B*N,1);xr1(1:((P-1)*N))];
      %
      if s==1
         w=VSmp(r,s);
      else
         w=feval(VSalg,r,s);       % ! Vector Selection                    **
      end
      %
      W(:,wl)=reshape(w,K,B);
      S(wl)=full(sum(W(:,wl)~=0));
   end
   Ssum=sum(S);
   disp([Mfile,': Number of weights is increased to ',int2str(Ssum),'.']);
end

return