function [K,erg,vK] = mdl(e,N,gamma,verbose);%MDL - find number of principal components such that it is optimal%      with respect to the minimum description length%      first in the input come the expected signal components %% [K,erg,vK] = mdl(e,N,gamma,verbose);%%  e       - signal strengths (row vector)%  N       - number of samples%  gamma   - gamma in MDL estimator (>0 mdl(vetter), =0 aic, <0 mdl)%%  erg     - array of mdl values%  K       - argmin(erg(i))%  vK      - erg(K)%% 1.8.03 - initial version% 1.9.04 - rewrite% Peter Grubererror(nargchk(2,4,nargin))     % check no. of input argsif (nargin<3); gamma=32; endif (nargin<4); verbose=0; end% check argumentsif ( (length(N)>1) || (length(gamma)>1) ) error('Wrong arguments to MDL!\n');endif (gamma>0); mdltype='p';elseif (gamma==0); mdltype='s'; else mdltype='r';end;gamma=abs(gamma);m=length(e);% in zeilenvektor umwandelnif (size(e,2)==1) e=e'; end% normalisiere the signal strengthse(e<0)=0;for i=1:size(e,1) e(i,:)=real(e(i,:)/norm(e(i,:)));end% trivial casesif (m==1)      K=1; erg=[K];       return; endif (sum(e)==0) K=m; erg=ones(m)*K; return; endmaxi=1; while ((maxi<m) && (all(e(:,maxi+1))>0)); maxi=maxi+1; ende=e(1:min(maxi,size(e,1)),1:maxi);% initialize mdlerg=ones(3,maxi)*NaN;% get mdli=[1:maxi];M=(i*maxi)-((i.*(i-1))/2)+1;t=[1:maxi-1];rt=[maxi-1:-1:1];erg(1,maxi)=0;if (size(e,1) == 1) erg(1,[maxi-1:-1:1])=-N*t.*(cumsum(log(e([maxi:-1:2])))./t-log(cumsum(e([maxi:-1:2]))./t));else x=triu(e,1); y=x; x(x==0)=1; erg(1,t)=-N*rt.*(log(prod(x(1:maxi-1,:),2)')./rt-log(sum(y(1:maxi-1,:),2)'./rt));enderg(2,:)=M*log(N)/2;switch mdltype case 'p'  if (size(e,1) == 1)   erg(3,:)=M.*(log(gamma)-(cumsum(log(e))./i));  else   x=tril(e);   x(x==0)=1;   erg(3,:)=M.*(log(gamma)-(log(prod(x,2)')./i));  end case 'r'  erg(3,:)=M.*log(gamma); otherwise  erg(3,:)=0;enderg=sum(erg(:,i),1);% find minimum[vK,K]=min(erg);if (verbose>1) warning(['MDL: ' sprintf('%d ',N,maxi) '-> ' sprintf('%f ',erg)]); if (verbose>7) cplot('split',erg1,erg2,erg3); title('MDL'); pause; end if (verbose>5) cplot('combined',e/max(e),erg/max(erg)); title('MDL'); pause; endendreturn