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