wolffd@0: function i = intersect(a,b,thr)
wolffd@0: 
wolffd@0: if nargin<3
wolffd@0:     thr = .01;
wolffd@0: end
wolffd@0: 
wolffd@0: va = get(a,'PeakVal');
wolffd@0: pa = get(a,'PeakPos');
wolffd@0: ma = get(a,'PeakMode');
wolffd@0: if isa(a,'mirscalar')
wolffd@0:     xa = get(a,'FramePos');
wolffd@0: else
wolffd@0:     xa = get(a,'Pos');
wolffd@0: end
wolffd@0: vb = get(b,'PeakVal');
wolffd@0: pb = get(b,'PeakPos');
wolffd@0: mb = get(b,'PeakMode');
wolffd@0: if isa(b,'mirscalar')
wolffd@0:     xb = get(b,'FramePos');
wolffd@0: else
wolffd@0:     xb = get(b,'Pos');
wolffd@0: end
wolffd@0: 
wolffd@0: for j = 1:length(va)
wolffd@0:     for k = 1:length(va{j})
wolffd@0:         [nl nc np] = size(va{j}{k});
wolffd@0:         for c = 1:nc
wolffd@0:             for p = 1:np
wolffd@0:                 [pacp,ix] = sort(pa{j}{k}{1,c,p});
wolffd@0:                 vacp = va{j}{k}{1,c,p}(ix);
wolffd@0:                 macp = ma{j}{k}{1,c,p}(ix);
wolffd@0:                 [pbcp,ix] = sort(pb{j}{k}{1,c,p});
wolffd@0:                 xajk = xa{j}{k};
wolffd@0:                 xbjk = xb{j}{k};
wolffd@0:                 if isa(a,'mirscalar')
wolffd@0:                     xajk = mean(xajk);
wolffd@0:                 end
wolffd@0:                 if isa(b,'mirscalar')
wolffd@0:                     xbjk = mean(xbjk);
wolffd@0:                 end
wolffd@0:                 ia = 1;
wolffd@0:                 ib = 1;
wolffd@0:                 ii = 1;
wolffd@0:                 while ia <= length(pacp) && ib <= length(pbcp)
wolffd@0:                     if abs(xajk(pacp(ia))-xbjk(pbcp(ib)))<thr
wolffd@0:                         picp(ii) = pacp(ia);
wolffd@0:                         vicp(ii) = vacp(ia);
wolffd@0:                         micp(ii) = macp(ia);
wolffd@0:                         ia = ia+1;
wolffd@0:                         ib = ib+1;
wolffd@0:                         ii = ii+1;
wolffd@0:                     elseif xajk(pacp(ia))<xbjk(pbcp(ib))
wolffd@0:                         ia = ia+1;
wolffd@0:                     else
wolffd@0:                         ib = ib+1;
wolffd@0:                     end
wolffd@0:                 end
wolffd@0:             pa{j}{k}{1,c,p} = picp;
wolffd@0:             va{j}{k}{1,c,p} = vicp;
wolffd@0:             ma{j}{k}{1,c,p} = micp;
wolffd@0:             end
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: 
wolffd@0: i = set(a,'PeakPos',pa,'PeakVal',va,'PeakMode',ma);