V_HORIZDIFF - Estimates the horizontal difference between two functions of x
Usage: z=horizdiff(y,v,x); % Approximately: y(x) = v(x+z)
Inputs: y(n,m) each column is a function of x
v(k,1) reference function of u
x(n,1) x values [default: x'=1:n]
u(k,1) x values for v [default: v=x]
q interpolation mode
'l' linear [default]
'p' pchip (n>=2)
's' spline (n>=4)
Outputs: z(n,m) horizontal difference between v and y
zm(1,m) MMSE horizontal difference that minimizes sum((y(x)-v(x+z)).^2)0001 function [z,zm]=v_horizdiff(y,v,x,u,q) 0002 %V_HORIZDIFF - Estimates the horizontal difference between two functions of x 0003 % 0004 % Usage: z=horizdiff(y,v,x); % Approximately: y(x) = v(x+z) 0005 % 0006 % Inputs: y(n,m) each column is a function of x 0007 % v(k,1) reference function of u 0008 % x(n,1) x values [default: x'=1:n] 0009 % u(k,1) x values for v [default: v=x] 0010 % q interpolation mode 0011 % 'l' linear [default] 0012 % 'p' pchip (n>=2) 0013 % 's' spline (n>=4) 0014 % 0015 % Outputs: z(n,m) horizontal difference between v and y 0016 % zm(1,m) MMSE horizontal difference that minimizes sum((y(x)-v(x+z)).^2) 0017 % 0018 0019 % Copyright (C) Mike Brookes 2012 0020 % Version: $Id: v_axisenlarge.m 10865 2018-09-21 17:22:45Z dmb $ 0021 % 0022 % VOICEBOX is a MATLAB toolbox for speech processing. 0023 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0024 % 0025 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0026 % This program is free software; you can redistribute it and/or modify 0027 % it under the terms of the GNU General Public License as published by 0028 % the Free Software Foundation; either version 2 of the License, or 0029 % (at your option) any later version. 0030 % 0031 % This program is distributed in the hope that it will be useful, 0032 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0033 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0034 % GNU General Public License for more details. 0035 % 0036 % You can obtain a copy of the GNU General Public License from 0037 % http://www.gnu.org/copyleft/gpl.html or by writing to 0038 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0039 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0040 % 0041 % 0042 % choose interpolation method 0043 % 0044 [n,m]=size(y); 0045 if nargin<5 || isempty(m) 0046 q=''; 0047 end 0048 if n>=4 && any(q=='s') 0049 im='spline'; 0050 elseif n>=2 && (any(q=='s') || any(q=='p')) 0051 im='pchip'; 0052 else 0053 im='linear'; 0054 end 0055 if nargin<3 || isempty(x) 0056 x=(1:n)'; 0057 end 0058 if nargin<4 || isempty(u) 0059 u=x; 0060 end 0061 % 0062 % Interpolate inverse function, u(v) 0063 % 0064 z=interp1(v,u,y,im,'extrap')-repmat(x,1,m); 0065 if nargout>1 0066 zm=zeros(1,m); 0067 ff=@(zm,c)sum((c-interp1(u,v,x+zm,im,'extrap')).^2); 0068 for i=1:m 0069 zm(i)=fminbnd(@(zm) ff(zm,y(:,i)),u(1)-x(end),u(end)-x(1)); 0070 end 0071 end 0072 if ~nargout 0073 subplot(212); 0074 plot(x,z); 0075 xlabel('x'); 0076 ylabel('x shift'); 0077 if all(abs(z-z(1))<=z(1)*1e-4) 0078 set(gca,'ylim',2*[min(z(1),0) max(z(1),0)]); 0079 end 0080 subplot(211); 0081 plot(x,y); 0082 hold on 0083 plot(u,v,':k'); 0084 hold off 0085 ylabel('y and v'); 0086 end