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