Home > chronux_1_50 > spikesort > alignwaveforms.m

alignwaveforms

PURPOSE ^

ALIGNWAVEFORMS Finds best alignment of waveforms to a reference.

SYNOPSIS ^

function [offsets, aligned] =alignwaveforms(template, waveforms, mx)

DESCRIPTION ^

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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Mon 09-Oct-2006 00:54:52 by m2html © 2003