function H = hurst_psd(x);
%
% Compute the Hurst parameter for a sequence using the 
% power spectral density method
%
%


%
% Compute log(|F(omega)|) and log(omega) sequences
%
%
fx = abs(fft(x));
y=log(fx(1:floor(length(fx)/2)));
y=y';                            % Shoot me now
x=log(1:length(y));
clear fx;

%
% Linear fit to those sequences
%
%

mx=mean(x);
my=mean(y);
n=length(y);

b1 = (sum(x.*y) - n*mx*my) / (sum(x.^2) - n*mx^2);
b0 = my - b1*mx;

%plot(x,y);
%line([x(1),x(n)],[b0,b0+b1*((x(n)-x(1)))]);

beta=b1;

H=(1-beta)/2;