Mercurial > hg > camir-aes2014

annotate toolboxes/MIRtoolbox1.3.2/MIRToolbox/mirstat.m @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 function varargout = mirstat(f,varargin)
wolffd@0 2 % stat = mirstat(f) returns basic statistics of the feature f as a
wolffd@0 3 % structure array stat such that:
wolffd@0 4 % stat.Mean is the mean of the feature;
wolffd@0 5 % stat.Std is the standard deviation;
wolffd@0 6 % stat.Slope is the linear slope of the curve;
wolffd@0 7 % stat.PeriodFreq is the frequency (in Hz) of the main periodicity;
wolffd@0 8 % stat.PeriodAmp is the normalized amplitude of the main periodicity;
wolffd@0 9 % stat.PeriodEntropy is the entropy of the periodicity curve.
wolffd@0 10 % The main periodicity and periodicity curve are evaluated through the
wolffd@0 11 % computation of the autocorrelation sequence related to f.
wolffd@0 12 %
wolffd@0 13 % f can be itself a structure array composed of features. In this case,
wolffd@0 14 % stat will be structured the same way.
wolffd@0 15 %
wolffd@0 16 % mirstat does not work for multi-channels objects.
wolffd@0 17
wolffd@0 18 if isa(f,'mirstruct')
wolffd@0 19 varargout = {set(f,'Stat',1)};
wolffd@0 20 elseif isstruct(f)
wolffd@0 21 if isdesign(f)
wolffd@0 22 f.Stat = 1;
wolffd@0 23 varargout = {f};
wolffd@0 24 else
wolffd@0 25 fields = fieldnames(f);
wolffd@0 26 for i = 1:length(fields)
wolffd@0 27 field = fields{i};
wolffd@0 28 ff = f.(field);
wolffd@0 29 if iscell(ff) && isa(ff{1},'mirdata')
wolffd@0 30 ff = ff{1};
wolffd@0 31 end
wolffd@0 32 if isa(ff,'mirdata')
wolffd@0 33 if i == 1
wolffd@0 34 la = get(ff,'Label');
wolffd@0 35 if not(isempty(la) || isempty(la{1}))
wolffd@0 36 stat.Class = la;
wolffd@0 37 end
wolffd@0 38 end
wolffd@0 39 ff = set(ff,'Label','');
wolffd@0 40 end
wolffd@0 41 stat.(field) = mirstat(ff,'FileNames',0);
wolffd@0 42 end
wolffd@0 43 if isempty(varargin)
wolffd@0 44 f0 = f;
wolffd@0 45 while isstruct(f0)
wolffd@0 46 fields = fieldnames(f0);
wolffd@0 47 f0 = f0.(fields{1});
wolffd@0 48 end
wolffd@0 49 if iscell(f0)
wolffd@0 50 f0 = f0{1};
wolffd@0 51 end
wolffd@0 52 stat.FileNames = get(f0,'Name');
wolffd@0 53 end
wolffd@0 54 varargout = {stat};
wolffd@0 55 end
wolffd@0 56 elseif iscell(f)
wolffd@0 57 if 0 %ischar(f{1})
wolffd@0 58 varargout = {f};
wolffd@0 59 else
wolffd@0 60 res = zeros(length(f),1);
wolffd@0 61 for i = 1:length(f)
wolffd@0 62 res(i) = mirstat(f{i});
wolffd@0 63 end
wolffd@0 64 varargout = {res};
wolffd@0 65 end
wolffd@0 66 elseif isnumeric(f)
wolffd@0 67 f(isnan(f)) = [];
wolffd@0 68 varargout = {mean(f)};
wolffd@0 69 else
wolffd@0 70 alongfiles.key = 'AlongFiles';
wolffd@0 71 alongfiles.type = 'Boolean';
wolffd@0 72 alongfiles.default = 0;
wolffd@0 73 option.alongfiles = alongfiles;
wolffd@0 74
wolffd@0 75 filenames.key = 'FileNames';
wolffd@0 76 filenames.type = 'Boolean';
wolffd@0 77 filenames.default = 1;
wolffd@0 78 option.filenames = filenames;
wolffd@0 79
wolffd@0 80 specif.option = option;
wolffd@0 81
wolffd@0 82 specif.nochunk = 1;
wolffd@0 83 varargout = mirfunction(@mirstat,f,varargin,nargout,specif,@init,@main);
wolffd@0 84 end
wolffd@0 85
wolffd@0 86
wolffd@0 87 function [x type] = init(x,option)
wolffd@0 88 type = '';
wolffd@0 89
wolffd@0 90
wolffd@0 91 function stat = main(f,option,postoption)
wolffd@0 92 if iscell(f)
wolffd@0 93 f = f{1};
wolffd@0 94 end
wolffd@0 95 if isa(f,'mirhisto')
wolffd@0 96 warning('WARNING IN MIRSTAT: histograms are not taken into consideration yet.')
wolffd@0 97 stat = struct;
wolffd@0 98 return
wolffd@0 99 end
wolffd@0 100 fp = get(f,'FramePos');
wolffd@0 101 if haspeaks(f)
wolffd@0 102 ppp = get(f,'PeakPrecisePos');
wolffd@0 103 if not(isempty(ppp)) && not(isempty(ppp{1}))
wolffd@0 104 stat = addstat(struct,ppp,fp,'PeakPos');
wolffd@0 105 stat = addstat(stat,get(f,'PeakPreciseVal'),fp,'PeakMag');
wolffd@0 106 else
wolffd@0 107 if isa(f,'mirkeystrength') || (isa(f,'mirchromagram') && get(f,'Wrap'))
wolffd@0 108 stat = struct; % This needs to be defined using some kind of circular statistics..
wolffd@0 109 else
wolffd@0 110 stat = addstat(struct,get(f,'PeakPosUnit'),fp,'PeakPos');
wolffd@0 111 end
wolffd@0 112 stat = addstat(stat,get(f,'PeakVal'),fp,'PeakMag');
wolffd@0 113 end
wolffd@0 114 else
wolffd@0 115 stat = addstat(struct,get(f,'Data'),fp,'',option.alongfiles,...
wolffd@0 116 get(f,'Name'),get(f,'Label'));
wolffd@0 117 end
wolffd@0 118 if option.filenames
wolffd@0 119 stat.FileNames = get(f,'Name');
wolffd@0 120 end
wolffd@0 121
wolffd@0 122
wolffd@0 123 function stat = addstat(stat,d,fp,field,alongfiles,filename,labels)
wolffd@0 124 l = length(d);
wolffd@0 125 if nargin<5
wolffd@0 126 alongfiles = 0;
wolffd@0 127 end
wolffd@0 128 if nargin<7
wolffd@0 129 labels = {};
wolffd@0 130 end
wolffd@0 131 if not(alongfiles) || l == 1
wolffd@0 132 anyframe = 0;
wolffd@0 133 for i = 1:l
wolffd@0 134 if not(iscell(d{i}))
wolffd@0 135 d{i} = {d{i}};
wolffd@0 136 end
wolffd@0 137 for k = 1:length(d{i})
wolffd@0 138 dd = d{i}{k};
wolffd@0 139 if iscell(dd)
wolffd@0 140 if length(dd)>1
wolffd@0 141 anyframe = 1;
wolffd@0 142 end
wolffd@0 143 dn = cell2array(dd);
wolffd@0 144 [m{i}{k},s{i}{k},sl{i}{k},pa{i}{k},pf{i}{k},pe{i}{k}] ...
wolffd@0 145 = computestat(dn,fp{i}{k});
wolffd@0 146 elseif size(dd,2) < 2
wolffd@0 147 nonan = find(not(isnan(dd)));
wolffd@0 148 dn = dd(nonan);
wolffd@0 149 if isempty(dn)
wolffd@0 150 m{i}{k} = NaN;
wolffd@0 151 else
wolffd@0 152 m{i}{k} = mean(dn,2);
wolffd@0 153 end
wolffd@0 154 s{i}{k} = NaN;
wolffd@0 155 sl{i}{k} = NaN;
wolffd@0 156 pa{i}{k} = NaN;
wolffd@0 157 pf{i}{k} = NaN;
wolffd@0 158 pe{i}{k} = NaN;
wolffd@0 159 else
wolffd@0 160 anyframe = 1;
wolffd@0 161 s1 = size(dd,1);
wolffd@0 162 s3 = size(dd,3);
wolffd@0 163 dd = mean(dd,4); %mean(dd,3),4);
wolffd@0 164 m{i}{k} = NaN(s1,1,s3);
wolffd@0 165 s{i}{k} = NaN(s1,1,s3);
wolffd@0 166 sl{i}{k} = NaN(s1,1,s3);
wolffd@0 167 pa{i}{k} = NaN(s1,1,s3);
wolffd@0 168 pf{i}{k} = NaN(s1,1,s3);
wolffd@0 169 pe{i}{k} = NaN(s1,1,s3);
wolffd@0 170 for h = 1:s3
wolffd@0 171 for j = 1:s1
wolffd@0 172 [m{i}{k}(j,1,h),s{i}{k}(j,1,h),sl{i}{k}(j,1,h),...
wolffd@0 173 pa{i}{k}(j,1,h),pf{i}{k}(j,1,h),pe{i}{k}(j,1,h)] = ...
wolffd@0 174 computestat(dd(j,:,h),fp{i}{k});
wolffd@0 175 end
wolffd@0 176 end
wolffd@0 177 end
wolffd@0 178 end
wolffd@0 179 end
wolffd@0 180 if anyframe
wolffd@0 181 fields = {'Mean','Std','Slope','PeriodFreq','PeriodAmp','PeriodEntropy'};
wolffd@0 182 stats = {m,s,sl,pf,pa,pe};
wolffd@0 183 else
wolffd@0 184 fields = {'Mean'};
wolffd@0 185 stats = {m};
wolffd@0 186 end
wolffd@0 187 for i = 1:length(stats)
wolffd@0 188 data = stats{i};
wolffd@0 189 data = uncell(data,NaN);
wolffd@0 190 stat.(strcat(field,fields{i})) = data;
wolffd@0 191 end
wolffd@0 192 else
wolffd@0 193 if iscell(d{1}{1})
wolffd@0 194 slash = strfind(filename{1},'/');
wolffd@0 195 nbfolders = 1;
wolffd@0 196 infolder = 0;
wolffd@0 197 foldername{1} = filename{1}(1:slash(end)-1);
wolffd@0 198 for i = 1:length(d)
wolffd@0 199 slash = strfind(filename{i},'/');
wolffd@0 200 if not(strcmpi(filename{i}(1:slash(end)-1),foldername))
wolffd@0 201 nbfolders = nbfolders + 1;
wolffd@0 202 infolder = 0;
wolffd@0 203 foldername{nbfolders} = filename{i}(1:slash(end)-1);
wolffd@0 204 end
wolffd@0 205 infolder = infolder+1;
wolffd@0 206 dd{nbfolders}(infolder,:) = cell2array(d{i}{1});
wolffd@0 207 end
wolffd@0 208 for i = 1:length(dd)
wolffd@0 209 figure
wolffd@0 210 plot(dd{i}')
wolffd@0 211 stat.Mean{i} = mean(dd{i});
wolffd@0 212 stat.Std{i} = std(dd{i});
wolffd@0 213 end
wolffd@0 214 end
wolffd@0 215 end
wolffd@0 216 if not(isempty(labels) || isempty(labels{1}))
wolffd@0 217 stat.Class = labels;
wolffd@0 218 end
wolffd@0 219
wolffd@0 220
wolffd@0 221 function dn = cell2array(dd)
wolffd@0 222 dn = zeros(1,length(dd));
wolffd@0 223 for j = 1:length(dd)
wolffd@0 224 if isempty(dd{j})
wolffd@0 225 dn(j) = NaN;
wolffd@0 226 else
wolffd@0 227 dn(j) = dd{j}(1);
wolffd@0 228 end
wolffd@0 229 end
wolffd@0 230
wolffd@0 231
wolffd@0 232 function [m,s,sl,pa,pf,pe] = computestat(d,fp)
wolffd@0 233 m = NaN;
wolffd@0 234 s = NaN;
wolffd@0 235 sl = NaN;
wolffd@0 236 pa = NaN;
wolffd@0 237 pf = NaN;
wolffd@0 238 pe = NaN;
wolffd@0 239 diffp = fp(1,2:end) - fp(1,1:end-1);
wolffd@0 240 if isempty(diffp) || sum(round((diffp(2:end)-diffp(1:end-1))*1000))
wolffd@0 241 % Not regular sampling (in mirattacktime for instance)
wolffd@0 242 framesampling = NaN;
wolffd@0 243 else
wolffd@0 244 framesampling = fp(1,2)-fp(1,1);
wolffd@0 245 end
wolffd@0 246 nonan = find(not(isnan(d)));
wolffd@0 247 if not(isempty(nonan))
wolffd@0 248 dn = d(nonan);
wolffd@0 249 m = mean(dn,2);
wolffd@0 250 s = std(dn,0,2);
wolffd@0 251 if not(isnan(s))
wolffd@0 252 if s
wolffd@0 253 dk = (dn-m)/s;
wolffd@0 254 tk = linspace(0,1,size(d,2));
wolffd@0 255 sl = dk(:)'/tk(nonan);
wolffd@0 256 elseif size(s,2) == 1
wolffd@0 257 s = NaN;
wolffd@0 258 end
wolffd@0 259 end
wolffd@0 260 if length(dn)>1
wolffd@0 261 cor = xcorr(dn',dn','coeff');
wolffd@0 262 cor = cor(ceil(length(cor)/2):end);
wolffd@0 263 % let's zero the first descending slope of the
wolffd@0 264 % autocorrelation function
wolffd@0 265 firstmin = find(cor(2:end)>cor(1:end-1));
wolffd@0 266 if not(isempty(firstmin) || isnan(framesampling))
wolffd@0 267 cor2 = cor;
wolffd@0 268 cor2(1:firstmin(1)) = 0;
wolffd@0 269 [pa,pfk] = max(cor2);
wolffd@0 270 if pfk > 1
wolffd@0 271 pf = 1/(pfk-1)/framesampling;
wolffd@0 272 end
wolffd@0 273 end
wolffd@0 274 cor = abs(cor);
wolffd@0 275 cor = cor/sum(cor);
wolffd@0 276 pe = -sum(cor.*log(cor+1e-12))./log(length(cor));
wolffd@0 277 end
wolffd@0 278 end
wolffd@0 279
wolffd@0 280
wolffd@0 281 function b = isdesign(f)
wolffd@0 282 if iscell(f)
wolffd@0 283 f = f{1};
wolffd@0 284 end
wolffd@0 285 if isa(f,'mirdesign') || isa(f,'mirstruct')
wolffd@0 286 b = 1;
wolffd@0 287 elseif isa(f,'mirdata') || not(isstruct(f))
wolffd@0 288 b = 0;
wolffd@0 289 else
wolffd@0 290 fields = fieldnames(f);
wolffd@0 291 b = isdesign(f.(fields{1}));
wolffd@0 292 end