V_TXALIGN Find the best alignment of two sets of time markers [KX,KY,N,M,S]=(X,Y,MAXT) x and y vectors contain a list of non-decreasing time values we find the alignment between them to minimize (x-y)^2 with a penalty of maxt^2 for every unmatched pair of entries from x and y. an x value can only match to the nearest y value in either the +ve or -ve direction If nsd is specified then the threshold is defined circularly to be nsd times the std deviation of the matches away from the mean
0001 function [kx,ky,nxy,mxy,sxy]=v_txalign(x,y,maxt,nsd) 0002 %V_TXALIGN Find the best alignment of two sets of time markers [KX,KY,N,M,S]=(X,Y,MAXT) 0003 % x and y vectors contain a list of non-decreasing time values 0004 % we find the alignment between them to minimize (x-y)^2 with a penalty of maxt^2 0005 % for every unmatched pair of entries from x and y. 0006 % an x value can only match to the nearest y value in either the +ve or -ve direction 0007 % If nsd is specified then the threshold is defined circularly to be nsd 0008 % times the std deviation of the matches away from the mean 0009 0010 % e.g. v_txalign([1 11 21 27 31 42 51],[2 12 15 22 25 32 41 52 61],1.1); 0011 0012 % 0013 % Outputs: kx is a column vector the same length as x. kx(i)=j if 0014 % x(i) is matched to y(j). kx(i)=0 if x(i) is unmatched. 0015 % ky is a column vector the same length as y. ky(j)=i if 0016 % y(j) is matched to x(i). kx(j)=0 if y(j) is unmatched. 0017 % nxy is the number of matched pairs 0018 % mxy is the mean of y(j)-x(i) for the matched pairs 0019 % sxy is the standard deviation (biassed) of y(j)-x(i) for matched pairs 0020 % 0021 % If no outputs are specified, v_txalign plots a graph showing y-x against mean(x,y) 0022 0023 % Copyright (C) Mike Brookes 2002 0024 % Version: $Id: v_txalign.m 10865 2018-09-21 17:22:45Z dmb $ 0025 % 0026 % VOICEBOX is a MATLAB toolbox for speech processing. 0027 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html% 0028 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0029 % This program is free software; you can redistribute it and/or modify 0030 % it under the terms of the GNU General Public License as published by 0031 % the Free Software Foundation; either version 2 of the License, or 0032 % (at your option) any later version. 0033 % 0034 % This program is distributed in the hope that it will be useful, 0035 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0036 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0037 % GNU General Public License for more details. 0038 % 0039 % You can obtain a copy of the GNU General Public License from 0040 % http://www.gnu.org/copyleft/gpl.html or by writing to 0041 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0042 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0043 0044 % d(*,i): i=1 cumulative cost if not the second of a pair 0045 % 2 cumulative cost if the second of a pair 0046 % 3 +1 if second of pair given that next is not second of pair 0047 % -1 if not allowed to be second of pair 0048 % 4 0 if this is an x value, 1 if a y value 0049 lx=length(x); 0050 ly=length(y); 0051 0052 if nargin<4 % no nsd specified 0053 x(lx+1)=2*abs(y(ly))+1; 0054 y(ly+1)=2*abs(x(lx))+1; 0055 lxy=lx+ly+1; 0056 d=zeros(lxy,4); 0057 0058 if lx>0 & ly>0 0059 c0=maxt^2; 0060 ix=1; 0061 iy=1; 0062 d(1,:)=[0 0 -1 y(1)<x(1)]; 0063 vp=0; 0064 for id=2:lxy 0065 if y(iy)<x(ix) % do y next 0066 v=y(iy); 0067 d(id,4)=1; 0068 iy=iy+1; 0069 else % do x next 0070 v=x(ix); 0071 ix=ix+1; 0072 end 0073 if d(id,4)==d(id-1,4) 0074 d(id,3)=-1; 0075 else 0076 d(id,2)=d(id-1,1)-c0+(v-vp)^2; 0077 end 0078 if ~d(id-1,3) & d(id-1,1)>=d(id-1,2) 0079 d(id,1)=d(id-1,2); 0080 d(id-1,3)=1; 0081 else 0082 d(id,1)=d(id-1,1); 0083 end 0084 vp=v; 0085 end 0086 if ~d(lxy,3) & d(lxy,1)>=d(lxy,2) 0087 d(lxy,3)=1; 0088 end 0089 0090 % now do the traceback 0091 0092 ix=lx; 0093 iy=ly; 0094 nxy=0; 0095 mxy=0; 0096 sxy=0; 0097 kx=zeros(lx,1); 0098 ky=zeros(ly,1); 0099 while (ix>0) & (iy>0) 0100 id=ix+iy+1; 0101 if d(id,3)>0 0102 ky(iy)=ix; % iy aligned with ix 0103 kx(ix)=iy; 0104 dt=y(iy)-x(ix); 0105 nxy=nxy+1; 0106 mxy=mxy+dt; 0107 sxy=sxy+dt^2; 0108 ix=ix-1; 0109 iy=iy-1; 0110 else 0111 ix = ix-1+d(id,4); 0112 iy = iy-d(id,4); 0113 end 0114 end 0115 if nxy 0116 mxy=mxy/nxy; 0117 sxy=sqrt(sxy/nxy-mxy^2); 0118 end 0119 % eliminate junk and plot results 0120 if ~nargout 0121 x(end)=[]; 0122 y(end)=[]; 0123 x=x(:); 0124 y=y(:); 0125 p=zeros(lxy-nxy-1,3); 0126 p=[x x -(kx==0); y(ky==0) y(ky==0) ones(ly-nxy,1)]; 0127 p(kx>0,2)=y(kx(kx>0)); 0128 p=sortrows(p); 0129 p(:,1:2)=p(:,1:2)*[0.5 -1; 0.5 1]; 0130 plot(p(p(:,3)==0,1),p(p(:,3)==0,2),'-*',p(p(:,3)<0,1),p(p(:,3)<0,2),'x',p(p(:,3)>0,1),p(p(:,3)>0,2),'o'); 0131 xlabel('Mean of x and y (x,o = unmatched x,y)'); 0132 ylabel('y-x difference'); 0133 end 0134 else % nsd specifies how many std deviations to exclude 0135 [kx,ky,nxy,mxy,sxy]=v_txalign(x,y,maxt); 0136 nxy0=nxy+1; 0137 while nxy<nxy0 0138 nxy0=nxy; 0139 mxy0=mxy; 0140 [kx,ky,nxy,mxy,sxy]=v_txalign(x+mxy0,y,nsd*sxy); 0141 end 0142 mxy=mxy+mxy0; 0143 if ~nargout 0144 v_txalign(x,y,nsd*sxy); % plot it 0145 end 0146 end 0147 else 0148 kx=zeros(lx,1); 0149 ky=zeros(ly,1); 0150 nxy=0; 0151 mxy=0; 0152 sxy=0; 0153 end