genlloyd.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

function F=GenLloyd(X,K,IsNormalized)
% GenLloyd  Design a codebook (frame) by the Generalized Lloyd Algorithm 
%           If initial codebook is not given it will be made by first assuming 
% that all vectors of X are in the same class, and the split the class with the
% largest error into two classes until K classes are made.
% IsNormalized is used to indicate if the vectors of X and F are normalized. 
% Normalized: norm(f)==1 and sum(f)>0 where f is a column of F 
% Then the codebook will be the shapes in a shape-gain quantizer
%
% F=GenLloyd(X,K,IsNormalized);        % K is a number
% F=GenLloyd(X,F0,IsNormalized);       % F0 is a matrix of size NxK
%----------------------------------------------------------------------
% Arguments:
%  F            - the codebook vectors, size is NxK
%  X            - the training vectors, a matrix of size NxL
%  K            - number of vectors in codebook
%  F0           - the initial codebook vectors, size is NxK
%  IsNormalized - could be 1 or 0, default is 0.
%----------------------------------------------------------------------

%----------------------------------------------------------------------
% 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  25.12.2002  KS: function made based on ..\Frames\GenLloyd.m
%----------------------------------------------------------------------

Mfile='GenLloyd';

if (nargin==2); IsNormalized=0; end;
if ((nargin==2) | (nargin==3))
    [N,L]=size(X);
    if prod(size(K))>1
        F=K;
        [n,K]=size(F);
        if (n~=N)
            error([Mfile,': size of X and F0 do not match, see help.']);
        end
        DoInitialization=0;
    else
        F=zeros(N,K);
        DoInitialization=1;
    end
else
   error([Mfile,': wrong number of arguments, see help.']);
end

BelongToClass=zeros(1,L);
NormOfErrorSq=zeros(1,L);
VectorsInClass=zeros(K,1);
ErrorsInClass=zeros(K,1);
sumX2=sum(sum(X.*X));
if DoInitialization
    f1=mean(X,2);
    if IsNormalized; f1=f1/sqrt(f1'*f1); end;
    BelongToClass=ones(1,L);
    NormOfErrorSq=ones(1,L);
    for j=1:L
        e=X(:,j)-f1;
        NormOfErrorSq(j)=e'*e;
    end
    VectorsInClass(1)=L;
    ErrorsInClass(1)=sum(NormOfErrorSq);
    k=1;   % number of classes now
    F(:,1)=f1;
    while k<K   % split one of the classes
        if 0      % split largest class
            [kmax,classmax]=max(VectorsInClass);
        else      % split class with largest errors into two
            [kmax,classmax]=max(ErrorsInClass);
        end
        if length(classmax)>1; classmax=classmax(1); end;
        I=find(BelongToClass==classmax);
        % find the most distant vector in this class
        Fd=X(:,I)-F(:,classmax)*ones(1,length(I));
        for i=1:length(I); Fd(1,i)=Fd(:,i)'*Fd(:,i); end;
        Fd=Fd(1,:);
        [Fdmax,imax]=max(Fd);
        k=k+1;    % a new class is found
        f1=F(:,classmax)+0.1*X(:,I(imax));
        if IsNormalized; f1=f1/sqrt(f1'*f1); end;
        F(:,k)=f1;
        f1=F(:,classmax)-0.1*X(:,I(imax));
        if IsNormalized; f1=f1/sqrt(f1'*f1); end;
        F(:,classmax)=f1;
        for GLAi=1:5    % now do five Lloyd iterations 
            % find class for each vector
            if IsNormalized
                [W,BelongToClass]=max(abs(F(:,1:k)'*X));
                NormOfErrorSq=1-W.*W;    % correct only when ||f||==1 and ||x||==1
            else
                for j=1:L
                    x=X(:,j);
                    NormOfErrorSq(j)=sumX2;
                    for kk=1:k
                        d=x-F(:,kk);
                        e2=d'*d;
                        if (e2<NormOfErrorSq(j))
                            NormOfErrorSq(j)=e2;
                            BelongToClass(j)=kk;
                        end
                    end
                end
            end
            for i=1:k     % update centroid in each class
                I=find(BelongToClass==i);
                VectorsInClass(i)=length(I);
                ErrorsInClass(i)=sum(NormOfErrorSq(I));
                if length(I)>1
                    f1=mean(X(:,I),2);
                elseif length(I)==1
                    f1=X(:,I);
                end
                if IsNormalized; f1=f1/sqrt(f1'*f1); end;
                F(:,i)=f1;
            end
        end
        disp([Mfile,': ',int2str(k),' vectors initialized in F.']);
    end
end

% start the GLA iterations, initial codebook vectors given in F
ErrorsIt=[];
F=F(1:N,1:K);
for GLAi=1:30   % many times
    % find class for each vector
    if IsNormalized
        [W,BelongToClass]=max(abs(F(:,1:k)'*X));
        NormOfErrorSq=1-W.*W;    % correct only when ||f||==1 and ||x||==1
    else
        for j=1:L
            x=X(:,j);
            NormOfErrorSq(j)=sumX2;
            for kk=1:K
                d=x-F(:,kk);
                e2=d'*d;
                if (e2<NormOfErrorSq(j))
                    NormOfErrorSq(j)=e2;
                    BelongToClass(j)=kk;
                end
            end
        end
    end
    ErrorsIt=[ErrorsIt,sum(NormOfErrorSq)];
    % SNR is for the (often wrongly) assumption that mean X is zero. 
    disp([Mfile,': classification in iteration ',int2str(GLAi),...
            ' gives total error ',num2str(ErrorsIt(GLAi),4),...
            ' (SNR=',num2str(10*log10(sumX2/ErrorsIt(GLAi)),4),')']);
    for i=1:K     % update centroid in each class
        I=find(BelongToClass==i);
        VectorsInClass(i)=length(I);
        ErrorsInClass(i)=sum(NormOfErrorSq(I));
        if length(I)>1
            f1=mean(X(:,I),2);
        elseif length(I)==1
            f1=X(:,I);
        end
        if IsNormalized; f1=f1/sqrt(f1'*f1); end;
        F(:,i)=f1;
    end
    if GLAi>5    % check if convergence is reached
        a=(ErrorsIt(GLAi-1)-ErrorsIt(GLAi))/ErrorsIt(GLAi);
        if a<1e-5; break; end;
    end
end

return