Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Higuchi's fractal dimension code

Subject: Higuchi's fractal dimension code

From: Tikkuhirvi Tietavainen

Date: 22 Apr, 2008 10:46:02

Message: 1 of 5

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

Subject: Higuchi's fractal dimension code

From: Steen Johansen

Date: 27 May, 2009 08:54:01

Message: 2 of 5

> 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.

Subject: Higuchi's fractal dimension code

From: Sarah

Date: 26 Feb, 2010 10:30:22

Message: 3 of 5

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

Subject: Higuchi's fractal dimension code

From: Aino

Date: 4 Jun, 2010 01:25:06

Message: 4 of 5

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.menetre@schiller.fr> wrote in message <hm87ru$g4g$1@fred.mathworks.com>...
> 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

Subject: Higuchi's fractal dimension code

From: Salai Selvam V

Date: 19 Jan, 2011 16:57:05

Message: 5 of 5

"Aino" <aino.tietavainen@removeThis.helsinki.fi> wrote in message <hu9kli$qi2$1@fred.mathworks.com>...
> 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.menetre@schiller.fr> wrote in message <hm87ru$g4g$1@fred.mathworks.com>...
> > 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
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: vsalaiselvam@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!!!****************

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us