designf.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

               end
            end
            clear Wt Wt1
         end
         EnergyR=R(:)'*R(:);
         er=er+EnergyR;
         %
         if Type~='s'
            TotWWT=TotWWT+WWT;
            TotXWT=TotXWT+XWT;
         end
         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)
            if Type=='g'
               % update of F must be done by using C and D matrices
               C=zeros(Q,Q);
               for i1=1:size(Ctab,1)
                  ci=Ctab(i1,1);ai=Ctab(i1,2);
                  C(ci)=C(ci)+sign(ai)*TotWWT(abs(ai));
               end
               C=C+C'-diag(diag(C));   % IMPORTANT!               
               D=zeros(Q,1);
               for i1=1:size(Dtab,1)
                  di=Dtab(i1,1);bi=Dtab(i1,2);
                  D(di)=D(di)+sign(bi)*TotXWT(abs(bi));
               end
               F=CalculateFg(C,D,G,10000);            % almost F=inv(C)*D
               [F,Fnormfact]=NormalizeF(F,G);
               % if (i==1); save('temp','C'); end;
            elseif Type=='s'
               % This is the special separable frame
               [N1,K1]=size(Fv);
               [N2,K2]=size(Fh);
               K=K1*K2;N=N1*N2;        % Fg is size NxK
               for i1=1:20
                  Fvhold=[Fv(:);Fh(:)]; 
                  % first update Fv
                  TotXWT=zeros(N1,K1);
                  TotWWT=zeros(K1,K1);
                  for m=1:Mdat
                     t1=int2str(m);t1=['0000',t1];t1=t1((length(t1)-2):length(t1));
                     DataFile=['DataX',t1]; 
                     load(DataFile); % saved in line 260, where X is NxL and W is KxL
                     [n,L]=size(X);  % we need the L variable
                     X=reshape(X,N1,N2*L);
                     W=reshape(full(W),K1,K2,L);
                     Wv=zeros(K1,N2,L);
                     for l=1:L; Wv(:,:,l)=W(:,:,l)*(Fh'); end;
                     Wv=reshape(Wv,K1,N2*L);
                     if 0  % test
                        R=X-Fv*Wv;
                        EnergyR=R(:)'*R(:);
                        disp(['Fv before iteration ',int2str(i1),' image ',Name,...
                              ' has SNR ',num2str(10*log10(ss2/EnergyR))]);
                     end
                     TotXWT=TotXWT+X*Wv';
                     TotWWT=TotWWT+Wv*Wv';
                  end
                  Fv=CalculateF(TotWWT,TotXWT,Fv,400);   % almost Fv=TotXWT*inv(TotWWT)
                  % then update Fh
                  TotXWT=zeros(N2,K2);
                  TotWWT=zeros(K2,K2);
                  for m=1:Mdat
                     t1=int2str(m);t1=['0000',t1];t1=t1((length(t1)-2):length(t1));
                     DataFile=['DataX',t1]; 
                     load(DataFile); % saved in line 260, where X is NxL and W is KxL
                     [n,L]=size(X);  % we need the L variable
                     X=reshape(X,N1,N2,L);
                     W=reshape(full(W),K1,K2,L);
                     Wh=zeros(K2,N1,L);   % actually Wh'
                     for l=1:L; Wh(:,:,l)=(Fv*W(:,:,l))'; end;
                     if 0  % test X (N1xN2xL) Wh(') (K2xN1xL) Fh (N2xK2)
                        R=X; for l=1:L; R(:,:,l)=X(:,:,l)-(Fh*Wh(:,:,l))'; end;
                        EnergyR=R(:)'*R(:);
                        disp(['Fh before iteration ',int2str(i1),' image ',Name,...
                              ' has SNR ',num2str(10*log10(ss2/EnergyR))]);
                     end
                     Wh=reshape(Wh,K2,N1*L);  
                     X=permute(X,[2,1,3]);    % change dim 1 and 2, X(:,:,l)=X(:,:,l)'
                     X=reshape(X,N2,N1*L);
                     TotXWT=TotXWT+X*Wh';
                     TotWWT=TotWWT+Wh*Wh';
                  end
                  Fh=CalculateF(TotWWT,TotXWT,Fh,400);   % almost Fh=TotXWT*inv(TotWWT)
                  change=norm(Fvhold-[Fv(:);Fh(:)]);
                  Fvhold=[Fv(:);Fh(:)];    
                  if Display
                     disp(['done Fv and Fh ',int2str(i1),...
                           ' times, change is ',num2str(change)]);
                  end
                  if (change < 1e-5); break; end;
                  if (change > 1000); 
                     disp(['done Fv and Fh ',int2str(i1),' times, change is ',...
                           num2str(change),' no convergence(?) !!']);
                     break 
                  end;
               end
               [Fh,Fhn]=NormalizeF(Fh);
               [Fv,Fvn]=NormalizeF(Fv);
               Fnormfact=kron(Fhn,Fvn);
               F{1}=Fv;
               F{2}=Fh;
               % the end for Type=='s'
            else         % Type=='b'  or Type=='o'  (if P>1)
               if (Mdim==2) & (P==4) & (L==4096)
                  % this is the special 2D overlapping case handled
                  Fstar=CalculateF(TotWWT,TotXWT,reshape(F,N,K*P),400);
                  F=reshape(Fstar,N,K,P);
                  clear Fstar
                  [F,Fnormfact]=NormalizeF(F);
               else
                  % the normal (block-oriented or 1D overlap) case
                  F=CalculateF(TotWWT,TotXWT,F,400);   % almost F=TotXWT*inv(TotWWT)
                  [F,Fnormfact]=NormalizeF(F);
               end
            end
         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);     
if P>1       % overlappping frame (and 1D signal)
   % input A is now different submatrices of the Cstar matrix
   % also B is submatrices of the D matrix (XW')
   % the current map from Cp to A is valid for 1D overlap only!!
   Cp=A;    % size should be KxKxP
   A=zeros(K*P,K*P);
   for p1=0:(P-1)
      for p2=0:(P-1)
         if p1>p2
            A((1:K)+p1*K,(1:K)+p2*K)=Cp(:,:,p1-p2+1)';
         else
            A((1:K)+p1*K,(1:K)+p2*K)=Cp(:,:,p2-p1+1);
         end
      end
   end
   B=reshape(B,N,K*P);   
   F=reshape(F,N,K*P);   
end
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))  
   if 1
      disp(['Special treatment for ',num2str(length(SmallAd)),' vectors.']);
   end
   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 
      [temp,Ima]=sort(-Ad);    % sort Ad with largest first
      if 0
         disp(['Update the small ones by adding factor of other vectors.']);
         for i1=1:length(SmallAd); 
            if SmallAd(i1)<=K
               f=(Adma+Ad(SmallAd(i1))*Ratio*0.98)/(2*Adma);  %  0.5 <= f <= 0.99
               t1=rem(Ima(i1)-1,K)+1;
               for p=0:(P-1)
                  F(:,SmallAd(i1)+p*K)=f*Ftemp(:,SmallAd(i1)+p*K)+(1-f)*Ftemp(:,t1); 
               end
            end
         end 
      else
         disp(['Update the small ones by adding noise to the largest ones.']);
         F(:,SmallAd)=F(:,Ima(1:length(SmallAd)))+randn(size(F,1),length(SmallAd))*0.01;
      end
   end
else
   % this is the standard way to update F
   F=B/A;
end
% now F is NxKP, we should resize F
if P>1; F=reshape(F,N,K,P); end;
%
return

%-----------------------------------------------------------------------------------
% F=CalculateFg(C,D,G,Ratio)    essentially this is: F=inv(C)*D
% (problem if C is non-invertible is handled), usually C is ivertible  
% arguments:
%   Ratio  the ratio of largest diagonal element on to small elements in C
%          possible values could be: 400, 1000, 160*K*Savg (Savg=1/16 f.ex) 
%   C      a QxQ matrix
%   D      a Qx1 vector
%   F      and the result is a Qx1 vector
function F=CalculateFg(C,D,G,Ratio)
d=diag(C);             % all these elements are non-negative
F=zeros(size(D));
[d_max,temp]=max(d);
Small_d=find((Ratio*d)<d_max);
if length(Small_d)
   disp(['Special treatment for ',num2str(length(Small_d)),' variables:']);
   % Small_d
   Large_d=find((Ratio*d)>=d_max);
   F(Large_d)=C(Large_d,Large_d)\D(Large_d);     % and update F
   % give som random values to the other values
   temp=mean(abs(F(Large_d)));
   if 1   % try 'low-pass'
      temp2=filter(fir1(9,0.1+0.4*rand),1,randn(length(Small_d)+10,1));
      F(Small_d)=temp*temp2(11:length(temp2));
   else   % just random
      F(Small_d)=0.5*temp*randn(length(Small_d),1);
   end   
   % save('temp','d');
else
   F=C\D;
end
%
return


% OLD code segments, not active now

% kode for VS, if strcmp(VecSel.Prog1,'BlockVS')
if VecSel.arg4
   if ~mod((i-1), VecSel.arg4)  
      [W,S]=GlobalMP(X,Fg,Savg);  
   end
else
   if (i==1)
      % perhaps this should be done for all i, by using the returned
      % S over and over again, it may become "corrupted"
      S=zeros(1,L);
      t1=0;
      for l=1:L
         t1=t1+N*Savg;
         S(l)=floor(t1);
         t1=t1-S(l);
      end
   end
end
if strcmp(VecSel.arg1,'VSab2')
   if (i==1)
      W=S;
   else
      W=(Fnormfact*ones(1,L)).*W;
   end
   [W,S]=BlockVS(X,Fg,W,VecSel.arg1,VecSel.arg2,VecSel.arg3);   
else
   % disp([Mfile,': start BlockVS, sum(S)=',int2str(sum(S)),'.']);
   [W,S]=BlockVS(X,Fg,S,VecSel.arg1,VecSel.arg2,VecSel.arg3);   
   % disp([Mfile,': ended BlockVS, sum(S)=',int2str(sum(S)),'.']);
end