wolffd@0: function varargout = mirautocor(orig,varargin)
wolffd@0: %   a = mirautocor(x) computes the autocorrelation function related to x.
wolffd@0: %   Optional parameters:
wolffd@0: %       mirautocor(...,'Min',mi) indicates the lowest delay taken into
wolffd@0: %           consideration. The unit can be precised:
wolffd@0: %               mirautocor(...,'Min',mi,'s') (default unit)
wolffd@0: %               mirautocor(...,'Min',mi,'Hz')
wolffd@0: %           Default value: 0 s.
wolffd@0: %       mirautocor(...,'Max',ma) indicates the highest delay taken into
wolffd@0: %           consideration. The unit can be specified as for 'Min'.
wolffd@0: %           Default value:
wolffd@0: %               if x is a signal, the highest delay is 0.05 s
wolffd@0: %                   (corresponding to a minimum frequency of 20 Hz).
wolffd@0: %               if x is an envelope, the highest delay is 2 s.
wolffd@0: %       mirautocor(...,'Resonance',r) multiplies the autocorrelation function
wolffd@0: %           with a resonance curve:
wolffd@0: %           Possible values:
wolffd@0: %               'Toiviainen' from (Toiviainen & Snyder, 2003)
wolffd@0: %               'vanNoorden' from (van Noorden & Moelants, 2001)
wolffd@0: %           mirautocor(...,'Center',c) assigns the center value of the
wolffd@0: %               resonance curve, in seconds.
wolffd@0: %               Works mainly with 'Toiviainen' option.
wolffd@0: %               Default value: c = 0.5
wolffd@0: %       mirautocor(...,'Enhanced',a) reduces the effect of subharmonics.
wolffd@0: %           The original autocorrelation function is half-wave rectified,
wolffd@0: %           time-scaled by the factor a (which can be a factor list as
wolffd@0: %           well), and substracted from the original clipped function.
wolffd@0: %           (Tolonen & Karjalainen)
wolffd@0: %               If the 'Enhanced' option is not followed by any value, 
wolffd@0: %                   default value is a = 2:10
wolffd@0: %       mirautocor(...,'Halfwave') performs a half-wave rectification on the
wolffd@0: %           result.
wolffd@0: %       mirautocor(...,'Freq') represents the autocorrelation function in the
wolffd@0: %           frequency domain.
wolffd@0: %       mirautocor(...,'NormalWindow',w): applies a window to the input 
wolffd@0: %           signal and divides the autocorrelation by the autocorrelation of  
wolffd@0: %           that window (Boersma 1993).
wolffd@0: %           Possible values: any windowing function proposed in the Signal
wolffd@0: %               Processing Toolbox (help window) plus 'rectangle' (no
wolffd@0: %               windowing)
wolffd@0: %           Default value:  w = 'hanning'
wolffd@0: %           mirautocor(...,'NormalWindow',0): toggles off this normalization
wolffd@0: %               (which is on by default).
wolffd@0: %   All the parameters described previously can be applied to an
wolffd@0: %       autocorrelation function itself, in order to arrange the results
wolffd@0: %       after the actual computation of the autocorrelation computations.
wolffd@0: %       For instance: a = mirautocor(a,'Resonance','Enhanced')
wolffd@0: %   Other optional parameter:
wolffd@0: %       mirautocor(...,'Compres',k) computes the autocorrelation in the
wolffd@0: %           frequency domain and includes a magnitude compression of the
wolffd@0: %           spectral representation. A normal autocorrelation corresponds
wolffd@0: %           to the value k=2, but values lower than 2 are suggested by
wolffd@0: %           (Tolonen & Karjalainen, 2000).
wolffd@0: %           Default value: k = 0.67
wolffd@0: %       mirautocor(...,'Normal',n) or simply mirautocor(...,n) specifies
wolffd@0: %           the normalization strategy. Accepted values are 'biased',
wolffd@0: %           'unbiased', 'coeff' (default  value) and 'none'.
wolffd@0: %           See help xcorr for an explanation. 
wolffd@0: 
wolffd@0:         min.key = 'Min';
wolffd@0:         min.type = 'Integer';
wolffd@0:         min.unit = {'s','Hz'};
wolffd@0:         if isamir(orig,'mirspectrum')
wolffd@0:             min.defaultunit = 'Hz';
wolffd@0:         else
wolffd@0:             min.defaultunit = 's';
wolffd@0:         end
wolffd@0:         min.default = 0;
wolffd@0:         min.opposite = 'max';
wolffd@0:     option.min = min;
wolffd@0:         
wolffd@0:         max.key = 'Max';
wolffd@0:         max.type = 'Integer';
wolffd@0:         max.unit = {'s','Hz'};
wolffd@0:         if isamir(orig,'mirspectrum')
wolffd@0:             max.defaultunit = 'Hz';
wolffd@0:         else
wolffd@0:             max.defaultunit = 's';
wolffd@0:         end
wolffd@0:         if isamir(orig,'mirenvelope') || isamir(orig,'mirdiffenvelope')
wolffd@0:             max.default = 2;             % for envelopes, longest period: 2 seconds.
wolffd@0:         elseif isamir(orig,'miraudio')  || ischar(orig)  % for audio signal,lowest frequency: 20 Hz.
wolffd@0:             max.default = 1/20; 
wolffd@0:         else
wolffd@0:             max.default = Inf;
wolffd@0:         end
wolffd@0:         max.opposite = 'min';
wolffd@0:     option.max = max;
wolffd@0:         
wolffd@0:         scaleoptbw.key = 'Normal'; %'Normal' keyword optional
wolffd@0:         scaleoptbw.key = 'Boolean';
wolffd@0:     option.scaleoptbw = scaleoptbw;
wolffd@0:         scaleopt.type = 'String';
wolffd@0:         scaleopt.choice = {'biased','unbiased','coeff','none'};
wolffd@0:         scaleopt.default = 'coeff';
wolffd@0:     option.scaleopt = scaleopt;
wolffd@0:             
wolffd@0:         gener.key = {'Generalized','Compres'};
wolffd@0:         gener.type = 'Integer';
wolffd@0:         gener.default = 2;
wolffd@0:         gener.keydefault = .67;
wolffd@0:     option.gener = gener;
wolffd@0:         
wolffd@0:         ni.key = 'NormalInput';  %% Normalize before frame or chunk??
wolffd@0:         ni.type = 'Boolean';
wolffd@0:         ni.default = 0;
wolffd@0:     option.ni = ni;
wolffd@0:     
wolffd@0:         reso.key = 'Resonance';
wolffd@0:         reso.type = 'String';
wolffd@0:         reso.choice = {'ToiviainenSnyder','Toiviainen','vanNoorden','no','off',0};
wolffd@0:         reso.keydefault = 'Toiviainen';
wolffd@0:         reso.when = 'After';
wolffd@0:         reso.default = 0;
wolffd@0:     option.reso = reso;
wolffd@0:         
wolffd@0:         resocenter.key = {'Center','Centre'};
wolffd@0:         resocenter.type = 'Integer';
wolffd@0:         resocenter.when = 'After';
wolffd@0:     option.resocenter = resocenter;
wolffd@0: 
wolffd@0:         h.key = 'Halfwave';
wolffd@0:         h.type = 'Boolean';
wolffd@0:         h.when = 'After';
wolffd@0:         h.default = 0;
wolffd@0:     option.h = h;
wolffd@0:         
wolffd@0:         e.key = 'Enhanced';
wolffd@0:         e.type = 'Integers';
wolffd@0:         e.default = [];
wolffd@0:         e.keydefault = 2:10;
wolffd@0:         e.when = 'After';
wolffd@0:     option.e = e;
wolffd@0:         
wolffd@0:         fr.key = 'Freq';
wolffd@0:         fr.type = 'Boolean';
wolffd@0:         fr.default = 0;
wolffd@0:         fr.when = 'After';
wolffd@0:     option.fr = fr;
wolffd@0:         
wolffd@0:         nw.key = 'NormalWindow';
wolffd@0:         nw.when = 'Both';
wolffd@0:         if isamir(orig,'mirspectrum')
wolffd@0:             nw.default = 0;
wolffd@0:         elseif isamir(orig,'mirenvelope')
wolffd@0:             nw.default = 'rectangular';
wolffd@0:         else
wolffd@0:             nw.default = 'hanning';
wolffd@0:         end
wolffd@0:     option.nw = nw;
wolffd@0:     
wolffd@0:         win.key = 'Window';
wolffd@0:         win.type = 'String';
wolffd@0:         win.default = NaN;
wolffd@0:     option.win = win;
wolffd@0:     
wolffd@0: specif.option = option;
wolffd@0: 
wolffd@0: specif.defaultframelength = 0.05;
wolffd@0: specif.defaultframehop = 0.5;
wolffd@0: specif.eachchunk = @eachchunk;
wolffd@0: specif.combinechunk = @combinechunk;
wolffd@0: 
wolffd@0: if isamir(orig,'mirscalar') || isamir(orig,'mirenvelope')
wolffd@0:     specif.nochunk = 1;
wolffd@0: end
wolffd@0: 
wolffd@0: varargout = mirfunction(@mirautocor,orig,varargin,nargout,specif,@init,@main);
wolffd@0: 
wolffd@0: 
wolffd@0: function [x type] = init(x,option)
wolffd@0: type = 'mirautocor';
wolffd@0: 
wolffd@0: 
wolffd@0: function a = main(orig,option,postoption)
wolffd@0: if iscell(orig)
wolffd@0:     orig = orig{1};
wolffd@0: end
wolffd@0: if isa(orig,'mirautocor')
wolffd@0:     a = orig;
wolffd@0:     if not(isempty(option)) && ...
wolffd@0:             (option.min || iscell(option.max) || option.max < Inf)
wolffd@0:         coeff = get(a,'Coeff');
wolffd@0:         delay = get(a,'Delay');
wolffd@0:         for h = 1:length(coeff)
wolffd@0:             if a.freq
wolffd@0:                 mi = 1/option.max;
wolffd@0:                 ma = 1/option.min;
wolffd@0:             else
wolffd@0:                 mi = option.min;
wolffd@0:                 ma = option.max;
wolffd@0:             end
wolffd@0:             for k = 1:length(coeff{h})
wolffd@0:                 range = find(and(delay{h}{k}(:,1,1) >= mi,...
wolffd@0:                                  delay{h}{k}(:,1,1) <= ma));
wolffd@0:                 coeff{h}{k} = coeff{h}{k}(range,:,:);
wolffd@0:                 delay{h}{k} = delay{h}{k}(range,:,:);
wolffd@0:             end
wolffd@0:         end
wolffd@0:         a = set(a,'Coeff',coeff,'Delay',delay);
wolffd@0:     end
wolffd@0:     if not(isempty(postoption)) && not(isequal(postoption,0))
wolffd@0:         a = post(a,postoption);
wolffd@0:     end
wolffd@0: elseif ischar(orig)
wolffd@0:     a = mirautocor(miraudio(orig),option,postoption);
wolffd@0: else
wolffd@0:     if nargin == 0
wolffd@0:         orig = [];
wolffd@0:     end
wolffd@0:     a.freq = 0;
wolffd@0:     a.ofspectrum = 0;
wolffd@0:     a.window = {};
wolffd@0:     a.normalwindow = 0;
wolffd@0:     a = class(a,'mirautocor',mirdata(orig));
wolffd@0:     a = purgedata(a);
wolffd@0:     a = set(a,'Ord','coefficients');
wolffd@0:     sig = get(orig,'Data');
wolffd@0:     if isa(orig,'mirspectrum')
wolffd@0:         a = set(a,'Title','Spectrum autocorrelation','OfSpectrum',1,...
wolffd@0:                   'Abs','frequency (Hz)');
wolffd@0:         pos = get(orig,'Pos');
wolffd@0:     else
wolffd@0:         if isa(orig,'mirscalar')
wolffd@0:             a = set(a,'Title',[get(orig,'Title') ' autocorrelation']);
wolffd@0:             pos = get(orig,'FramePos');
wolffd@0:             for k = 1:length(sig)
wolffd@0:                 for l = 1:length(sig{k})
wolffd@0:                     sig{k}{l} = sig{k}{l}';
wolffd@0:                     pos{k}{l} = pos{k}{l}(1,:,:)';
wolffd@0:                 end
wolffd@0:             end
wolffd@0:         else
wolffd@0:             if isa(orig,'mirenvelope')
wolffd@0:                 a = set(a,'Title','Envelope autocorrelation');
wolffd@0:             elseif not(isa(orig,'mirautocor'))
wolffd@0:                 a = set(a,'Title','Waveform autocorrelation');
wolffd@0:             end
wolffd@0:             pos = get(orig,'Pos');
wolffd@0:         end
wolffd@0:         a = set(a,'Abs','lag (s)');
wolffd@0:     end
wolffd@0:     f = get(orig,'Sampling');
wolffd@0:     
wolffd@0:     if isnan(option.win) 
wolffd@0:         if isequal(option.nw,0) || ...
wolffd@0:                 strcmpi(option.nw,'Off') || strcmpi(option.nw,'No')
wolffd@0:             option.win = 0;
wolffd@0:         elseif isequal(option.nw,1) || strcmpi(option.nw,'On') || ...
wolffd@0:                                      strcmpi(option.nw,'Yes')
wolffd@0:             option.win = 'hanning';
wolffd@0:         else
wolffd@0:             option.win = postoption.nw;
wolffd@0:         end
wolffd@0:     end
wolffd@0: 
wolffd@0:     coeff = cell(1,length(sig));
wolffd@0:     lags = cell(1,length(sig));
wolffd@0:     wind = cell(1,length(sig));
wolffd@0:     for k = 1:length(sig)
wolffd@0:         s = sig{k};
wolffd@0:         p = pos{k};
wolffd@0:         fk = f{k};
wolffd@0:         if iscell(option.max)
wolffd@0:             mi = option.min{k};
wolffd@0:             ma = option.max{k};
wolffd@0:         else
wolffd@0:             mi = option.min;
wolffd@0:             ma = option.max;
wolffd@0:         end
wolffd@0:         coeffk = cell(1,length(s));
wolffd@0:         lagsk = cell(1,length(s));
wolffd@0:         windk = cell(1,length(s));
wolffd@0:         for l = 1:length(s)
wolffd@0:             sl = s{l};
wolffd@0:             sl(isnan(sl)) = 0;
wolffd@0:             if option.ni
wolffd@0:                 mxsl = repmat(max(sl),[size(sl,1),1,1]);
wolffd@0:                 mnsl = repmat(min(sl),[size(sl,1),1,1]);
wolffd@0:                 sl = (sl-mnsl)./(mxsl-mnsl);
wolffd@0:             end
wolffd@0:             pl = p{l};
wolffd@0:             pl = pl-repmat(pl(1,:,:),[size(pl,1),1,1]);
wolffd@0:             ls = size(sl,1);
wolffd@0:             if mi
wolffd@0:                 misp = find(pl(:,1,1)>=mi);
wolffd@0:                 if isempty(misp)
wolffd@0:                     warning('WARNING IN MIRAUTOCOR: The specified range of delays exceeds the temporal length of the signal.');
wolffd@0:                     disp('Minimum delay set to zero.')
wolffd@0:                     misp = 1;  % misp is the lowest index of the lag range
wolffd@0:                     mi = 0;
wolffd@0:                 else
wolffd@0:                     misp = misp(1);
wolffd@0:                 end
wolffd@0:             else
wolffd@0:                 misp = 1;
wolffd@0:             end
wolffd@0:             if ma
wolffd@0:                 masp = find(pl(:,1,1)>=ma);
wolffd@0:                 if isempty(masp)
wolffd@0:                     masp = Inf;
wolffd@0:                 else
wolffd@0:                     masp = masp(1);
wolffd@0:                 end
wolffd@0:             else
wolffd@0:                 masp = Inf;
wolffd@0:             end
wolffd@0:             masp = min(masp,ceil(ls/2));
wolffd@0:             if masp <= misp
wolffd@0:                 if size(sl,2) > 1
wolffd@0:                     warning('WARNING IN MIRAUTOCOR: Frame length is too small.');    
wolffd@0:                 else
wolffd@0:                     warning('WARNING IN MIRAUTOCOR: The audio sequence is too small.');    
wolffd@0:                 end
wolffd@0:                 display('The autocorrelation is not defined for this range of delays.');
wolffd@0:             end
wolffd@0:             sl = center(sl);
wolffd@0:             if not(ischar(option.win)) || strcmpi(option.win,'Rectangular')
wolffd@0:                 kw = ones(size(sl));
wolffd@0:             else
wolffd@0:                 N = size(sl,1);
wolffd@0:                 winf = str2func(option.win);
wolffd@0:                 try
wolffd@0:                     w = window(winf,N);
wolffd@0:                 catch
wolffd@0:                     if strcmpi(option.win,'hamming')
wolffd@0:                         disp('Signal Processing Toolbox does not seem to be installed. Recompute the hamming window manually.');
wolffd@0:                         w = 0.54 - 0.46 * cos(2*pi*(0:N-1)'/(N-1));
wolffd@0:                     else
wolffd@0:                         warning(['WARNING in MIRAUTOCOR: Unknown windowing function ',option.win,' (maybe Signal Processing Toolbox is not installed).']);
wolffd@0:                         disp('No windowing performed.')
wolffd@0:                         w = ones(size(sl,1),1);
wolffd@0:                     end
wolffd@0:                 end
wolffd@0:                 kw = repmat(w,[1,size(sl,2),size(sl,3)]);
wolffd@0:                 sl = sl.* kw;
wolffd@0:             end
wolffd@0: 
wolffd@0:             if strcmpi(option.scaleopt,'coeff')
wolffd@0:                 scaleopt = 'none';
wolffd@0:             else
wolffd@0:                 scaleopt = option.scaleopt;
wolffd@0:             end
wolffd@0:             c = zeros(masp,size(sl,2),size(sl,3));
wolffd@0:             for i = 1:size(sl,2)
wolffd@0:                 for j = 1:size(sl,3)
wolffd@0:                     if option.gener == 2
wolffd@0:                         cc = xcorr(sl(:,i,j),masp-1,scaleopt);
wolffd@0:                         c(:,i,j) = cc(masp:end);
wolffd@0:                     else
wolffd@0:                         ss = abs(fft(sl(:,i,j)));
wolffd@0:                         ss = ss.^option.gener;
wolffd@0:                         cc = ifft(ss);
wolffd@0:                         ll = (0:masp-1);
wolffd@0:                         c(:,i,j) = cc(ll+1);
wolffd@0:                     end
wolffd@0:                 end
wolffd@0:                 if strcmpi(option.scaleopt,'coeff') && option.gener == 2
wolffd@0:                     % to be adapted to generalized autocor
wolffd@0:                     c(:,i,:) = c(:,i,:)/xcorr(sum(sl(:,i,:),3),0);
wolffd@0:                     % This is a kind of generalization of the 'coeff'
wolffd@0:                     % normalization for multi-channels signals. In the
wolffd@0:                     % original 'coeff' option, the autocorrelation at zero
wolffd@0:                     % lag is identically 1.0. In this multi-channels
wolffd@0:                     % version, the autocorrelation at zero lag is such that
wolffd@0:                     % the sum over channels becomes identically 1.0. 
wolffd@0:                 end
wolffd@0:             end
wolffd@0:             coeffk{l} = c(misp:end,:,:);
wolffd@0:             pl = pl(find(pl(:,1,1) >=mi),:,:);
wolffd@0:             lagsk{l} = pl(1:min(size(coeffk{l},1),size(pl,1)),:,:);
wolffd@0:             windk{l} = kw;
wolffd@0:         end
wolffd@0:         coeff{k} = coeffk;
wolffd@0:         lags{k} = lagsk;
wolffd@0:         wind{k} = windk;
wolffd@0:     end
wolffd@0:     a = set(a,'Coeff',coeff,'Delay',lags,'Window',wind);
wolffd@0:     if not(isempty(postoption))
wolffd@0:         a = post(a,postoption);
wolffd@0:     end
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: function a = post(a,option)
wolffd@0: debug = 0;
wolffd@0: coeff = get(a,'Coeff');
wolffd@0: lags = get(a,'Delay');
wolffd@0: wind = get(a,'Window');
wolffd@0: freq = option.fr && not(get(a,'FreqDomain'));
wolffd@0: if isequal(option.e,1)
wolffd@0:     option.e = 2:10;
wolffd@0: end
wolffd@0: if max(option.e) > 1
wolffd@0:     pa = mirpeaks(a,'NoBegin','NoEnd','Contrast',.01);
wolffd@0:     va = mirpeaks(a,'Valleys','Contrast',.01);
wolffd@0:     pv = get(pa,'PeakVal');
wolffd@0:     vv = get(va,'PeakVal');
wolffd@0: end
wolffd@0: for k = 1:length(coeff)
wolffd@0:     for l = 1:length(coeff{k})
wolffd@0:         c = coeff{k}{l};  % Coefficients of autocorrelation
wolffd@0:         t = lags{k}{l};   % Delays of autocorrelation
wolffd@0:         if not(isempty(c))
wolffd@0:             if not(isequal(option.nw,0) || strcmpi(option.nw,'No') || ...
wolffd@0:                    strcmpi(option.nw,'Off') || a.normalwindow) % 'NormalWindow' option
wolffd@0:                 xw = zeros(size(c));
wolffd@0:                 lc = size(c,1);
wolffd@0:                 for j = 1:size(c,3)
wolffd@0:                     for i = 1:size(c,2)
wolffd@0:                         xwij = xcorr(wind{k}{l}(:,i,j),lc,'coeff');
wolffd@0:                         xw(:,i,j) = xwij(lc+2:end);
wolffd@0:                     end
wolffd@0:                 end
wolffd@0:                 c = c./ xw;
wolffd@0:                 a.normalwindow = 1;
wolffd@0:             end
wolffd@0:             if ischar(option.reso) && ...
wolffd@0:                     (strcmpi(option.reso,'ToiviainenSnyder') || ...
wolffd@0:                     strcmpi(option.reso,'Toiviainen') || ...
wolffd@0:                     strcmpi(option.reso,'vanNoorden'))
wolffd@0:                 if isa(a,'mirautocor') && get(a,'FreqDomain')
wolffd@0:                     ll = 1./t;
wolffd@0:                 else
wolffd@0:                     ll = t;
wolffd@0:                 end
wolffd@0:                 if not(option.resocenter)
wolffd@0:                     option.resocenter = .5;
wolffd@0:                 end
wolffd@0:                 if strcmpi(option.reso,'ToiviainenSnyder') || ...
wolffd@0:                     strcmpi(option.reso,'Toiviainen')
wolffd@0:                     w = max(1 - 0.25*(log2(max(ll,1e-12)/option.resocenter)).^2, 0);
wolffd@0:                 elseif strcmpi(option.reso,'vanNoorden')
wolffd@0:                     f0=2.193; b=option.resocenter; 
wolffd@0:                     f=1./ll; a1=(f0*f0-f.*f).^2+b*f.^2; a2=f0^4+f.^4;
wolffd@0:                     w=(1./sqrt(a1))-(1./sqrt(a2));
wolffd@0:                 end
wolffd@0:                 if max(w) == 0
wolffd@0:                     warning('The resonance curve, not defined for this range of delays, will not be applied.')
wolffd@0:                 else
wolffd@0:                     w = w/max(w);
wolffd@0:                     c = c.* repmat(w,[1,size(c,2),size(c,3)]);
wolffd@0:                 end
wolffd@0:             end
wolffd@0:             if option.h
wolffd@0:                 c = hwr(c);
wolffd@0:             end
wolffd@0:             if max(option.e) > 1
wolffd@0:                 if a.freq
wolffd@0:                     freq = 1;
wolffd@0:                     for i = 1:size(c,3)
wolffd@0:                         c(:,:,i) = flipud(c(:,:,i));
wolffd@0:                     end
wolffd@0:                     t = flipud(1./t);
wolffd@0:                 end
wolffd@0:                 
wolffd@0:                 for g = 1:size(c,2)
wolffd@0:                     for h = 1:size(c,3)
wolffd@0:                         cgh = c(:,g,h);
wolffd@0:                         if length(cgh)>1
wolffd@0:                             pvk = pv{k}{l}{1,g,h};
wolffd@0:                             mv = [];
wolffd@0:                             if not(isempty(pvk))
wolffd@0:                                 mp = min(pv{k}{l}{1,g,h}); %Lowest peak
wolffd@0:                                 vvv = vv{k}{l}{1,g,h}; %Valleys
wolffd@0:                                 mv = vvv(find(vvv<mp,1,'last'));
wolffd@0:                                     %Highest valley below the lowest peak
wolffd@0: 
wolffd@0:                                 if not(isempty(mv))
wolffd@0:                                     cgh = cgh-mv;
wolffd@0:                                 end
wolffd@0:                             end
wolffd@0:                             cgh2 = cgh;
wolffd@0:                             tgh2 = t(:,g,1);
wolffd@0:                             coef = cgh(2)-cgh(1); % initial slope of the autocor curve
wolffd@0:                             tcoef = tgh2(2)-tgh2(1);
wolffd@0:                             deter = 0;
wolffd@0:                             inter = 0;
wolffd@0: 
wolffd@0:                             repet = find(not(diff(tgh2)));  % Avoid bug if repeated x-values
wolffd@0:                             if repet
wolffd@0:                                 warning('WARNING in MIRAUTOCOR: Two successive samples have exactly same temporal position.');
wolffd@0:                                 tgh2(repet+1) = tgh2(repet)+1e-12;
wolffd@0:                             end
wolffd@0: 
wolffd@0:                             if coef < 0
wolffd@0:                                 % initial descending slope removed
wolffd@0:                                 deter = find(diff(cgh2)>0,1)-1;
wolffd@0:                                     % number of removed points
wolffd@0:                                 if isempty(deter)
wolffd@0:                                     deter = 0;
wolffd@0:                                 end
wolffd@0:                                 cgh2(1:deter) = [];
wolffd@0:                                 tgh2(1:deter) = [];
wolffd@0:                                 coef = cgh2(2)-cgh2(1);
wolffd@0:                             end
wolffd@0: 
wolffd@0:                             if coef > 0
wolffd@0:                                 % initial ascending slope prolonged to the left
wolffd@0:                                 % until it reaches the x-axis
wolffd@0:                                 while cgh2(1) > 0
wolffd@0:                                     coef = coef*1.1;
wolffd@0:                                         % the further to the left, ...
wolffd@0:                                         % the more ascending is the slope
wolffd@0:                                         % (not sure it always works, though...)
wolffd@0:                                     inter = inter+1;
wolffd@0:                                         % number of added points
wolffd@0:                                     cgh2 = [cgh2(1)-coef; cgh2];
wolffd@0:                                     tgh2 = [tgh2(1)-tcoef; tgh2];
wolffd@0:                                 end
wolffd@0:                                 cgh2(1) = 0;
wolffd@0:                             end
wolffd@0: 
wolffd@0:                             for i = option.e  % Enhancing procedure
wolffd@0:                                 % option.e is the list of scaling factors
wolffd@0:                                 % i is the scaling factor
wolffd@0:                                 if i
wolffd@0:                                     be = find(tgh2 & tgh2/i >= tgh2(1),1);
wolffd@0:                                         % starting point of the substraction
wolffd@0:                                         % on the X-axis
wolffd@0: 
wolffd@0:                                     if not(isempty(be))
wolffd@0:                                         ic = interp1(tgh2,cgh2,tgh2/i);
wolffd@0:                                             % The scaled autocorrelation
wolffd@0:                                         ic(1:be-1) = 0;
wolffd@0:                                         ic(find(isnan(ic))) = Inf;
wolffd@0:                                             % All the NaN values are changed
wolffd@0:                                             % into 0 in the resulting curve
wolffd@0:                                         ic = max(ic,0);
wolffd@0: 
wolffd@0:                                         if debug
wolffd@0:                                            hold off,plot(tgh2,cgh2)
wolffd@0:                                         end
wolffd@0: 
wolffd@0:                                         cgh2 = cgh2 - ic;       
wolffd@0:                                             % The scaled autocorrelation
wolffd@0:                                             % is substracted to the initial one
wolffd@0: 
wolffd@0:                                         cgh2 = max(cgh2,0);
wolffd@0:                                             % Half-wave rectification
wolffd@0: 
wolffd@0:                                         if debug
wolffd@0:                                            hold on,plot(tgh2,ic,'r')
wolffd@0:                                            hold on,plot(tgh2,cgh2,'g')
wolffd@0:                                            drawnow
wolffd@0:                                            figure
wolffd@0:                                         end
wolffd@0:                                     end
wolffd@0:                                 end
wolffd@0:                             end
wolffd@0: 
wolffd@0:                             % The  temporary modifications are
wolffd@0:                             % removed from the final curve
wolffd@0:                             if inter>=deter
wolffd@0:                                 c(:,g,h) = cgh2(inter-deter+1:end);
wolffd@0:                                 if not(isempty(mv))
wolffd@0:                                     c(:,g,h) = c(:,g,h) + mv;
wolffd@0:                                 end
wolffd@0:                             else
wolffd@0:                                 c(:,g,h) = [zeros(deter-inter,1);cgh2];
wolffd@0:                             end
wolffd@0:                         end
wolffd@0:                     end
wolffd@0:                 end
wolffd@0:             end
wolffd@0:             if freq
wolffd@0:                 if t(1,1) == 0
wolffd@0:                     c = c(2:end,:,:);
wolffd@0:                     t = t(2:end,:,:);
wolffd@0:                 end
wolffd@0:                 for i = 1:size(c,3)
wolffd@0:                     c(:,:,i) = flipud(c(:,:,i));
wolffd@0:                 end
wolffd@0:                 t = flipud(1./t);
wolffd@0:             end
wolffd@0:             coeff{k}{l} = c;
wolffd@0:             lags{k}{l} = t;
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: a = set(a,'Coeff',coeff,'Delay',lags,'Freq');
wolffd@0: if freq
wolffd@0:     a = set(a,'FreqDomain',1,'Abs','frequency (Hz)');
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: function [y orig] = eachchunk(orig,option,missing,postchunk)
wolffd@0: option.scaleopt = 'none';
wolffd@0: y = mirautocor(orig,option);
wolffd@0: 
wolffd@0: 
wolffd@0: function y = combinechunk(old,new)
wolffd@0: do = get(old,'Data');
wolffd@0: do = do{1}{1};
wolffd@0: dn = get(new,'Data');
wolffd@0: dn = dn{1}{1};
wolffd@0: if abs(size(dn,1)-size(do,1)) <= 2 % Probleme of border fluctuation
wolffd@0:     mi = min(size(dn,1),size(do,1));
wolffd@0:     dn = dn(1:mi,:,:);
wolffd@0:     do = do(1:mi,:,:);
wolffd@0: elseif length(dn) < length(do)
wolffd@0:     dn(length(do),:,:) = 0; % Zero-padding
wolffd@0: end
wolffd@0: y = set(old,'ChunkData',do+dn);