frameproperties.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

    end
    % The minimum lower frame bound for all combinations of 3 vectors
    I=intersect(out,23);
    if length(I) 
        I=nchoosek(1:K,3);      % this is a large matrix of indexes!
        A3=1;
        for i=1:length(I)
            Fi=F(:,I(i,:));
            e=eig(Fi'*Fi);
            J=find(e>1e-14);
            A3=min([A3;e(J)]);
        end
        Prop{23}=A3;
    end
    % The minimum lower frame bound for all combinations of 4 vectors
    I=intersect(out,24);
    if length(I) 
        I=nchoosek(1:K,4);      % this is a large matrix of indexes!
        A4=1;
        for i=1:length(I)
            Fi=F(:,I(i,:));
            e=eig(Fi'*Fi);
            J=find(e>1e-14);
            A4=min([A4;e(J)]);
        end
        Prop{23}=A4;
    end
else                       % overlapping frame, P>1
    % The following properties are based on the eigenvalues of the frame operator
    I=intersect(out,[1:6,9]);
    if length(I)   
        j=sqrt(-1);
        theta=0:(1/ThetaLength):(1-1/ThetaLength);   
        A=zeros(1,ThetaLength);
        B=zeros(1,ThetaLength);
        lambda=zeros(N,ThetaLength);
        for i=1:ThetaLength;
            z=exp(j*2*pi*theta(i));
            R=F(:,:,1);
            for p=2:P; R=R+F(:,:,p)*z^(1-p); end;
            S=R*R';             % generally S is not real
            e=real(eig(S));     % but the eigenvalues are real
            [temp,I]=sort(-e); e=e(I);  % sort eigenvalue with largest first
            lambda(:,i)=e;
            I=find(e>1e-12);e=e(I);     % non-zero values only
            A(i)=min(e);
            B(i)=max(e);
        end
        Prop{1}=min(A);
        Prop{2}=max(B);
        Prop{3}=lambda(:,1);
        Prop{4}=A;
        Prop{5}=B;
        Prop{6}=lambda;
        % the betamse property is calculated from eigenvalues!
        % but this is a difficult measure, what should be used
        %    % this may give values out of range
        %    betamse=K/(K-1)-sum(sum(lambda.*lambda))/(K*(K-1)*ThetaLength);
        %    % also this is perhaps difficult ???
        %    temp=zeros(1,ThetaLength);
        %    for i=1:ThetaLength;
        %        temp(i)=(lambda(:,i)'*lambda(:,i))/sum(lambda(:,i));
        %    end
        %    betamse=(K-sum(temp)/ThetaLength)/(K-1);     % sine squared
        %    % this is perhaps the better variant 
        temp=zeros(1,ThetaLength);
        for i=1:ThetaLength;
            temp(i)=(lambda(:,i)'*lambda(:,i))/(sum(lambda(:,i)))^2;
        end
        betamse=K*(1-sum(temp)/ThetaLength)/(K-1);     % sine squared
        %
        betamse=asin(sqrt(betamse))*180/pi;  % degrees
        Prop{9}=betamse;
    end
    % The following properties are for angles between frame vectors 
    I=intersect(out,7:8);
    if length(I)   
        FF2=zeros(N*P,K*(2*P-1));
        for p1=1:P
            k=(p1-2+p)*K+(1:K);
            for p=1:P
                n=(p1+p-2);
                if n>(P-1); n=n-P; k2=k-K*P; else k2=k; end;
                n=n*N+(1:N);
                FF2(n,k2)=F(:,:,p);
            end
        end
        F1=reshape(permute(F,[1,3,2]),N*P,K);
        C=F1'*FF2;
        C(:,K*(P-1)+(1:K))=C(:,K*(P-1)+(1:K))-eye(K);
        beta=max(abs(C),[],2);    % 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;        
    end
    % find some angles
    I=intersect(out,[29,32]);
    if length(I)   
        Fhat=zeros(N*(2*P-1),K*P);
        for p1=1:P
            for p2=1:P
                Fhat(((p1+p2-2)*N)+(1:N),((p2-1)*K)+(1:K))=F(:,:,p1);
            end
        end
        cA=Fhat'*Fhat;   % cosine of angles 
        dA=diag(cA);
        if (max(dA)>1) | (min(dA)<1)
            dA=sqrt(dA);
            cA=cA.*(dA*dA');
            I=find(cA>1); cA(I)=1;
            I=find(cA<(-1)); cA(I)=(-1);
        end
        Ang=acos(cA)*180/pi;    % degrees
        Ang=triu(Ang);          % upper right triangular matrix
        for k=1:(P*K)
            Ang(k,k)=0;
        end
        Prop{29}=Ang;
        clear cA dA 
    end
end

% the figures are made
%             31 - Plot frame vectors in figure 1
%             32 - Plot frame angles in figure 2
%             33 - Plot frequency response of frame vectors in figure 3
if length(find(out==31)) 
    figure(1);clf;
    if K>32
        PlotF(F,1,0,0,2);
    elseif K>16
        PlotF(F,2,0,0,3);
    else
        PlotF(F,0,0,0,1);
    end
end
if length(find(out==32)) 
    figure(2);clf;
    hold on;
    A1=triu((Ang<5)+(Ang>175));
    A2=((Ang>5).*(Ang<15)+(Ang<175).*(Ang>165));
    A3=((Ang>15).*(Ang<30)+(Ang<165).*(Ang>150));
    A4=((Ang>30).*(Ang<45)+(Ang<150).*(Ang>135));
    A5=((Ang>45).*(Ang<60)+(Ang<135).*(Ang>120));
    A6=((Ang>60).*(Ang<88)+(Ang<120).*(Ang>92));
    A7=((Ang>88).*(Ang<92));
    Leg='First line';
    [I,J]=find(A1);plot(I,J,'or');
    if length(I); Leg=char(Leg,'0-5 and 175-180 degrees'); end;
    [I,J]=find(A2);plot(I,J,'*r');
    if length(I); Leg=char(Leg,'5-15 and 165-175 degrees'); end;
    [I,J]=find(A3);plot(I,J,'^b');
    if length(I); Leg=char(Leg,'15-30 and 150-165 degrees'); end;
    [I,J]=find(A4);plot(I,J,'sb');
    if length(I); Leg=char(Leg,'30-45 and 135-150 degrees'); end;
    [I,J]=find(A5);plot(I,J,'hb');
    if length(I); Leg=char(Leg,'45-60 and 120-135 degrees'); end;
    [I,J]=find(A6);plot(I,J,'.g');
    if length(I); Leg=char(Leg,'60-88 and 92-120 degrees'); end;
    [I,J]=find(A7);plot(I,J,'+g');
    if length(I); Leg=char(Leg,'88-92 degrees'); end;
    Leg=Leg(2:size(Leg,1),:);
    if ((K*P)<20)
        I=find(Ang>0.001);
        T=int2str(round(Ang(I)));
        H=text(rem(I-1,K*P)+1.2,floor((I-1)/(K*P))+1,T);
        set(H,'FontSize',8);
        H=text(4,2,'Angle between two vectors (degrees) is printed in plot.');
        set(H,'FontSize',8);
    else
        legend(Leg,1);
    end
    axis([0,K*P+1,0,K*P+1]);
    hold off;
    set(gca,'ydir','reverse');
    title('Image of angles between vectors in F.');
    xlabel(['Vector number (up to KP=',int2str(K),'x',int2str(P),')']);
    ylabel(['Vector number (up to KP=',int2str(K),'x',int2str(P),')']);
end
if length(find(out==33)) 
    figure(3);clf;
    if P>1
        Fhat=zeros(N*(2*P-1),K*P);
        for p1=1:P
            for p2=1:P
                Fhat(((p1+p2-2)*N)+(1:N),((p2-1)*K)+(1:K))=F(:,:,p1);
            end
        end
    else
        Fhat=F;
    end
    Ff=fft(Fhat(:,1:K),256);
    Ff=Ff/(max(abs(Ff(1,:))));
    Ff=Ff(1:129,:);
    w=linspace(0,1,129);
    Ff=20*log10(abs(Ff)+1e-10);
    plot(w,Ff);
    V=[0,1,-40,1];axis(V);
    grid on;
    xlabel('Normalized frequency, pi=1.');
    ylabel('Magnitude in desibel');
    title('Frequency response for vectors in Fhat.');
end

t=[Mfile,' returns: '];
for i=1:NoOut
    varargout{i}=Prop{out(i)}; 
    switch out(i)
    case 1 
        t=[t,'A, '];
        if Display; disp(['Lower frame bound, A=',num2str(Prop{out(i)})]); end;
    case 2 
        t=[t,'B, '];
        if Display; disp(['Upper frame bound, B=',num2str(Prop{out(i)})]); end;
    case 3 
        t=[t,'lambda, '];
    case 4 
        t=[t,'A(th), '];
    case 5 
        t=[t,'B(th), '];
    case 6 
        t=[t,'lambda(th), '];
    case 7 
        t=[t,'betamin, '];
        if Display; disp(['betamin=',num2str(Prop{out(i)})]); end;
    case 8 
        t=[t,'betaavg, '];
        if Display; disp(['betaavg=',num2str(Prop{out(i)})]); end;
    case 9
        t=[t,'betamse, '];
        if Display; disp(['betamse=',num2str(Prop{out(i)})]); end;
    case 10
        t=[t,'r1, '];
    case 11
        t=[t,'r1max, '];
        if Display; disp(['r1max=',num2str(Prop{out(i)})]); end;
    case 12
        t=[t,'r1avg, '];
        if Display; disp(['r1avg=',num2str(Prop{out(i)})]); end;
    case 13
        t=[t,'r1mse, '];
        if Display; disp(['r1mse=',num2str(Prop{out(i)})]); end;
    case 14
        t=[t,'r2, '];
    case 15
        t=[t,'r2max, '];
        if Display; disp(['r2max=',num2str(Prop{out(i)})]); end;
    case 16
        t=[t,'r2avg, '];
        if Display; disp(['r2avg=',num2str(Prop{out(i)})]); end;
    case 17
        t=[t,'r2mse, '];
        if Display; disp(['r2mse=',num2str(Prop{out(i)})]); end;
    case 18
        t=[t,'r3, '];
    case 19
        t=[t,'r3max, '];
        if Display; disp(['r3max=',num2str(Prop{out(i)})]); end;
    case 20
        t=[t,'r3avg, '];
        if Display; disp(['r3avg=',num2str(Prop{out(i)})]); end;
    case 21
        t=[t,'r3mse, '];
        if Display; disp(['r3mse=',num2str(Prop{out(i)})]); end;
    case 22
        t=[t,'A2, '];
        if Display; disp(['A2=',num2str(Prop{out(i)})]); end;
    case 23
        t=[t,'A3, '];
        if Display; disp(['A3=',num2str(Prop{out(i)})]); end;
    case 24
        t=[t,'A4, '];
        if Display; disp(['A4=',num2str(Prop{out(i)})]); end;
    case 25
        t=[t,'fA, '];
    case 26
        t=[t,'fB, '];
    case 27
        t=[t,'betaavg2, '];
        if Display; disp(['betaavg2=',num2str(Prop{out(i)})]); end;
    case 28
        t=[t,'betaavg3, '];
        if Display; disp(['betaavg3=',num2str(Prop{out(i)})]); end;
    case 29
        t=[t,'Ang, '];
    otherwise 
        t=[t,'?, '];
    end
end
t=t(1:(end-2));
disp(t);

return