designfb.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

function ErrCode=DesignFb(FrameFile,TotIt,MaxDays)
% DesignFb  This is a stripped version of DesignF, only for block-oriented frame
%           DesignF is the complete program, and works also for block-oriented
% frames. This much more simple version also only use VSblock for vector selection.
% DesignFb is used exactly as DesignF.
% 
% ErrCode=DesignFb(FrameFile,TotIt,MaxDays);
%-----------------------------------------------------------------------------------
% arguments:
%  ErrCode   - 0 is returned if the function execute without error
%              1 an 'unexpected' error occurred, lasterr may explain it
%              2 or larger, another error occurred, 
%  FrameFile - the name of the mat-file used to store the frame
%  TotIt     - total number of iterations to do 
%              If more than TotIt iterations already has been done on this frame
%              then no iterations will be done now, the function just return and
%              do not change the frame in FrameFile.
%  MaxDays   - how many days the design process should be allowed to run
%              if MaxDays=0.25 it is 6 hours, 0.01 is 14 minutes and 24 seconds
%-----------------------------------------------------------------------------------

%----------------------------------------------------------------------
% 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.3  07.11.2001  KS: function made based on DesignF version 1.3
% Ver. 1.4  25.11.2002  KS: moved from ..\Frames to ..\FrameTools
%----------------------------------------------------------------------

Mfile='DesignFb';
ErrCode=0; 
Display=1;
UpdateFforEachSignal=0;

if (nargin<3)
   disp([Mfile,': wrong number of arguments, see help.']);
   ErrCode=2;
   return
end
StopTime=now+MaxDays;

if ~exist([FrameFile,'.mat'])
   disp([Mfile,': the FrameFile ',FrameFile,'.mat does not exist.']);
   ErrCode=3;
   return
end
try
   % load the following variables from FrameFile: Class, Type, Mdim, F, SizeF
   % G, Ctab, Dtab, Fbest, Savg, Mdat, PreProc, VecSel, InitialF, History, SNRtot
   load(FrameFile);
catch
   disp([Mfile,': the FrameFile ',FrameFile,'.mat could not be read.']);
   ErrCode=3;
   return
end
History0=char(History,[Mfile,': Training of ',FrameFile,' started ',datestr(now,0)]);
S=[];  % will be overwritten when datafile is loaded

% The data-files used for training of the frame should be ok and
% stored in the correct mat-files, DataXnnn.mat

% check and display the frame type
if Type=='b'
   disp([Mfile,': ',FrameFile,' is a block oriented frame.']);
else
   disp([Mfile,': ',FrameFile,' is type ',Type,' which is not allowed in ',Mfile,'.']);
   disp('You may try the complete program, DesignF.');
   ErrCode=5;
   return
end
if (length(VecSel.arg4)==0); VecSel.arg4=0; end;
   
% now do the main iterations
try
   K=SizeF(Mdim+1);
   N=prod(SizeF(1:Mdim));
   P=prod(SizeF((Mdim+2):(2*Mdim+1)));
   if (P>1) 
      disp([Mfile,': Type (b) and value of P (from SizeF) does not correspond.']);
      disp([Mfile,': P=',int2str(P),'>1 which is not allowed in ',Mfile,'.']);
      disp('For overlapping frames you may try the complete program, DesignF.');
      ErrCode=5;
      return
   end
   TotWWT=zeros(K,K);
   TotXWT=zeros(N,K);
   ItNo=size(SNRtot,2);
   Fnormfact=ones(K,1);
   if ItNo>0
      SNRbest=max(SNRtot(Mdat+1,:));
   else
      SNRbest=0;
   end
   if Display
      disp([Mfile,': Training started ',datestr(now,0),', ',int2str(TotIt-ItNo),' it.']);
      t3=num2str(SNRbest,4);t3=[blanks(8-length(t3)),t3];
      disp(['Best SNR is now ',t3]);
   end
   %
   for i=1:(TotIt-ItNo)
      if (now>StopTime); break; end;
      if ~UpdateFforEachSignal
         TotWWT=zeros(K,K);
         TotXWT=zeros(N,K);
      end
      ex=0;er=0;
      for m=1:Mdat
         if (Mdat>1) | (i==1) 
            % read the data file
            t1=int2str(m);t1=['0000',t1];t1=t1((length(t1)-2):length(t1));
            DataFile=['DataX',t1]; 
            if ~exist([DataFile,'.mat'])
               disp([Mfile,': the DataFile ',DataFile,'.mat does not exist.']);
               ErrCode=3;
               return
            end
            load(DataFile);    
            [n,L]=size(X);     % we need the L variable
         end
         if UpdateFforEachSignal & (i>1)
            TotWWT=TotWWT-WWT;
            TotXWT=TotXWT-XWT;
         end
         % do vector selection
         if strcmp(VecSel.Prog1,'VSblock')
            %if (i==1); S=Savg; end;
            if (length(S)==0); S=Savg; end;
            if strcmp(VecSel.arg1,'VSab2') & (i>1)
               W=(Fnormfact(:)*ones(1,L)).*W;   % W was stored instead of S
               S=W;              % and also used as third argument in VSblock
            end
            W=VSblock(X,F,S,VecSel.arg1,VecSel.arg2,VecSel.arg3,VecSel.arg4);   
            S=full(sum(W~=0));
         else
            disp([Mfile,': Try vector selection using ',VecSel.Prog1,'.']);
            W=feval(VecSel.Prog1,X,F,floor(Savg*N),VecSel.arg1);  %% OBS 16.may.2003
         end
         % calculate WWT and XWT
         WWT=full(W*W');   
         XWT=full(X*W');   
         % store DataFile with S/W and WWT and XWT
         if strcmp(VecSel.arg1,'VSab2')
            save(DataFile,'Name','X_DC','X_LP','X','SizeX','ss2','W','WWT','XWT');
         else
            save(DataFile,'Name','X_DC','X_LP','X','SizeX','ss2','S','WWT','XWT');
         end
         %
         ex=ex+ss2;    
         R=X-full(F*W);
         EnergyR=R(:)'*R(:);
         er=er+EnergyR;
         %
         TotWWT=TotWWT+WWT;
         TotXWT=TotXWT+XWT;
         %
         SNRtot(m,ItNo+i)=10*log10(ss2/EnergyR);
         if m==Mdat
            SNRtot(Mdat+1,ItNo+i)=10*log10(ex/er);
         end
         %
         if Display
            t1=int2str(i);t1=[blanks(4-length(t1)),t1];
            t2=int2str(m);t2=[blanks(4-length(t2)),t2];
            t3=num2str(SNRtot(m,ItNo+i),4);t3=[blanks(8-length(t3)),t3];
            disp(['Find W, iteration ',t1,'  signal ',t2,' (',Name,') gives SNR =',t3]);
         end
         %   
         if (UpdateFforEachSignal & i>1) | (m==Mdat)
            F=CalculateF(TotWWT,TotXWT,F,400);   % almost F=TotXWT*inv(TotWWT)
            [F,Fnormfact]=NormalizeF(F);
         end
         if (now>StopTime); break; end;
      end   % for m
      %
      if (m==Mdat)    % this is the normal situation      
         if SNRtot(Mdat+1,ItNo+i)>SNRbest
            SNRbest=SNRtot(Mdat+1,ItNo+i);
            Fbest=F;
         end
         %
         History=char(History0,[Mfile,': ',datestr(now,0),' done ',...
               int2str(i),' iterations.']);
         save(FrameFile,'Class','Type','Mdim','F','SizeF','G','Ctab','Dtab',...
            'Fbest','Savg','Mdat','PreProc','VecSel','InitialF','History','SNRtot');
         %
         if Display
            t1=int2str(i);t1=[blanks(4-length(t1)),t1];
            t3=num2str(SNRtot(Mdat+1,ItNo+i),4);t3=[blanks(8-length(t3)),t3];
            disp([FrameFile,': Iteration ',t1,'  for all signals gives SNR =',t3]);
         end
      end
      %       
   end   % for i
catch
   ErrCode=1;
end

return;

%-----------------------------------------------------------------------------------
% F=CalculateF(A,B,F,Ratio)    essentially this is: F=B*inv(A) 
% (problem if A is non-invertible is handled), usually A is ivertible  
% arguments:
%   Ratio  the ratio of largest diagonal element on to small elements in A
%          possible values could be: 400, 1000, 160*K*Savg (Savg=1/16 f.ex) 
function F=CalculateF(A,B,F,Ratio)
[N,K,P]=size(F);     
Ad=diag(A);     % all these elements are non-negative
[Adma,Ima]=max(Ad);
SmallAd=find((Ratio*Ad)<Adma);
if (length(SmallAd) | (rcond(A)<1e-6))  
   disp(['Special treatment for ',num2str(length(SmallAd)),' vectors.']);
   Ftemp=F;
   % there may be different ways to deal with singular A matrix
   if length(SmallAd)==0
      % we add values to the diagonal to avoid sigular matrix A
      disp(['Add a little bit to the diagonal of A.']);
      A=A+eye(length(Ad))*Adma/Ratio;
      F=B/A;     % and update F the usual way
   else
      % we use only the larger ones
      Ad=diag(A); 
      LargeAd=find((Ratio*Ad)>=Adma);
      disp(['Update the large ones in a normal way.']);
      F(:,LargeAd)=B(:,LargeAd)/A(LargeAd,LargeAd);     % and update F
      % the rarly used ones are changed in another way 
      if 1
         disp(['Update the small ones by adding other frame vectors.']);
         for i1=1:length(SmallAd); 
            i2=ceil(rand(1)*length(LargeAd))
            t1=0.5+rand(1);   % 0.5 < t1 < 1.5
            F(:,SmallAd(i1))=F(:,SmallAd(i1))+t1*F(:,LargeAd(i2));
         end 
      else
         disp(['Update the small ones by adding noise to the larger ones.']);
         t1=1/sqrt(10*N);
         for i1=1:length(SmallAd); 
            i2=ceil(rand(1)*length(LargeAd))
            F(:,SmallAd(i1))=F(:,LargeAd(i2))+randn(N,1)*t1;
         end 
      end
   end
else
   % this is the standard way to update F
   F=B/A;
end
%
return