function [H_fft,sdh1p,sdh1mm,D_fft] = frac_fft(p)

%% License:
% "This work is licensed under the Creative Commons Attribution 4.0
% International (CC BY 4.0) License. To view a copy of this license, visit:
% http://creativecommons.org/licenses/by/4.0/"

%% please cite 
% --> Kristofer Marsch (2020)
% --> DOI for this data available at https://orcid.org/0000-0002-1857-506X


%% calculating Hurst Exponent, sdh1p, sdh1mm and Fractal Dimension using FFT
% based on Stigsson & Mas Ivars (2019)

%% --parameters-- %%
% p: array of profile, first column horizontal values, second column heights
% profile must be de-trended
% profile must have equal point distances, e.g. constant sampling interval

%% calculate H and D
yP = p(:,2);
transyP = fft(yP);

freq = linspace(0,length(transyP)/2-1,length(transyP)/2); %half the population
powersub = abs(transyP).^2;
power = powersub(1:length(freq),:); %powers according to frequency vector

X = log10((freq(:,2:end-1))'); %don't use the first and Nyquist frequency
Y = log10(power(2:end-1,:));
X2 = [ones(length(X),1) X];
c = X2\Y;

H_fft = ((-c(2))-1)/2;
D_fft = 2-H_fft;

%% calculate sdh1p and sdh1mm
intersec = 10^(c(1));
sub1 = (freq.^(-(H_fft+0.5)).*sin(pi.*freq/length(yP))).^2;
sub2 = sum(sub1(:,2:end-1));
sdh1p = 2*(2^0.5)/length(yP)*(intersec^0.5)*(sub2^0.5);

sampstep = p(2,1)-p(1,1);
sdh1mm = sdh1p/(sampstep/1)^H_fft;

