Higuchi's fractal dimension code

[Tikkuhirvi Tietavainen]

Hello all,

I have tried to calculate a fractal dimension (FD) trace
using Higuchi's algorithm. My problem is very simple: the
code doesn't work, since it doesn't provide a FD between 1-2
(stays under 1). There are two M-files underneath, first is
the Higushi code and the second is a fractal trace
(Weierstrass, FD=1.5), which can be used to test the Higuchi
code.

Could someone please tell me what seems to be the problem
with the code? Thank you!

function Higuchi_v4(trace)
%This function calculates a fractal dimension trace using
Higuchi's
%[Higuchi 1988] method.
%First windowsize/2 and last windowsize/2 datapoinst are zeros

windowsize=200;%From how many datapoints the FD is calculated
format long
for v=1:length(trace)-windowsize;

w=v+windowsize;
N=windowsize;
COP=trace(v:w);%Takes a piece of the trace, w-v datapoints.
kmax=50;
L_m=zeros(1,kmax);
L_m_length=zeros(kmax,kmax);

for k=1:kmax
for m=1:k;
for i=1:fix((N-m)/k)
L_m(i)=abs(COP(m+i*k)-COP(m+(i-1)*k));
end
a=(N-1)/(fix((N-m)/k)*k);
L_m_length(m,k)=(sum(L_m)*a)/k;
end
end

for j=1:kmax

L_m_length_mean(j)=mean(nonzeros(L_m_length(:,j)));%Extra
zeros from the matrix
L_m_length_std(j)=std(nonzeros(L_m_length(:,j)));
end

if v==length(trace)-windowsize;
figure(1);plot(1:kmax,L_m_length_mean,'*')

figure(2);errorbar(1:kmax,L_m_length_mean,L_m_length_std)
end

k=1:kmax;
p=polyfit(log(1./k), log(L_m_length_mean), 1);
%FD(N)=p(1);%Fractal dimension is the fitted slope

L_m_windowed(v+(w-v)/2)=p(1);%When FD is calculated
between w and v, the FD value will be between them.
end

figure;plot(L_m_windowed)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [W] = Wei()
%Weierstrass fractal trace

H=0.5; %FD should be 2-H.
y=5;

for t=1:1000;
for i=1:100
W_apu(i,t)=y^(-i*H)*cos(2*3.1416*y^i*t);
end
W(t)=sum(W_apu(:,t));
end

[Steen Johansen]

> Could someone please tell me what seems to be the problem
> with the code? Thank you!

Hello Tikkuhirvi

You may want to review your formula for calculation of the Higuichi dimension.
I think the dimension is the slope of the inverse function you are using.
(Biomedical signal and image processing; Kayvan Najarian, Robert Splinter).

Try using:
p=polyfit(log(L_m_length_mean),log(1./k), 1);

This will give you values between one and two.

[Sarah]

Hello Tikkuhirvi,

I would be interested to know if you have solved your problem exposed here about the Higuchi's fractal dimension code. Was the proposition from Steen Johansen the solution you adopted ?
Because I found nowhere (neither in the Biomedical signal and image processing book from Kayvan Najarian, Robert Splinter nor on the websites), something like "dimension is the slope of the inverse function you are using" !

Thank you in advance,

Sarah

[Aino]

function [FD]=Higuchi(COP)
%This function calculates a fractal dimension of a COP trace using
%Higuchi's [Higuchi 1988] method.

N=length(COP);
kmax=100;
L_m=zeros(1,kmax);
L_m_length=zeros(kmax,kmax);

for k=1:kmax
for m=1:k;

L_m=zeros(1,kmax);

for i=1:fix((N-m)/k)
L_m(i)=abs(COP(m+i*k)-COP(m+(i-1)*k));
end

a=((N-1))/(fix((N-m)/k)*k);
L_m_length(m,k)=(sum(L_m)*a)/k;
end
L_m_length_mean(k)=mean(nonzeros(L_m_length(:,k))); %Extra zeros from the matrix
L_m_length_std(k)=std(nonzeros(L_m_length(:,k)));
end

k=1:kmax;
p=polyfit(log(1./k), log(L_m_length_mean), 1);

FD=p(1);%FD is the slope

________________________________________

I sure hope it does what it is supposed to do. But I now have a more resent problem concerning it. The kmax is supposed to be found by plotting FDs that are acquired with different kmaxs against kmaxs. The suitable kmax should be the one at what point the plot plateaus. First of all, does that mean that the FD increases as the kmax increases? And second of all, what if (as in my case) the FD get all the way to 2 before it plateaus? The longer the trace is, the longer the kmax needed, it would seem..

Thank you for your answers.

-Aino

"Sarah " <sarah....@schiller.fr> wrote in message <hm87ru$g4g$1...@fred.mathworks.com>...

[Salai Selvam V]

"Aino" <aino.tie...@removeThis.helsinki.fi> wrote in message <hu9kli$qi2$1...@fred.mathworks.com>...

Dear Aino, Sarah and all others, The mistake you all made was the wrong algorithm. For exact algorithm please refer to Higuchi's original work (even the mistake is there even in IEEE paper by Rosana Esteller et al but the paper by Polychronaki et al explains the Higuchi FD algorithm correctly). The mistake is the curve length formula, which contains the term, k dividing the entire value. Here is the correct MATLAB Code: Please download it from MATLAB Central. If you think that I have solved your problem, please mail me: vsalai...@yahoo.com. Thanks to all.
***********COMPLETE MATLAB CODE FOR HIGUCHI FD!!!****************
function [varargout]=hfd(x,kmax)

if ~exist('kmax','var')||isempty(kmax),
kmax=5;
end;

N=length(x);

Lmk=zeros(kmax,kmax);
for k=1:kmax,
for m=1:k,
Lmki=0;
for i=1:fix((N-m)/k),
Lmki=Lmki+abs(x(m+i*k)-x(m+(i-1)*k));
end;
Ng=(N-1)/(fix((N-m)/k)*k);
Lmk(m,k)=(Lmki*Ng)/k;%Here is the problem in your code, Mr. Tikkuhirvi & Mr. Aino
end;
end;

Lk=zeros(1,kmax);
for k=1:kmax,
Lk(1,k)=sum(Lmk(1:k,k))/k;
end;

lnLk=log(Lk);
lnk=log(1./[1:kmax]);

b=polyfit(lnk,lnLk,1);
xhfd=b(1);

varargout={xhfd,lnk,lnLk,Lk,Lmk};
***********COMPLETE MATLAB CODE FOR HIGUCHI FD!!!****************

[anadi.j...@gmail.com]

please sir, give me the codes to find kmax