dcorr.m

biomedical stuffs...

%	dcorr (a modification of corrdim) uses:
%	x = data "vector" of size p x n
%	n = dimension of each element (i.e., row) of x
%	n is the embedding dimension
%	p = length of x (number of data points)
%	calculates regression for several possible ranges
%	and prints results on screen (capture with diary)
%	saving of metafile has been commented out
%
%	This program is distributed as a supplement to the book
%	"Biomedical Signal Processing and Signal Modeling" by E. N. Bruce,
%	published by Wiley, 2000.  It is provided for educational use only.
%  While every effort has been made to insure its suitability to illustrate
%  principles described in the above book, no specific feature or capability 
%  is implied or guaranteed.
	for j=1:p
a(j)=norm(x(j,:));
	end
rng=max(a)
clear a
levels=zeros(1,50);
%eps=2^k*rng
%kmax=1 -> counts(50)
%kmin=-48 -> counts(1)
levels(1)=2^(-48)*rng;
	for j=2:50
levels(j)=2*levels(j-1);
	end
ceps=zeros(1,50);
	for j=1:p-1
nj=p-j;
ntot=0;
w1=x(j+1:p,:);
w=(w1-ones(nj,1)*x(j,:)).^2;
w2=[w zeros(nj,1)];
d=sum(w2');
		for k=1:50
over = d <= (levels(k))^2*ones(1,nj);
ntot=sum(over);
ceps(k)=ceps(k)+ntot;
		end
	end
psq=p*(p-1)/2
ceps=ceps/psq;
loglog(levels,ceps,'o')
grid
xlabel('levels')
ylabel('normalized counts')
%meta dcorr
%	find range where C(eps) is >0 and <1
levs=1:1:50;
g=levs(ceps>0&ceps<1);
%	find where C(eps) starts to increase
lg=length(g);
ming=1; i=1;
	while i<lg
if ceps(g(i)+1)>ceps(g(i)) ming=i; i=lg; end
i=i+1;
	end % of i loop
%  calculate regressions for various combinations of points
disp('   minlev   maxlev   slope   intercept')
minlev=min(g(ming)+2,g(lg)-2);
maxlev=max(g(lg)-2,minlev+1);
v1=[(log(levels(minlev:maxlev)))',(ones(size(levels(minlev:maxlev))))'];
v2=(log(ceps(minlev:maxlev)))';
dc=v1\v2; out=[minlev maxlev dc']; disp(out)