function MMA(recordName)
% % Multiscale Multifractal Analysis
% % Copyright (C) 2014 Jan Gieraltowski
% %
% % This program is free software; you can redistribute it and/or modify it under
% % the terms of the GNU General Public License as published by the Free Software
% % Foundation; either version 2 of the License, or (at your option) any later
% % version.
% %
% % This program is distributed in the hope that it will be useful, but WITHOUT ANY
% % WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
% % PARTICULAR PURPOSE. See the GNU General Public License for more details.
% %
% % You should have received a copy of the GNU General Public License along with
% % this program; if not, write to the Free Software Foundation, Inc., 59 Temple
% % Place - Suite 330, Boston, MA 02111-1307, USA.
% %
% % Author: Jan Gieraltowski
% % Warsaw University of Technology
% % Faculty of Physics
% % gieraltowski@if.pw.edu.pl
% % http://gieraltowski.fizyka.pw.edu.pl/
% %
% % Method was first proposed in:
% % J. Gieraltowski, J. J. Zebrowski, and R. Baranowski,
% % Multiscale multifractal analysis of heart rate variability recordings with a large number of occurrences of arrhythmia,
% % Phys. Rev. E 85, 021915 (2012).
% % http://dx.doi.org/10.1103/PhysRevE.85.021915
% %
% % Please cite the above publication when referencing this material,
% % and also include the standard citation for PhysioNet:
% % Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG, Mietus JE, Moody GB, Peng C-K, Stanley HE.
% % PhysioBank, PhysioToolkit, and PhysioNet: Components of a New Research Resource for Complex Physiologic Signals.
% % Circulation 101(23):e215-e220
% % [Circulation Electronic Pages; http://circ.ahajournals.org/cgi/content/full/101/23/e215]; 2000 (June 13).
overlapping = 0; %0 - time series is partitioned into no overlapping windows of analysis, 1- time series is partitioned into overlapping windows of analysis, step between consecutive windows is = 1 (much longer calculations)
smin = 10; %minimal s scale used, when calculating Fq(s) functions family (default 10)
smax = 600; %maximal s scale used, when calculating Fq(s) functions family, has to be multiple of 5 (default 600; in general should be near to N/50, where N is a time series length)
qmin = -5; %minimal multifractal parameter q used (default -5)
qmax = 5; %maximal multifractal parameter q used (deafault 5)
qlist = qmin:0.1:qmax; qlist(qlist == 0) = 0.0001;
precisionMode = 0; % 0 (default) better looking plot, smaller files, faster calculations; set 1 for enhanced precision (smaller q and s steps)
[pathstr,name] = fileparts(recordName);
warning('off','MATLAB:MKDIR:DirectoryExists');
if isdir(recordName) == 1
dirtemp = dir([recordName filesep '*.txt']);
filenames = {dirtemp.name}';
filenum = size(filenames,1);
dirname = [pathstr filesep name];
else
filenames = [name '.txt'];
filenum = 1;
dirname = pathstr;
end
for fileit = 1:filenum
signal = load([dirname filesep char(filenames(fileit,:))]);
prof = cumsum(signal);
slength = size(prof);
fqs = [];
for s = smin:1:smax
if overlapping == 1
vec = [0:s-1];
ind = [1:slength-s+1]';
coordinates = bsxfun(@plus, vec, ind);
else
ind2 = [1:size(prof,1)];
coordinates = reshape(ind2(1:(size(prof,1)-mod(size(prof,1),s))),s,(size(prof,1)-mod(size(prof,1),s))/s)';
end
segments = prof(coordinates);
xbase = [1:1:s];
f2nis = [];
for ni = 1:size(segments,1)
seg = segments(ni,:);
fit = polyfit(xbase,seg,2);
variance = mean((seg - polyval(fit,xbase)).^2);
f2nis(end+1) = variance;
end
for q = qlist
fqs = [fqs; q s (mean(f2nis.^(q/2)))^(1/q)];
end
end
fqsll = [fqs(:,1) fqs(:,2) log(fqs(:,2)) log(fqs(:,3))];
hqs = [];
if precisionMode == 0
sspacing = ((smax/5)-smin)/11;
qlist = qmin:1:qmax; qlist(qlist == 0) = 0.0001;
else
sspacing = 1;
end
for sit = smin:sspacing:(smax/5)
for qit = qlist
fittemp = fqsll(fqsll(:,1) == qit & fqsll(:,2) >= sit & fqsll(:,2) <= 5*sit,:);
htemp = polyfit(fittemp(:,3),fittemp(:,4),1);
hqs = [hqs; qit 3*sit htemp(1)];
end
end
mkdir([dirname filesep 'MMA_results']);
[p1,n1,e1] = fileparts(char(filenames(fileit,:)));
csvwrite([dirname filesep 'MMA_results' filesep 'MMA_' n1 '.txt'],hqs);
hqsplotdata = reshape(hqs(:,3),size(qlist,2),size(smin:sspacing:(smax/5),2));
if max(max(hqsplotdata)) < 1.5
hlim = 1.5;
elseif max(max(hqsplotdata)) < 2.5
hlim = 2.5;
else
hlim = ceil((max(max(hqsplotdata))*10))/10;
end
hqsplot = surf(unique(hqs(:,2)),unique(hqs(:,1)),hqsplotdata);
colormap(jet);
colorbar;
caxis([0,hlim]);
set(gca,'YDir','reverse');
view(-62,50);
axis([3*smin,0.6*smax,qmin,qmax,0,hlim]);
saveas(hqsplot, [dirname filesep 'MMA_results' filesep 'MMA_' n1 '.jpg']);
end
end