Mercurial > hg > camir-aes2014

annotate toolboxes/MIRtoolbox1.3.2/MIRToolbox/@mirautocor/mirautocor.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 = mirautocor(orig,varargin)
wolffd@0 2 % a = mirautocor(x) computes the autocorrelation function related to x.
wolffd@0 3 % Optional parameters:
wolffd@0 4 % mirautocor(...,'Min',mi) indicates the lowest delay taken into
wolffd@0 5 % consideration. The unit can be precised:
wolffd@0 6 % mirautocor(...,'Min',mi,'s') (default unit)
wolffd@0 7 % mirautocor(...,'Min',mi,'Hz')
wolffd@0 8 % Default value: 0 s.
wolffd@0 9 % mirautocor(...,'Max',ma) indicates the highest delay taken into
wolffd@0 10 % consideration. The unit can be specified as for 'Min'.
wolffd@0 11 % Default value:
wolffd@0 12 % if x is a signal, the highest delay is 0.05 s
wolffd@0 13 % (corresponding to a minimum frequency of 20 Hz).
wolffd@0 14 % if x is an envelope, the highest delay is 2 s.
wolffd@0 15 % mirautocor(...,'Resonance',r) multiplies the autocorrelation function
wolffd@0 16 % with a resonance curve:
wolffd@0 17 % Possible values:
wolffd@0 18 % 'Toiviainen' from (Toiviainen & Snyder, 2003)
wolffd@0 19 % 'vanNoorden' from (van Noorden & Moelants, 2001)
wolffd@0 20 % mirautocor(...,'Center',c) assigns the center value of the
wolffd@0 21 % resonance curve, in seconds.
wolffd@0 22 % Works mainly with 'Toiviainen' option.
wolffd@0 23 % Default value: c = 0.5
wolffd@0 24 % mirautocor(...,'Enhanced',a) reduces the effect of subharmonics.
wolffd@0 25 % The original autocorrelation function is half-wave rectified,
wolffd@0 26 % time-scaled by the factor a (which can be a factor list as
wolffd@0 27 % well), and substracted from the original clipped function.
wolffd@0 28 % (Tolonen & Karjalainen)
wolffd@0 29 % If the 'Enhanced' option is not followed by any value,
wolffd@0 30 % default value is a = 2:10
wolffd@0 31 % mirautocor(...,'Halfwave') performs a half-wave rectification on the
wolffd@0 32 % result.
wolffd@0 33 % mirautocor(...,'Freq') represents the autocorrelation function in the
wolffd@0 34 % frequency domain.
wolffd@0 35 % mirautocor(...,'NormalWindow',w): applies a window to the input
wolffd@0 36 % signal and divides the autocorrelation by the autocorrelation of
wolffd@0 37 % that window (Boersma 1993).
wolffd@0 38 % Possible values: any windowing function proposed in the Signal
wolffd@0 39 % Processing Toolbox (help window) plus 'rectangle' (no
wolffd@0 40 % windowing)
wolffd@0 41 % Default value: w = 'hanning'
wolffd@0 42 % mirautocor(...,'NormalWindow',0): toggles off this normalization
wolffd@0 43 % (which is on by default).
wolffd@0 44 % All the parameters described previously can be applied to an
wolffd@0 45 % autocorrelation function itself, in order to arrange the results
wolffd@0 46 % after the actual computation of the autocorrelation computations.
wolffd@0 47 % For instance: a = mirautocor(a,'Resonance','Enhanced')
wolffd@0 48 % Other optional parameter:
wolffd@0 49 % mirautocor(...,'Compres',k) computes the autocorrelation in the
wolffd@0 50 % frequency domain and includes a magnitude compression of the
wolffd@0 51 % spectral representation. A normal autocorrelation corresponds
wolffd@0 52 % to the value k=2, but values lower than 2 are suggested by
wolffd@0 53 % (Tolonen & Karjalainen, 2000).
wolffd@0 54 % Default value: k = 0.67
wolffd@0 55 % mirautocor(...,'Normal',n) or simply mirautocor(...,n) specifies
wolffd@0 56 % the normalization strategy. Accepted values are 'biased',
wolffd@0 57 % 'unbiased', 'coeff' (default value) and 'none'.
wolffd@0 58 % See help xcorr for an explanation.
wolffd@0 59
wolffd@0 60 min.key = 'Min';
wolffd@0 61 min.type = 'Integer';
wolffd@0 62 min.unit = {'s','Hz'};
wolffd@0 63 if isamir(orig,'mirspectrum')
wolffd@0 64 min.defaultunit = 'Hz';
wolffd@0 65 else
wolffd@0 66 min.defaultunit = 's';
wolffd@0 67 end
wolffd@0 68 min.default = 0;
wolffd@0 69 min.opposite = 'max';
wolffd@0 70 option.min = min;
wolffd@0 71
wolffd@0 72 max.key = 'Max';
wolffd@0 73 max.type = 'Integer';
wolffd@0 74 max.unit = {'s','Hz'};
wolffd@0 75 if isamir(orig,'mirspectrum')
wolffd@0 76 max.defaultunit = 'Hz';
wolffd@0 77 else
wolffd@0 78 max.defaultunit = 's';
wolffd@0 79 end
wolffd@0 80 if isamir(orig,'mirenvelope') || isamir(orig,'mirdiffenvelope')
wolffd@0 81 max.default = 2; % for envelopes, longest period: 2 seconds.
wolffd@0 82 elseif isamir(orig,'miraudio') || ischar(orig) % for audio signal,lowest frequency: 20 Hz.
wolffd@0 83 max.default = 1/20;
wolffd@0 84 else
wolffd@0 85 max.default = Inf;
wolffd@0 86 end
wolffd@0 87 max.opposite = 'min';
wolffd@0 88 option.max = max;
wolffd@0 89
wolffd@0 90 scaleoptbw.key = 'Normal'; %'Normal' keyword optional
wolffd@0 91 scaleoptbw.key = 'Boolean';
wolffd@0 92 option.scaleoptbw = scaleoptbw;
wolffd@0 93 scaleopt.type = 'String';
wolffd@0 94 scaleopt.choice = {'biased','unbiased','coeff','none'};
wolffd@0 95 scaleopt.default = 'coeff';
wolffd@0 96 option.scaleopt = scaleopt;
wolffd@0 97
wolffd@0 98 gener.key = {'Generalized','Compres'};
wolffd@0 99 gener.type = 'Integer';
wolffd@0 100 gener.default = 2;
wolffd@0 101 gener.keydefault = .67;
wolffd@0 102 option.gener = gener;
wolffd@0 103
wolffd@0 104 ni.key = 'NormalInput'; %% Normalize before frame or chunk??
wolffd@0 105 ni.type = 'Boolean';
wolffd@0 106 ni.default = 0;
wolffd@0 107 option.ni = ni;
wolffd@0 108
wolffd@0 109 reso.key = 'Resonance';
wolffd@0 110 reso.type = 'String';
wolffd@0 111 reso.choice = {'ToiviainenSnyder','Toiviainen','vanNoorden','no','off',0};
wolffd@0 112 reso.keydefault = 'Toiviainen';
wolffd@0 113 reso.when = 'After';
wolffd@0 114 reso.default = 0;
wolffd@0 115 option.reso = reso;
wolffd@0 116
wolffd@0 117 resocenter.key = {'Center','Centre'};
wolffd@0 118 resocenter.type = 'Integer';
wolffd@0 119 resocenter.when = 'After';
wolffd@0 120 option.resocenter = resocenter;
wolffd@0 121
wolffd@0 122 h.key = 'Halfwave';
wolffd@0 123 h.type = 'Boolean';
wolffd@0 124 h.when = 'After';
wolffd@0 125 h.default = 0;
wolffd@0 126 option.h = h;
wolffd@0 127
wolffd@0 128 e.key = 'Enhanced';
wolffd@0 129 e.type = 'Integers';
wolffd@0 130 e.default = [];
wolffd@0 131 e.keydefault = 2:10;
wolffd@0 132 e.when = 'After';
wolffd@0 133 option.e = e;
wolffd@0 134
wolffd@0 135 fr.key = 'Freq';
wolffd@0 136 fr.type = 'Boolean';
wolffd@0 137 fr.default = 0;
wolffd@0 138 fr.when = 'After';
wolffd@0 139 option.fr = fr;
wolffd@0 140
wolffd@0 141 nw.key = 'NormalWindow';
wolffd@0 142 nw.when = 'Both';
wolffd@0 143 if isamir(orig,'mirspectrum')
wolffd@0 144 nw.default = 0;
wolffd@0 145 elseif isamir(orig,'mirenvelope')
wolffd@0 146 nw.default = 'rectangular';
wolffd@0 147 else
wolffd@0 148 nw.default = 'hanning';
wolffd@0 149 end
wolffd@0 150 option.nw = nw;
wolffd@0 151
wolffd@0 152 win.key = 'Window';
wolffd@0 153 win.type = 'String';
wolffd@0 154 win.default = NaN;
wolffd@0 155 option.win = win;
wolffd@0 156
wolffd@0 157 specif.option = option;
wolffd@0 158
wolffd@0 159 specif.defaultframelength = 0.05;
wolffd@0 160 specif.defaultframehop = 0.5;
wolffd@0 161 specif.eachchunk = @eachchunk;
wolffd@0 162 specif.combinechunk = @combinechunk;
wolffd@0 163
wolffd@0 164 if isamir(orig,'mirscalar') || isamir(orig,'mirenvelope')
wolffd@0 165 specif.nochunk = 1;
wolffd@0 166 end
wolffd@0 167
wolffd@0 168 varargout = mirfunction(@mirautocor,orig,varargin,nargout,specif,@init,@main);
wolffd@0 169
wolffd@0 170
wolffd@0 171 function [x type] = init(x,option)
wolffd@0 172 type = 'mirautocor';
wolffd@0 173
wolffd@0 174
wolffd@0 175 function a = main(orig,option,postoption)
wolffd@0 176 if iscell(orig)
wolffd@0 177 orig = orig{1};
wolffd@0 178 end
wolffd@0 179 if isa(orig,'mirautocor')
wolffd@0 180 a = orig;
wolffd@0 181 if not(isempty(option)) && ...
wolffd@0 182 (option.min || iscell(option.max) || option.max < Inf)
wolffd@0 183 coeff = get(a,'Coeff');
wolffd@0 184 delay = get(a,'Delay');
wolffd@0 185 for h = 1:length(coeff)
wolffd@0 186 if a.freq
wolffd@0 187 mi = 1/option.max;
wolffd@0 188 ma = 1/option.min;
wolffd@0 189 else
wolffd@0 190 mi = option.min;
wolffd@0 191 ma = option.max;
wolffd@0 192 end
wolffd@0 193 for k = 1:length(coeff{h})
wolffd@0 194 range = find(and(delay{h}{k}(:,1,1) >= mi,...
wolffd@0 195 delay{h}{k}(:,1,1) <= ma));
wolffd@0 196 coeff{h}{k} = coeff{h}{k}(range,:,:);
wolffd@0 197 delay{h}{k} = delay{h}{k}(range,:,:);
wolffd@0 198 end
wolffd@0 199 end
wolffd@0 200 a = set(a,'Coeff',coeff,'Delay',delay);
wolffd@0 201 end
wolffd@0 202 if not(isempty(postoption)) && not(isequal(postoption,0))
wolffd@0 203 a = post(a,postoption);
wolffd@0 204 end
wolffd@0 205 elseif ischar(orig)
wolffd@0 206 a = mirautocor(miraudio(orig),option,postoption);
wolffd@0 207 else
wolffd@0 208 if nargin == 0
wolffd@0 209 orig = [];
wolffd@0 210 end
wolffd@0 211 a.freq = 0;
wolffd@0 212 a.ofspectrum = 0;
wolffd@0 213 a.window = {};
wolffd@0 214 a.normalwindow = 0;
wolffd@0 215 a = class(a,'mirautocor',mirdata(orig));
wolffd@0 216 a = purgedata(a);
wolffd@0 217 a = set(a,'Ord','coefficients');
wolffd@0 218 sig = get(orig,'Data');
wolffd@0 219 if isa(orig,'mirspectrum')
wolffd@0 220 a = set(a,'Title','Spectrum autocorrelation','OfSpectrum',1,...
wolffd@0 221 'Abs','frequency (Hz)');
wolffd@0 222 pos = get(orig,'Pos');
wolffd@0 223 else
wolffd@0 224 if isa(orig,'mirscalar')
wolffd@0 225 a = set(a,'Title',[get(orig,'Title') ' autocorrelation']);
wolffd@0 226 pos = get(orig,'FramePos');
wolffd@0 227 for k = 1:length(sig)
wolffd@0 228 for l = 1:length(sig{k})
wolffd@0 229 sig{k}{l} = sig{k}{l}';
wolffd@0 230 pos{k}{l} = pos{k}{l}(1,:,:)';
wolffd@0 231 end
wolffd@0 232 end
wolffd@0 233 else
wolffd@0 234 if isa(orig,'mirenvelope')
wolffd@0 235 a = set(a,'Title','Envelope autocorrelation');
wolffd@0 236 elseif not(isa(orig,'mirautocor'))
wolffd@0 237 a = set(a,'Title','Waveform autocorrelation');
wolffd@0 238 end
wolffd@0 239 pos = get(orig,'Pos');
wolffd@0 240 end
wolffd@0 241 a = set(a,'Abs','lag (s)');
wolffd@0 242 end
wolffd@0 243 f = get(orig,'Sampling');
wolffd@0 244
wolffd@0 245 if isnan(option.win)
wolffd@0 246 if isequal(option.nw,0) || ...
wolffd@0 247 strcmpi(option.nw,'Off') || strcmpi(option.nw,'No')
wolffd@0 248 option.win = 0;
wolffd@0 249 elseif isequal(option.nw,1) || strcmpi(option.nw,'On') || ...
wolffd@0 250 strcmpi(option.nw,'Yes')
wolffd@0 251 option.win = 'hanning';
wolffd@0 252 else
wolffd@0 253 option.win = postoption.nw;
wolffd@0 254 end
wolffd@0 255 end
wolffd@0 256
wolffd@0 257 coeff = cell(1,length(sig));
wolffd@0 258 lags = cell(1,length(sig));
wolffd@0 259 wind = cell(1,length(sig));
wolffd@0 260 for k = 1:length(sig)
wolffd@0 261 s = sig{k};
wolffd@0 262 p = pos{k};
wolffd@0 263 fk = f{k};
wolffd@0 264 if iscell(option.max)
wolffd@0 265 mi = option.min{k};
wolffd@0 266 ma = option.max{k};
wolffd@0 267 else
wolffd@0 268 mi = option.min;
wolffd@0 269 ma = option.max;
wolffd@0 270 end
wolffd@0 271 coeffk = cell(1,length(s));
wolffd@0 272 lagsk = cell(1,length(s));
wolffd@0 273 windk = cell(1,length(s));
wolffd@0 274 for l = 1:length(s)
wolffd@0 275 sl = s{l};
wolffd@0 276 sl(isnan(sl)) = 0;
wolffd@0 277 if option.ni
wolffd@0 278 mxsl = repmat(max(sl),[size(sl,1),1,1]);
wolffd@0 279 mnsl = repmat(min(sl),[size(sl,1),1,1]);
wolffd@0 280 sl = (sl-mnsl)./(mxsl-mnsl);
wolffd@0 281 end
wolffd@0 282 pl = p{l};
wolffd@0 283 pl = pl-repmat(pl(1,:,:),[size(pl,1),1,1]);
wolffd@0 284 ls = size(sl,1);
wolffd@0 285 if mi
wolffd@0 286 misp = find(pl(:,1,1)>=mi);
wolffd@0 287 if isempty(misp)
wolffd@0 288 warning('WARNING IN MIRAUTOCOR: The specified range of delays exceeds the temporal length of the signal.');
wolffd@0 289 disp('Minimum delay set to zero.')
wolffd@0 290 misp = 1; % misp is the lowest index of the lag range
wolffd@0 291 mi = 0;
wolffd@0 292 else
wolffd@0 293 misp = misp(1);
wolffd@0 294 end
wolffd@0 295 else
wolffd@0 296 misp = 1;
wolffd@0 297 end
wolffd@0 298 if ma
wolffd@0 299 masp = find(pl(:,1,1)>=ma);
wolffd@0 300 if isempty(masp)
wolffd@0 301 masp = Inf;
wolffd@0 302 else
wolffd@0 303 masp = masp(1);
wolffd@0 304 end
wolffd@0 305 else
wolffd@0 306 masp = Inf;
wolffd@0 307 end
wolffd@0 308 masp = min(masp,ceil(ls/2));
wolffd@0 309 if masp <= misp
wolffd@0 310 if size(sl,2) > 1
wolffd@0 311 warning('WARNING IN MIRAUTOCOR: Frame length is too small.');
wolffd@0 312 else
wolffd@0 313 warning('WARNING IN MIRAUTOCOR: The audio sequence is too small.');
wolffd@0 314 end
wolffd@0 315 display('The autocorrelation is not defined for this range of delays.');
wolffd@0 316 end
wolffd@0 317 sl = center(sl);
wolffd@0 318 if not(ischar(option.win)) || strcmpi(option.win,'Rectangular')
wolffd@0 319 kw = ones(size(sl));
wolffd@0 320 else
wolffd@0 321 N = size(sl,1);
wolffd@0 322 winf = str2func(option.win);
wolffd@0 323 try
wolffd@0 324 w = window(winf,N);
wolffd@0 325 catch
wolffd@0 326 if strcmpi(option.win,'hamming')
wolffd@0 327 disp('Signal Processing Toolbox does not seem to be installed. Recompute the hamming window manually.');
wolffd@0 328 w = 0.54 - 0.46 * cos(2*pi*(0:N-1)'/(N-1));
wolffd@0 329 else
wolffd@0 330 warning(['WARNING in MIRAUTOCOR: Unknown windowing function ',option.win,' (maybe Signal Processing Toolbox is not installed).']);
wolffd@0 331 disp('No windowing performed.')
wolffd@0 332 w = ones(size(sl,1),1);
wolffd@0 333 end
wolffd@0 334 end
wolffd@0 335 kw = repmat(w,[1,size(sl,2),size(sl,3)]);
wolffd@0 336 sl = sl.* kw;
wolffd@0 337 end
wolffd@0 338
wolffd@0 339 if strcmpi(option.scaleopt,'coeff')
wolffd@0 340 scaleopt = 'none';
wolffd@0 341 else
wolffd@0 342 scaleopt = option.scaleopt;
wolffd@0 343 end
wolffd@0 344 c = zeros(masp,size(sl,2),size(sl,3));
wolffd@0 345 for i = 1:size(sl,2)
wolffd@0 346 for j = 1:size(sl,3)
wolffd@0 347 if option.gener == 2
wolffd@0 348 cc = xcorr(sl(:,i,j),masp-1,scaleopt);
wolffd@0 349 c(:,i,j) = cc(masp:end);
wolffd@0 350 else
wolffd@0 351 ss = abs(fft(sl(:,i,j)));
wolffd@0 352 ss = ss.^option.gener;
wolffd@0 353 cc = ifft(ss);
wolffd@0 354 ll = (0:masp-1);
wolffd@0 355 c(:,i,j) = cc(ll+1);
wolffd@0 356 end
wolffd@0 357 end
wolffd@0 358 if strcmpi(option.scaleopt,'coeff') && option.gener == 2
wolffd@0 359 % to be adapted to generalized autocor
wolffd@0 360 c(:,i,:) = c(:,i,:)/xcorr(sum(sl(:,i,:),3),0);
wolffd@0 361 % This is a kind of generalization of the 'coeff'
wolffd@0 362 % normalization for multi-channels signals. In the
wolffd@0 363 % original 'coeff' option, the autocorrelation at zero
wolffd@0 364 % lag is identically 1.0. In this multi-channels
wolffd@0 365 % version, the autocorrelation at zero lag is such that
wolffd@0 366 % the sum over channels becomes identically 1.0.
wolffd@0 367 end
wolffd@0 368 end
wolffd@0 369 coeffk{l} = c(misp:end,:,:);
wolffd@0 370 pl = pl(find(pl(:,1,1) >=mi),:,:);
wolffd@0 371 lagsk{l} = pl(1:min(size(coeffk{l},1),size(pl,1)),:,:);
wolffd@0 372 windk{l} = kw;
wolffd@0 373 end
wolffd@0 374 coeff{k} = coeffk;
wolffd@0 375 lags{k} = lagsk;
wolffd@0 376 wind{k} = windk;
wolffd@0 377 end
wolffd@0 378 a = set(a,'Coeff',coeff,'Delay',lags,'Window',wind);
wolffd@0 379 if not(isempty(postoption))
wolffd@0 380 a = post(a,postoption);
wolffd@0 381 end
wolffd@0 382 end
wolffd@0 383
wolffd@0 384
wolffd@0 385 function a = post(a,option)
wolffd@0 386 debug = 0;
wolffd@0 387 coeff = get(a,'Coeff');
wolffd@0 388 lags = get(a,'Delay');
wolffd@0 389 wind = get(a,'Window');
wolffd@0 390 freq = option.fr && not(get(a,'FreqDomain'));
wolffd@0 391 if isequal(option.e,1)
wolffd@0 392 option.e = 2:10;
wolffd@0 393 end
wolffd@0 394 if max(option.e) > 1
wolffd@0 395 pa = mirpeaks(a,'NoBegin','NoEnd','Contrast',.01);
wolffd@0 396 va = mirpeaks(a,'Valleys','Contrast',.01);
wolffd@0 397 pv = get(pa,'PeakVal');
wolffd@0 398 vv = get(va,'PeakVal');
wolffd@0 399 end
wolffd@0 400 for k = 1:length(coeff)
wolffd@0 401 for l = 1:length(coeff{k})
wolffd@0 402 c = coeff{k}{l}; % Coefficients of autocorrelation
wolffd@0 403 t = lags{k}{l}; % Delays of autocorrelation
wolffd@0 404 if not(isempty(c))
wolffd@0 405 if not(isequal(option.nw,0) || strcmpi(option.nw,'No') || ...
wolffd@0 406 strcmpi(option.nw,'Off') || a.normalwindow) % 'NormalWindow' option
wolffd@0 407 xw = zeros(size(c));
wolffd@0 408 lc = size(c,1);
wolffd@0 409 for j = 1:size(c,3)
wolffd@0 410 for i = 1:size(c,2)
wolffd@0 411 xwij = xcorr(wind{k}{l}(:,i,j),lc,'coeff');
wolffd@0 412 xw(:,i,j) = xwij(lc+2:end);
wolffd@0 413 end
wolffd@0 414 end
wolffd@0 415 c = c./ xw;
wolffd@0 416 a.normalwindow = 1;
wolffd@0 417 end
wolffd@0 418 if ischar(option.reso) && ...
wolffd@0 419 (strcmpi(option.reso,'ToiviainenSnyder') || ...
wolffd@0 420 strcmpi(option.reso,'Toiviainen') || ...
wolffd@0 421 strcmpi(option.reso,'vanNoorden'))
wolffd@0 422 if isa(a,'mirautocor') && get(a,'FreqDomain')
wolffd@0 423 ll = 1./t;
wolffd@0 424 else
wolffd@0 425 ll = t;
wolffd@0 426 end
wolffd@0 427 if not(option.resocenter)
wolffd@0 428 option.resocenter = .5;
wolffd@0 429 end
wolffd@0 430 if strcmpi(option.reso,'ToiviainenSnyder') || ...
wolffd@0 431 strcmpi(option.reso,'Toiviainen')
wolffd@0 432 w = max(1 - 0.25*(log2(max(ll,1e-12)/option.resocenter)).^2, 0);
wolffd@0 433 elseif strcmpi(option.reso,'vanNoorden')
wolffd@0 434 f0=2.193; b=option.resocenter;
wolffd@0 435 f=1./ll; a1=(f0*f0-f.*f).^2+b*f.^2; a2=f0^4+f.^4;
wolffd@0 436 w=(1./sqrt(a1))-(1./sqrt(a2));
wolffd@0 437 end
wolffd@0 438 if max(w) == 0
wolffd@0 439 warning('The resonance curve, not defined for this range of delays, will not be applied.')
wolffd@0 440 else
wolffd@0 441 w = w/max(w);
wolffd@0 442 c = c.* repmat(w,[1,size(c,2),size(c,3)]);
wolffd@0 443 end
wolffd@0 444 end
wolffd@0 445 if option.h
wolffd@0 446 c = hwr(c);
wolffd@0 447 end
wolffd@0 448 if max(option.e) > 1
wolffd@0 449 if a.freq
wolffd@0 450 freq = 1;
wolffd@0 451 for i = 1:size(c,3)
wolffd@0 452 c(:,:,i) = flipud(c(:,:,i));
wolffd@0 453 end
wolffd@0 454 t = flipud(1./t);
wolffd@0 455 end
wolffd@0 456
wolffd@0 457 for g = 1:size(c,2)
wolffd@0 458 for h = 1:size(c,3)
wolffd@0 459 cgh = c(:,g,h);
wolffd@0 460 if length(cgh)>1
wolffd@0 461 pvk = pv{k}{l}{1,g,h};
wolffd@0 462 mv = [];
wolffd@0 463 if not(isempty(pvk))
wolffd@0 464 mp = min(pv{k}{l}{1,g,h}); %Lowest peak
wolffd@0 465 vvv = vv{k}{l}{1,g,h}; %Valleys
wolffd@0 466 mv = vvv(find(vvv<mp,1,'last'));
wolffd@0 467 %Highest valley below the lowest peak
wolffd@0 468
wolffd@0 469 if not(isempty(mv))
wolffd@0 470 cgh = cgh-mv;
wolffd@0 471 end
wolffd@0 472 end
wolffd@0 473 cgh2 = cgh;
wolffd@0 474 tgh2 = t(:,g,1);
wolffd@0 475 coef = cgh(2)-cgh(1); % initial slope of the autocor curve
wolffd@0 476 tcoef = tgh2(2)-tgh2(1);
wolffd@0 477 deter = 0;
wolffd@0 478 inter = 0;
wolffd@0 479
wolffd@0 480 repet = find(not(diff(tgh2))); % Avoid bug if repeated x-values
wolffd@0 481 if repet
wolffd@0 482 warning('WARNING in MIRAUTOCOR: Two successive samples have exactly same temporal position.');
wolffd@0 483 tgh2(repet+1) = tgh2(repet)+1e-12;
wolffd@0 484 end
wolffd@0 485
wolffd@0 486 if coef < 0
wolffd@0 487 % initial descending slope removed
wolffd@0 488 deter = find(diff(cgh2)>0,1)-1;
wolffd@0 489 % number of removed points
wolffd@0 490 if isempty(deter)
wolffd@0 491 deter = 0;
wolffd@0 492 end
wolffd@0 493 cgh2(1:deter) = [];
wolffd@0 494 tgh2(1:deter) = [];
wolffd@0 495 coef = cgh2(2)-cgh2(1);
wolffd@0 496 end
wolffd@0 497
wolffd@0 498 if coef > 0
wolffd@0 499 % initial ascending slope prolonged to the left
wolffd@0 500 % until it reaches the x-axis
wolffd@0 501 while cgh2(1) > 0
wolffd@0 502 coef = coef*1.1;
wolffd@0 503 % the further to the left, ...
wolffd@0 504 % the more ascending is the slope
wolffd@0 505 % (not sure it always works, though...)
wolffd@0 506 inter = inter+1;
wolffd@0 507 % number of added points
wolffd@0 508 cgh2 = [cgh2(1)-coef; cgh2];
wolffd@0 509 tgh2 = [tgh2(1)-tcoef; tgh2];
wolffd@0 510 end
wolffd@0 511 cgh2(1) = 0;
wolffd@0 512 end
wolffd@0 513
wolffd@0 514 for i = option.e % Enhancing procedure
wolffd@0 515 % option.e is the list of scaling factors
wolffd@0 516 % i is the scaling factor
wolffd@0 517 if i
wolffd@0 518 be = find(tgh2 & tgh2/i >= tgh2(1),1);
wolffd@0 519 % starting point of the substraction
wolffd@0 520 % on the X-axis
wolffd@0 521
wolffd@0 522 if not(isempty(be))
wolffd@0 523 ic = interp1(tgh2,cgh2,tgh2/i);
wolffd@0 524 % The scaled autocorrelation
wolffd@0 525 ic(1:be-1) = 0;
wolffd@0 526 ic(find(isnan(ic))) = Inf;
wolffd@0 527 % All the NaN values are changed
wolffd@0 528 % into 0 in the resulting curve
wolffd@0 529 ic = max(ic,0);
wolffd@0 530
wolffd@0 531 if debug
wolffd@0 532 hold off,plot(tgh2,cgh2)
wolffd@0 533 end
wolffd@0 534
wolffd@0 535 cgh2 = cgh2 - ic;
wolffd@0 536 % The scaled autocorrelation
wolffd@0 537 % is substracted to the initial one
wolffd@0 538
wolffd@0 539 cgh2 = max(cgh2,0);
wolffd@0 540 % Half-wave rectification
wolffd@0 541
wolffd@0 542 if debug
wolffd@0 543 hold on,plot(tgh2,ic,'r')
wolffd@0 544 hold on,plot(tgh2,cgh2,'g')
wolffd@0 545 drawnow
wolffd@0 546 figure
wolffd@0 547 end
wolffd@0 548 end
wolffd@0 549 end
wolffd@0 550 end
wolffd@0 551
wolffd@0 552 % The temporary modifications are
wolffd@0 553 % removed from the final curve
wolffd@0 554 if inter>=deter
wolffd@0 555 c(:,g,h) = cgh2(inter-deter+1:end);
wolffd@0 556 if not(isempty(mv))
wolffd@0 557 c(:,g,h) = c(:,g,h) + mv;
wolffd@0 558 end
wolffd@0 559 else
wolffd@0 560 c(:,g,h) = [zeros(deter-inter,1);cgh2];
wolffd@0 561 end
wolffd@0 562 end
wolffd@0 563 end
wolffd@0 564 end
wolffd@0 565 end
wolffd@0 566 if freq
wolffd@0 567 if t(1,1) == 0
wolffd@0 568 c = c(2:end,:,:);
wolffd@0 569 t = t(2:end,:,:);
wolffd@0 570 end
wolffd@0 571 for i = 1:size(c,3)
wolffd@0 572 c(:,:,i) = flipud(c(:,:,i));
wolffd@0 573 end
wolffd@0 574 t = flipud(1./t);
wolffd@0 575 end
wolffd@0 576 coeff{k}{l} = c;
wolffd@0 577 lags{k}{l} = t;
wolffd@0 578 end
wolffd@0 579 end
wolffd@0 580 end
wolffd@0 581 a = set(a,'Coeff',coeff,'Delay',lags,'Freq');
wolffd@0 582 if freq
wolffd@0 583 a = set(a,'FreqDomain',1,'Abs','frequency (Hz)');
wolffd@0 584 end
wolffd@0 585
wolffd@0 586
wolffd@0 587 function [y orig] = eachchunk(orig,option,missing,postchunk)
wolffd@0 588 option.scaleopt = 'none';
wolffd@0 589 y = mirautocor(orig,option);
wolffd@0 590
wolffd@0 591
wolffd@0 592 function y = combinechunk(old,new)
wolffd@0 593 do = get(old,'Data');
wolffd@0 594 do = do{1}{1};
wolffd@0 595 dn = get(new,'Data');
wolffd@0 596 dn = dn{1}{1};
wolffd@0 597 if abs(size(dn,1)-size(do,1)) <= 2 % Probleme of border fluctuation
wolffd@0 598 mi = min(size(dn,1),size(do,1));
wolffd@0 599 dn = dn(1:mi,:,:);
wolffd@0 600 do = do(1:mi,:,:);
wolffd@0 601 elseif length(dn) < length(do)
wolffd@0 602 dn(length(do),:,:) = 0; % Zero-padding
wolffd@0 603 end
wolffd@0 604 y = set(old,'ChunkData',do+dn);