Daniel@0: function varargout = mirstat(f,varargin)
Daniel@0: % stat = mirstat(f) returns basic statistics of the feature f as a
Daniel@0: %   structure array stat such that:
Daniel@0: %       stat.Mean is the mean of the feature;
Daniel@0: %       stat.Std is the standard deviation;
Daniel@0: %       stat.Slope is the linear slope of the curve;
Daniel@0: %       stat.PeriodFreq is the frequency (in Hz) of the main periodicity;
Daniel@0: %       stat.PeriodAmp is the normalized amplitude of the main periodicity;
Daniel@0: %       stat.PeriodEntropy is the entropy of the periodicity curve.
Daniel@0: %   The main periodicity and periodicity curve are evaluated through the
Daniel@0: %       computation of the autocorrelation sequence related to f.
Daniel@0: %
Daniel@0: %   f can be itself a structure array composed of features. In this case,
Daniel@0: %       stat will be structured the same way.
Daniel@0: %
Daniel@0: % mirstat does not work for multi-channels objects.
Daniel@0: 
Daniel@0: if isa(f,'mirstruct')
Daniel@0:     varargout = {set(f,'Stat',1)};
Daniel@0: elseif isstruct(f)
Daniel@0:     if isdesign(f)
Daniel@0:         f.Stat = 1;
Daniel@0:         varargout = {f};
Daniel@0:     else
Daniel@0:         fields = fieldnames(f);
Daniel@0:         for i = 1:length(fields)
Daniel@0:             field = fields{i};
Daniel@0:             ff = f.(field);
Daniel@0:             if iscell(ff) && isa(ff{1},'mirdata')
Daniel@0:                 ff = ff{1};
Daniel@0:             end
Daniel@0:             if isa(ff,'mirdata')
Daniel@0:                 if i == 1 
Daniel@0:                     la = get(ff,'Label');
Daniel@0:                     if not(isempty(la) || isempty(la{1}))
Daniel@0:                         stat.Class = la;
Daniel@0:                     end
Daniel@0:                 end
Daniel@0:                 ff = set(ff,'Label','');
Daniel@0:             end
Daniel@0:             stat.(field) = mirstat(ff,'FileNames',0);
Daniel@0:         end
Daniel@0:         if isempty(varargin)
Daniel@0:             f0 = f;
Daniel@0:             while isstruct(f0)
Daniel@0:                 fields = fieldnames(f0);
Daniel@0:                 f0 = f0.(fields{1});
Daniel@0:             end
Daniel@0:             if iscell(f0)
Daniel@0:                 f0 = f0{1};
Daniel@0:             end
Daniel@0:             stat.FileNames = get(f0,'Name');
Daniel@0:         end
Daniel@0:         varargout = {stat};
Daniel@0:     end
Daniel@0: elseif iscell(f)
Daniel@0:     if 0 %ischar(f{1})
Daniel@0:         varargout = {f};
Daniel@0:     else
Daniel@0:         res = zeros(length(f),1);
Daniel@0:         for i = 1:length(f)
Daniel@0:             res(i) = mirstat(f{i});
Daniel@0:         end
Daniel@0:         varargout = {res};
Daniel@0:     end
Daniel@0: elseif isnumeric(f)
Daniel@0:     f(isnan(f)) = [];
Daniel@0:     varargout = {mean(f)};
Daniel@0: else
Daniel@0:             alongfiles.key = 'AlongFiles';
Daniel@0:             alongfiles.type = 'Boolean';
Daniel@0:             alongfiles.default = 0;
Daniel@0:         option.alongfiles = alongfiles;
Daniel@0: 
Daniel@0:             filenames.key = 'FileNames';
Daniel@0:             filenames.type = 'Boolean';
Daniel@0:             filenames.default = 1;
Daniel@0:         option.filenames = filenames;
Daniel@0:         
Daniel@0:     specif.option = option;
Daniel@0: 
Daniel@0:     specif.nochunk = 1;
Daniel@0:     varargout = mirfunction(@mirstat,f,varargin,nargout,specif,@init,@main);
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: function [x type] = init(x,option)
Daniel@0: type = '';
Daniel@0: 
Daniel@0: 
Daniel@0: function stat = main(f,option,postoption)
Daniel@0: if iscell(f)
Daniel@0:     f = f{1};
Daniel@0: end
Daniel@0: if isa(f,'mirhisto')
Daniel@0:     warning('WARNING IN MIRSTAT: histograms are not taken into consideration yet.')
Daniel@0:     stat = struct;
Daniel@0:     return
Daniel@0: end
Daniel@0: fp = get(f,'FramePos');
Daniel@0: if haspeaks(f)
Daniel@0:     ppp = get(f,'PeakPrecisePos');
Daniel@0:     if not(isempty(ppp)) && not(isempty(ppp{1}))
Daniel@0:         stat = addstat(struct,ppp,fp,'PeakPos');
Daniel@0:         stat = addstat(stat,get(f,'PeakPreciseVal'),fp,'PeakMag');
Daniel@0:     else
Daniel@0:         if isa(f,'mirkeystrength') || (isa(f,'mirchromagram') && get(f,'Wrap'))
Daniel@0:             stat = struct; % This needs to be defined using some kind of circular statistics..
Daniel@0:         else
Daniel@0:             stat = addstat(struct,get(f,'PeakPosUnit'),fp,'PeakPos');
Daniel@0:         end
Daniel@0:         stat = addstat(stat,get(f,'PeakVal'),fp,'PeakMag');
Daniel@0:     end
Daniel@0: else
Daniel@0:     stat = addstat(struct,get(f,'Data'),fp,'',option.alongfiles,...
Daniel@0:                           get(f,'Name'),get(f,'Label'));
Daniel@0: end
Daniel@0: if option.filenames
Daniel@0:     stat.FileNames = get(f,'Name');
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: function stat = addstat(stat,d,fp,field,alongfiles,filename,labels)
Daniel@0: l = length(d);
Daniel@0: if nargin<5
Daniel@0:     alongfiles = 0;
Daniel@0: end
Daniel@0: if nargin<7
Daniel@0:     labels = {};
Daniel@0: end
Daniel@0: if not(alongfiles) || l == 1
Daniel@0:     anyframe = 0;
Daniel@0:     for i = 1:l
Daniel@0:         if not(iscell(d{i}))
Daniel@0:             d{i} = {d{i}};
Daniel@0:         end
Daniel@0:         for k = 1:length(d{i})
Daniel@0:             dd = d{i}{k};
Daniel@0:             if iscell(dd)
Daniel@0:                 if length(dd)>1
Daniel@0:                     anyframe = 1;
Daniel@0:                 end
Daniel@0:                 dn = cell2array(dd);
Daniel@0:                 [m{i}{k},s{i}{k},sl{i}{k},pa{i}{k},pf{i}{k},pe{i}{k}] ...
Daniel@0:                     = computestat(dn,fp{i}{k});
Daniel@0:             elseif size(dd,2) < 2
Daniel@0:                 nonan = find(not(isnan(dd)));
Daniel@0:                 dn = dd(nonan);
Daniel@0:                 if isempty(dn)
Daniel@0:                     m{i}{k} = NaN;
Daniel@0:                 else
Daniel@0:                     m{i}{k} = mean(dn,2);
Daniel@0:                 end
Daniel@0:                 s{i}{k} = NaN;
Daniel@0:                 sl{i}{k} = NaN;
Daniel@0:                 pa{i}{k} = NaN;
Daniel@0:                 pf{i}{k} = NaN;
Daniel@0:                 pe{i}{k} = NaN;
Daniel@0:             else
Daniel@0:                 anyframe = 1;
Daniel@0:                 s1 = size(dd,1);
Daniel@0:                 s3 = size(dd,3);
Daniel@0:                 dd = mean(dd,4); %mean(dd,3),4);
Daniel@0:                 m{i}{k} = NaN(s1,1,s3);
Daniel@0:                 s{i}{k} = NaN(s1,1,s3);
Daniel@0:                 sl{i}{k} = NaN(s1,1,s3);
Daniel@0:                 pa{i}{k} = NaN(s1,1,s3);
Daniel@0:                 pf{i}{k} = NaN(s1,1,s3);
Daniel@0:                 pe{i}{k} = NaN(s1,1,s3);
Daniel@0:                 for h = 1:s3
Daniel@0:                     for j = 1:s1
Daniel@0:                         [m{i}{k}(j,1,h),s{i}{k}(j,1,h),sl{i}{k}(j,1,h),...
Daniel@0:                          pa{i}{k}(j,1,h),pf{i}{k}(j,1,h),pe{i}{k}(j,1,h)] = ...
Daniel@0:                             computestat(dd(j,:,h),fp{i}{k});
Daniel@0:                     end
Daniel@0:                 end
Daniel@0:             end
Daniel@0:         end
Daniel@0:     end
Daniel@0:     if anyframe
Daniel@0:         fields = {'Mean','Std','Slope','PeriodFreq','PeriodAmp','PeriodEntropy'};
Daniel@0:         stats = {m,s,sl,pf,pa,pe};   
Daniel@0:     else
Daniel@0:         fields = {'Mean'};
Daniel@0:         stats = {m};   
Daniel@0:     end
Daniel@0:     for i = 1:length(stats)
Daniel@0:         data = stats{i};
Daniel@0:         data = uncell(data,NaN);
Daniel@0:         stat.(strcat(field,fields{i})) = data;
Daniel@0:     end
Daniel@0: else
Daniel@0:     if iscell(d{1}{1})
Daniel@0:         slash = strfind(filename{1},'/');
Daniel@0:         nbfolders = 1;
Daniel@0:         infolder = 0;
Daniel@0:         foldername{1} = filename{1}(1:slash(end)-1);
Daniel@0:         for i = 1:length(d)
Daniel@0:             slash = strfind(filename{i},'/');
Daniel@0:             if not(strcmpi(filename{i}(1:slash(end)-1),foldername))
Daniel@0:                 nbfolders = nbfolders + 1;
Daniel@0:                 infolder = 0;
Daniel@0:                 foldername{nbfolders} = filename{i}(1:slash(end)-1);
Daniel@0:             end
Daniel@0:             infolder = infolder+1;
Daniel@0:             dd{nbfolders}(infolder,:) = cell2array(d{i}{1});
Daniel@0:         end
Daniel@0:         for i = 1:length(dd)
Daniel@0:             figure
Daniel@0:             plot(dd{i}')
Daniel@0:             stat.Mean{i} = mean(dd{i});
Daniel@0:             stat.Std{i} = std(dd{i});
Daniel@0:         end
Daniel@0:     end
Daniel@0: end
Daniel@0: if not(isempty(labels) || isempty(labels{1}))
Daniel@0:     stat.Class = labels;
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: function dn = cell2array(dd)
Daniel@0: dn = zeros(1,length(dd));
Daniel@0: for j = 1:length(dd)
Daniel@0:     if isempty(dd{j})
Daniel@0:         dn(j) = NaN;
Daniel@0:     else
Daniel@0:         dn(j) = dd{j}(1);
Daniel@0:     end
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: function [m,s,sl,pa,pf,pe] = computestat(d,fp)
Daniel@0: m = NaN;
Daniel@0: s = NaN;
Daniel@0: sl = NaN;
Daniel@0: pa = NaN;
Daniel@0: pf = NaN;
Daniel@0: pe = NaN;
Daniel@0: diffp = fp(1,2:end) - fp(1,1:end-1);
Daniel@0: if isempty(diffp) || sum(round((diffp(2:end)-diffp(1:end-1))*1000))
Daniel@0:     % Not regular sampling (in mirattacktime for instance)
Daniel@0:     framesampling = NaN;
Daniel@0: else
Daniel@0:     framesampling = fp(1,2)-fp(1,1);
Daniel@0: end
Daniel@0: nonan = find(not(isnan(d)));
Daniel@0: if not(isempty(nonan))
Daniel@0:     dn = d(nonan);
Daniel@0:     m = mean(dn,2);
Daniel@0:     s = std(dn,0,2);
Daniel@0:     if not(isnan(s))
Daniel@0:         if s
Daniel@0:             dk = (dn-m)/s;
Daniel@0:             tk = linspace(0,1,size(d,2));
Daniel@0:             sl = dk(:)'/tk(nonan);
Daniel@0:         elseif size(s,2) == 1
Daniel@0:             s = NaN;
Daniel@0:         end
Daniel@0:     end
Daniel@0:     if length(dn)>1
Daniel@0:         cor = xcorr(dn',dn','coeff');
Daniel@0:         cor = cor(ceil(length(cor)/2):end);
Daniel@0:         % let's zero the first descending slope of the
Daniel@0:         % autocorrelation function
Daniel@0:         firstmin = find(cor(2:end)>cor(1:end-1));
Daniel@0:         if not(isempty(firstmin) || isnan(framesampling))
Daniel@0:             cor2 = cor;
Daniel@0:             cor2(1:firstmin(1)) = 0;
Daniel@0:             [pa,pfk] = max(cor2);
Daniel@0:             if pfk > 1
Daniel@0:                 pf = 1/(pfk-1)/framesampling;
Daniel@0:             end
Daniel@0:         end
Daniel@0:         cor = abs(cor);
Daniel@0:         cor = cor/sum(cor);
Daniel@0:         pe = -sum(cor.*log(cor+1e-12))./log(length(cor));
Daniel@0:     end
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: function b = isdesign(f)
Daniel@0: if iscell(f)
Daniel@0:     f = f{1};
Daniel@0: end
Daniel@0: if isa(f,'mirdesign') || isa(f,'mirstruct')
Daniel@0:     b = 1;
Daniel@0: elseif isa(f,'mirdata') || not(isstruct(f))
Daniel@0:     b = 0;
Daniel@0: else
Daniel@0:     fields = fieldnames(f);
Daniel@0:     b = isdesign(f.(fields{1}));
Daniel@0: end