frameproperties.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

function varargout=FrameProperties(F,out)
% FrameProperties   Find, display and return some properties of a frame F
%                   Many frame properties are only defined for
% block oriented frames. The number of out arguments must correspond
% to the number of entries in 'outargno' (that return output arguments).
% For more information on the different frame properties you should read 
% chapter 7 in the thesis "Sparse Signal Representation using Overlapping 
% Frames" by Karl Skretting, available at: http://www.ux.his.no/~karlsk/
%
% outarg=FrameProperties(F,out);    % the general way of calling this function
% outarg=FrameProperties(FrameFile,out);   
% A=FrameProperties(F,1);            % the lower frame bound A
% B=FrameProperties(F,2);            % the upper frame bound B
% [A,B]=FrameProperties('FrameEx1s20',[1,2,31,32,33]);  % more examples
% [A,B,betamin,r1max]=FrameProperties('FrameEx1s20',[1,2,7,11,31,32,33]);  
% [lambda,r2max,r3max,betaavg,betamse,betaavg2]=FrameProperties('FrameEx1s20',[3,15,19,8,9,27]);  
% [a,b,c,d,e,f,g,Ang]=FrameProperties('FrameEx1s20',[2,3,7,8,9,27,28,29,31,32,33]);  
% [a,b]=FrameProperties('FrameEx1s20',[11,15]);  
%-----------------------------------------------------------------------------------
% arguments:
%  F         - The frame, size is NxKxP 
%  FrameFile - the name of the mat-file used to store the frame
%  out       - An array telling which properties to return
%            The first six values are based on the eigenvalues of the frame operator
%              1 - A, lower frame bound 
%              2 - B, upper frame bound 
%              3 - lambda, eigenvalues of frame operator
%             25 - fA, eigenvector corresponding to smallest eigenvalue
%             26 - fB, eigenvector corresponding to largest eigenvalue
%                  For overlapping frames, P>1, the following properties may be relevant
%              4 - Atheta, lower frame bounds for different values of theta
%              5 - Btheta, upper frame bounds, theta= 0:0.002:0.998  (when length=500) 
%              6 - lambdatheta, eigenvalues of frame operator, size is (N x length)
%            The following properties are for angles (in degrees) between frame vectors 
%              7 - betamin, the smallest angle between two frame vectors, 
%              8 - betaavg, average for all frame vectors to closest neighbor 
%              9 - betamse, average of all angles in mean square error sense 
%             27 - betaavg2, average of angles between fB and the frame vectors 
%             28 - betaavg3, average of angles between f-mean and the frame vectors 
%             29 - Ang,  the angle in degrees between two vectors in F.
%            The representation error when selectiong 1, 2 or 3 frame vectors
%            Now we use L=10000 gaussian test vectors and estimate the properties
%            These properties are only returned for block-oriented frames, P=1
%             10 - r1, returns all the L errors r1(l)= || x(l)-F*w(l) ||
%             11 - r1max, estimate for max error 
%             12 - r1avg, estimate for mean (average) error
%             13 - r1mse, estimate for square root of the mean square error
%             14,15,16,17  - r2, r2max, r2avg and r2mse
%             18,19,20,21  - r3, r3max, r3avg and r3mse
%            The minimum lower frame bound for all combinations of s vectors
%            These properties are only returned for block-oriented frames, P=1
%             22,23,24 - A2, A3 and A4
%            The following numbers are allowed in the end of outargno, they
%            do not give an output argument but produce a figure
%             31 - Plot frame vectors in figure 1
%             32 - Plot frame angles in figure 2
%             33 - Plot frequency response of frame vectors in figure 3
%  outarg    - the respective output arguments
%-----------------------------------------------------------------------------------

%-----------------------------------------------------------------------------------
% Copyright (c) 2002.  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  19.03.2002  KS  made function
% Ver. 1.1  10.04.2002  KS  added some of the eigenvectors, fA, fB, and betaav2
% Ver. 1.2  26.04.2002  KS  property 13, 17 and 21 changed from: estimate for the mean square error
% Ver. 1.3  30.04.2002  KS  added betaav3 property
% Ver. 1.4  30.04.2002  KS  swapped the (definitions for) betaav2 and betaav3 properties
% Ver. 1.5  21.06.2002  KS  swapped betaavg2 and betaavg3 properties again
% Ver. 1.6  04.12.2002  KS: moved from ..\Frames\ to ..\FrameTools, and added figures
%-----------------------------------------------------------------------------------

Mfile='FrameProperties';
TotalProperties=29;
ThetaLength=500;      % length of theta
L=10000;              % number of (gaussian) test vectors
Display=1;

if (nargin~=2)
    error([Mfile,' should have 2 input arguments, see help.']);
end
t=[Mfile,' using given frame'];
if ischar(F)
    FrameFile=F;
    if ~exist([FrameFile,'.mat'])
        disp([Mfile,': can not open FrameFile.']);
        return;
    end
    % load the following variables from FrameFile: Class, Type, Mdim, F, SizeF
    % G, Ctab, Dtab, Fbest, Savg, Mdat, PreProc, VecSel, InitialF, History, SNRtot
    load(FrameFile);
    if (Type == 'g')
        F=BuildFg(F,G);
    end
    clear Class Type Mdim SizeF G Ctab Dtab Fbest Savg Mdat PreProc VecSel InitialF History SNRtot
    t=[Mfile,' using frame in ',FrameFile];
end

[N,K,P]=size(F);
NoOut=sum(out<=TotalProperties);
if nargout~=NoOut
    error([Mfile,' number of output arguments is wrong.']);
end
t=[t,' of size ',int2str(N),'x',int2str(K),'x',int2str(P),'.'];
if Display; disp(t); end;

Prop=cell(TotalProperties,1);

if P==1
    % The following properties are based on the eigenvalues of the frame operator
    I=intersect(out,[1:6,9,25,26,27]);
    if length(I)   
        S=F*F';            % the frame operator
        [V,lambda]=eig(S);      % eigendecomposition of S
        lambda=diag(lambda);         
        [temp,I]=sort(-lambda); lambda=lambda(I);  % sort eigenvalue with largest first
        V=V(:,I);                                  % sort eigenvectors with largest first
        fB=V(:,1);
        Prop{25}=V(:,end);
        Prop{26}=fB;
        Prop{3}=lambda;
        Prop{6}=lambda;
        I=find(lambda>1e-12);lambda=lambda(I);     % non-zero values only
        Prop{1}=min(lambda);
        Prop{2}=max(lambda);
        Prop{4}=min(lambda);
        Prop{5}=max(lambda);
        % the betamse property is calculated from eigenvalues!
        betamse=(K*K-lambda'*lambda)/(K*(K-1));     % sine squared
        betamse=asin(sqrt(betamse))*180/pi;         % angle in degrees
        Prop{9}=betamse;
    end
    % Find betaavg2 
    I=intersect(out,27);
    if length(I)   
        C=F'*fB;              % cosine of angles
        beta=acos(abs(C));     % angles in radians
        betaavg2=mean(beta)*180/pi;  % degrees
        Prop{27}=betaavg2;
    end
    % Find betaavg3
    I=intersect(out,28);
    if length(I)   
        Ftemp=F;
        while length(I)
            fm=sum(Ftemp,2);
            t=fm'*fm; if t>0; fm=fm/sqrt(t); end;
            C=Ftemp'*fm;              % cosine of angles
            I=find(C<0);
            for i=I; Ftemp(:,i)=-Ftemp(:,i); end;
        end
        beta=acos(C);     % angles in radians
        betaavg3=mean(beta)*180/pi;   % degrees
        Prop{28}=betaavg3;
    end
    % The following properties are for angles between frame vectors 
    I=intersect(out,[7,8,29,32]);
    if length(I)   
        C=F'*F;              % cosine of angles
        dA=diag(C);
        if (max(dA)>1) | (min(dA)<1)   % if frame vectors were not normalized (uniform frame)
            dA=sqrt(dA);
            C=C.*(dA*dA');
            J=find(C>1); C(J)=1;
            J=find(C<(-1)); C(J)=(-1);
        end
        Ang=acos(C)*180/pi;
        Ang=triu(Ang);          % upper right triangular matrix
        for k=1:K; Ang(k,k)=0; end;
        Prop{29}=Ang;    
        C=C-eye(K);          % remove the diagonal
        beta=max(abs(C));    % largest cosine is smallest angle, size 1xK
        beta=acos(beta);     % angles in radians
        betamin=min(beta);
        betaavg=mean(beta);
        Prop{7}=betamin*180/pi;  % angles in degrees
        Prop{8}=betaavg*180/pi;
        clear dA 
    end
    % The following properties are representation error selectiong 1 frame vector
    I=intersect(out,10:13);
    if length(I)   
        X=randn(N,L);
        X=X./(ones(N,1)*sqrt(sum(X.*X)));     % normalize each column vector
        r1=MCsimF(X,F,1);   % L random on the unit ball   
        [temp,I]=sort(-r1);
        X=X(:,I(1:5));      % the most difficult vectors
        % use optimization to find the worst case;  r1max=MCsimF(X,F,1)   
        warning off;
        for i=1:ceil(5-N/4)
            x=X(:,i);
            try
                x=fminunc('FramePropertiesObjFun',x,[],F,1);
            catch
                x=X(:,i);
            end
            X(:,i)=x;
        end
        warning on;
        X=NormalizeF(X);
        r1max=MCsimF(X,F,1);
        r1max=max(r1max);
        r1avg=mean(r1);
        r1mse=sqrt(mean(r1.*r1));
        Prop{10}=r1;
        Prop{11}=r1max;
        Prop{12}=r1avg;
        Prop{13}=r1mse;
    end
    % The following properties are representation error selectiong 2 frame vectors
    I=intersect(out,14:17);
    if length(I) 
        if N>2
            X=randn(N,L);
            X=X./(ones(N,1)*sqrt(sum(X.*X)));     % normalize each column vector
            r2=MCsimF(X,F,2);   % L random on the unit ball   
            [temp,I]=sort(-r2);
            X=X(:,I(1:5));      % the most difficult vectors
            % use optimization to find the worst case; 
            warning off;
            for i=1:ceil(5-N/4)
                x=X(:,i);
                try
                    x=fminunc('FramePropertiesObjFun',x,[],F,2);
                catch
                    x=X(:,i);
                end
                X(:,i)=x;
            end
            warning on;
            X=NormalizeF(X);
            r2max=MCsimF(X,F,2);
            r2max=max(r2max);
            r2avg=mean(r2);
            r2mse=sqrt(mean(r2.*r2));
            Prop{14}=r2;
            Prop{15}=r2max;
            Prop{16}=r2avg;
            Prop{17}=r2mse;
        else
            Prop{14}=0;
            Prop{15}=0;
            Prop{16}=0;
            Prop{17}=0;
        end
    end
    % The following properties are representation error selectiong 3 frame vectors
    I=intersect(out,18:21);
    if length(I) 
        if N>3
            X=randn(N,L);
            X=X./(ones(N,1)*sqrt(sum(X.*X)));     % normalize each column vector
            r3=MCsimF(X,F,3);   % L random on the unit ball   
            [temp,I]=sort(-r3);
            X=X(:,I(1:5));      % the most difficult vectors
            % use optimization to find the worst case; 
            warning off;
            for i=1:ceil(5-N/4)
                x=X(:,i);
                try
                    x=fminunc('FramePropertiesObjFun',x,[],F,3);
                catch
                    x=X(:,i);
                end
                X(:,i)=x;
            end
            warning on;
            X=NormalizeF(X);
            r3max=MCsimF(X,F,3);
            r3max=max(r3max);
            r3avg=mean(r3);
            r3mse=sqrt(mean(r3.*r3));
            Prop{18}=r3;
            Prop{19}=r3max;
            Prop{20}=r3avg;
            Prop{21}=r3mse;
        else
            Prop{18}=0;
            Prop{19}=0;
            Prop{20}=0;
            Prop{21}=0;
        end
    end
    % The minimum lower frame bound for all combinations of 2 vectors
    I=intersect(out,22);
    if length(I) 
        I=nchoosek(1:K,2);      % this is a large matrix of indexes!
        A2=1;
        for i=1:length(I)
            Fi=F(:,I(i,:));
            e=eig(Fi'*Fi);
            J=find(e>1e-14);
            A2=min([A2;e(J)]);
        end
        Prop{22}=A2;