vsfs.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

function w=VSfs(F,x,S,maxM)
% VSfs      Vector Selection algorithm doing Full Search
%
% w=VSfs(F,x,S);
%-----------------------------------------------------------------------------------
% arguments:
%   F    - the normalized dictionary (or F matrix), size NxK 
%   x    - the vector to approximate, size Nx1
%   S    - number of vectors to select or non-zero weights, 1<=S<=N 
%   w    - the weights where at most S of the elements are non-zero, size Kx1
%-----------------------------------------------------------------------------------
% Note that F is not global in this function as in many other VSxxx functions


%----------------------------------------------------------------------
% 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:
% Ver. 1.0  30.11.2001  KS: function made 
% Ver. 1.1  04.12.2002  KS: moved from ..\Frames\ to ..\FrameTools
%----------------------------------------------------------------------

Mfile='VSfs';

if nargin<3
   error([Mfile,': wrong number of arguments, see help.']);
end
if nargin<4
   maxM=500000;   % maximal number of conbinations to consider
end

[N,K]=size(F);
x=x(:);
if length(x)~=N
   error([Mfile,': size of x and F do not match, see help.']);
end

S=floor(S(1));
if S<1; S=1; end;
if S>N; S=N; end;

M=nchoosek(K,S);      % the number of combinations to do
if M>maxM
   disp([Mfile,': number of combinations to do is M=',int2str(M),...
         ', and that is too large.']);
   w=F(:,1:S)\x;    
   return
end
if S>16
   disp([Mfile,': number of vectors to select is S=',int2str(S),...
         ', and that is too large.']);
   w=F(:,1:S)\x;    
   return
end

% improvements should be larger than a small number to count!
xx=x'*x;
xeps=2*eps*xx;   % a small number
rrbest=xx;

% start the search
wwbest=inf;
II=nchoosek(1:K,S);      % this is a large matrix of indexes!
for m=1:M
   Fm=F(:,II(m,:));
   if rank(Fm)==S
      wm=Fm\x;               % only one soultion, find this
   else      % what now, rank(Fm)<S
      % wm=VSfs(Fm,x,S-1);   % select fewer non-zero weights or
      wm=pinv(Fm)*x;         % use pseudoinverse to get minimum norm solution
   end
   rm=x-Fm*wm;
   rr=rm'*rm;
   if (rr+eps)<rrbest    % large enough improvement to count
      rrbest=rr;
      w=zeros(K,1);w(II(m,:))=wm;
      wwbest=w'*w;
   elseif (rr-eps)<rrbest   % approximate equal error
      if (wm'*wm)<wwbest      % check norm of weights
         w=zeros(K,1);w(II(m,:))=wm;
         wwbest=w'*w;
      end
   end
end

return;


% test of function
clear all
N=8;K=16;
% want a uniform distribution within the unit ball 
F=zeros(N,K);
k=1;
rand('state',1000);
while k <= K-2
   x=2*rand(N,1)-1;    % uniform within unit cube (-1, 1)
   if x'*x<1           % is vector within unit ball ?
      F(:,k)=x;        % then use it
      k=k+1;
   end
end
F(:,k)=F(:,1)+F(:,2);        % some are linear dependent of only some few
k=k+1;F(:,k)=F(:,3)+F(:,4)+F(:,5);
F=NormalizeF(F);
%
w=zeros(K,1);
w(1)=0.8;w(2)=0.8;w(4)=0.5;w(15)=0.5;
x=F*w;
%
tic;
w1=VSfs(F,x,1);
w2=VSfs(F,x,2);
w3=VSfs(F,x,3);
w4=VSfs(F,x,4);
w5=VSfs(F,x,5);
toc;
%
disp(['Length (norm) of the signal is   ',num2str(norm(x))]);
disp(['Error when selecting 1 vector is ',num2str(norm(x-F*w1))]);
disp(['Error when selecting 2 vector is ',num2str(norm(x-F*w2))]);
disp(['Error when selecting 3 vector is ',num2str(norm(x-F*w3))]);
disp(['Error when selecting 4 vector is ',num2str(norm(x-F*w4))]);
disp(['Error when selecting 5 vector is ',num2str(norm(x-F*w5))]);
disp('The weights:');
disp('         w        w1        w2        w3        w4       w5');
disp([w,w1,w2,w3,w4,w5]);
disp('And the norm of the weights.');
disp([sqrt(w'*w),sqrt(w1'*w1),sqrt(w2'*w2),sqrt(w3'*w3),sqrt(w4'*w4),sqrt(w5'*w5)]);
%
% even if 4 non zero weights are used, since they are linearly dependent
% only 3 vectors are needed for the perfect recostruction.
% But using more may reduce the norm of w.
%