


ALIGNWAVEFORMS Finds best alignment of waveforms to a reference. OFFSETS = ALIGNWAVEFORMS(TEMPLATE,WAVEFORMS,MAXSHIFT), where TEMPLATE is an 1 x T vector and WAVEFORMS is an N x T matrix, returns an N x 1 vector OFFSETS giving a shift in the range [-MAXSHIFT,MAXSHIFT] for each row of WAVEFORMS. If OFFSETS(k) is X, linear interpolation of WAVEFORMS(k,:) on an abcissa shifted by X samples will minimize the mean squared error between WAVEFORMS(k,:) and TEMPLATE. If TEMPLATE is an 1 x T x D array and WAVEFORMS is an N x T x D array, ALIGNWAVEFORMS finds OFFSETS that minimize the mean squared error of simultaneously shifting all pages of WAVEFORMS for a given row by the same amount. OFFSETS = ALIGNWAVEFORMS(TEMPLATE,WAVEFORMS) uses MAXSHIFT = 1. [OFFSETS,ALIGNED] = ALIGNWAVEFORMS(TEMPLATE,WAVEFORMS,...) also returns an N x (T-2*MAXSHIFT) x D matrix ALIGNED in which the rows of WAVEFORMS have been shifted using linear interpolation by the amounts given in OFFSETS. Note that ALIGNED has 2*MAXSHIFT fewer columns than WAVEFORMS since the interpolation operation invalidates edge samples. The algorithm used in ALIGNWAVEFORMS actually shifts the TEMPLATE waveform rather than the WAVEFORMS. This significantly speeds up the estimation but can provide erroneous results if TEMPLATE and WAVEFORMS do not oversample the underlying signal.


0001 function [offsets, aligned] = ... 0002 alignwaveforms(template, waveforms, mx) 0003 %ALIGNWAVEFORMS Finds best alignment of waveforms to a reference. 0004 % OFFSETS = ALIGNWAVEFORMS(TEMPLATE,WAVEFORMS,MAXSHIFT), where TEMPLATE 0005 % is an 1 x T vector and WAVEFORMS is an N x T matrix, returns an N x 1 0006 % vector OFFSETS giving a shift in the range [-MAXSHIFT,MAXSHIFT] for 0007 % each row of WAVEFORMS. If OFFSETS(k) is X, linear interpolation of 0008 % WAVEFORMS(k,:) on an abcissa shifted by X samples will minimize the 0009 % mean squared error between WAVEFORMS(k,:) and TEMPLATE. 0010 % 0011 % If TEMPLATE is an 1 x T x D array and WAVEFORMS is an N x T x D array, 0012 % ALIGNWAVEFORMS finds OFFSETS that minimize the mean squared error of 0013 % simultaneously shifting all pages of WAVEFORMS for a given row by the 0014 % same amount. 0015 % 0016 % OFFSETS = ALIGNWAVEFORMS(TEMPLATE,WAVEFORMS) uses MAXSHIFT = 1. 0017 % 0018 % [OFFSETS,ALIGNED] = ALIGNWAVEFORMS(TEMPLATE,WAVEFORMS,...) also 0019 % returns an N x (T-2*MAXSHIFT) x D matrix ALIGNED in which the rows of 0020 % WAVEFORMS have been shifted using linear interpolation by the amounts 0021 % given in OFFSETS. Note that ALIGNED has 2*MAXSHIFT fewer columns than 0022 % WAVEFORMS since the interpolation operation invalidates edge samples. 0023 % 0024 % The algorithm used in ALIGNWAVEFORMS actually shifts the TEMPLATE 0025 % waveform rather than the WAVEFORMS. This significantly speeds up the 0026 % estimation but can provide erroneous results if TEMPLATE and WAVEFORMS 0027 % do not oversample the underlying signal. 0028 0029 %%%%%%%%%%%%%%%%%%%%%%%%%%%% Parse Inputs %%%%%%%%%%%%%%%%%%%%%%%%%%%% 0030 if (nargin < 3), mx = 1; end; 0031 if (~isnumeric(mx) || mx <= 0 || mx~=round(mx)) 0032 error('MAXSHIFT must be a positive integer.'); 0033 end 0034 0035 [N,T,D] = size(waveforms); 0036 if (~isequal([1,T,D],size(template))) 0037 error('The second and third dimensions of TEMPLATE must match WAVEFORMS.'); 0038 end 0039 0040 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Setup %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0041 testoff = (-mx:1/3:mx)'; 0042 A = vander(testoff); A = A(:,end-2:end); 0043 Q = inv(A' * A) * A'; % Q fits a quadratic to errors at the different offsets 0044 0045 valid = (mx+1:T-mx); 0046 for off = 1:length(testoff) 0047 temp = permute(template, [2,3,1]); 0048 templateShift(off,:,:) = ... 0049 permute(interp1(1:T,temp,valid+testoff(off),'linear'), [3,1,2]); 0050 end 0051 0052 %%%%%%%%%%%%%%%%%%%%%%%%%% Find Alignments %%%%%%%%%%%%%%%%%%%%%%%%%% 0053 scores = zeros(length(testoff),N); offsets = zeros(N,1); 0054 for n = 1:N 0055 for off = 1:length(testoff) 0056 temp = templateShift(off,:,:); test = waveforms(n,valid,:); 0057 scores(off,n) = sum((temp(:) - test(:)).^2); 0058 end 0059 p = Q * scores(:,n); 0060 offsets(n) = -p(2)/(2*p(1)); 0061 end 0062 0063 offsets = -offsets; % change to waveforms shifts instead of template shifts 0064 offsets(offsets > mx) = mx; offsets(offsets < -mx) = -mx; % clip to maxshift 0065 0066 0067 %%%%%%%%%%%%%%%%%%%%%%%%% Realign waveforms %%%%%%%%%%%%%%%%%%%%%%%%%% 0068 if (nargout < 2), return; end; 0069 0070 [X,Y] = meshgrid(1:T, 1:N); 0071 aligned = zeros(N, T-2*mx, D); 0072 for d = 1:D 0073 Xoff = X(:,valid) + repmat(offsets, [1,T-2*mx]); 0074 aligned(:,:,d) = interp2(X,Y,waveforms(:,:,d), Xoff, Y(:,valid), 'linear'); 0075 end 0076 0077 0078 %%%%% The following code demonstrates the near-but-not-quite similar 0079 %%%%% effects of shifting a and comparing to b vs shifting b and comparing 0080 %%%%% to a. When the sampling rate of a and b is oversampled relative to 0081 %%%%% Nyquist for the underlying waveform, the difference is small. 0082 % D = 31; offsets = [-1:0.1:1]'; U = 100; 0083 % up = conv2(randn(1,D*U), gausskernel([0 2*U], U) , 'same'); 0084 % b = up(35:U:end); a = up(2:U:end); 0085 % Sab = zeros(length(offsets),1); Sba = zeros(length(offsets),1); 0086 % for k = 1:length(offsets) 0087 % at = interp1(1:D,a,[2:D-1]+offsets(k),'linear'); bt = b(2:end-1); 0088 % Sab(k) = sum([at(:)-bt(:)].^2); 0089 % 0090 % bt = interp1(1:D,b,[2:D-1]+offsets(k),'linear'); at = a(2:end-1); 0091 % Sba(k) = sum([at(:)-bt(:)].^2); 0092 % end 0093 % subplot(2,1,1); plot([a;b]'); 0094 % subplot(2,1,2); plot(offsets, [Sab flipud(Sba)]); 0095 0096 % Last Modified By: sbm on Wed Mar 1 21:18:50 2006 0097 0098 0099 % Last Modified By: sbm on Wed Mar 1 21:19:15 2006 0100 0101 0102 % Last Modified By: sbm on Wed Mar 1 21:19:46 2006 0103