Finding optimal waveform overlap in Matlab

Question

I'm trying to see some examples of how to find the optimal overlap between two waveforms. Here is some example data.

x1 = [108.1 108.2 108.3 108.4 108.5 108.6 108.7 108.8 108.9 109.0 109.1 109.2 109.3 109.4 109.5 109.6];
y1 = [0 0 2 6 7 6 2 -5 -6 -5 0 8 9 8 0 0];

x2 = [-1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20];
y2 = [0 0 0 3 6 9 8 7 6 5 3 -3 -7 -3 0 1 10 9 4 0 0 0];

Note the waveforms have different x-values and ranges for x-values. Specifically, I'd like to ONLY modify x2 and y2. The x values must maintain their relative spacing (ie, the distance between any two consecutive x-values (in the above x2 example, that's 1) must be the same, but in finding the optimal overlap, that new distance can be different from the original distance). In other words, the "shape" of the 2nd waveform must remain the same.

Generally, two steps need to be taken:

1. Assign new x2-values
2. Scale the y2-values by some global factor

I'd like to define optimal overlap as the minimized subtraction between the two waveforms. That is, from the first x-value (lowest of x1 and x2) to the last x-value (highest of x1 and x2), the subtraction of the two waveforms is a minimum. Note that the waveforms can be subtracted even if they have different x-values by interpolating between points. Where data is nonexistent, the subtraction should basically just include a 0 for one wave or the other.

Any ideas? Thanks in advance!

Answer 1

Method 1: Works well if overall shape is the same but there are no absolute features which are reliably comparable (like extrema).

function optim = minimize_distance()

x1 = [108.1 108.2 108.3 108.4 108.5 108.6 108.7 108.8 108.9 109.0 109.1 109.2 109.3 109.4 109.5 109.6];
y1 = [0 0 2 6 7 6 2 -5 -6 -5 0 8 9 8 0 0];
x2 = [-1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20];
y2 = [0 0 0 3 6 9 8 7 6 5 3 -3 -7 -3 0 1 10 9 4 0 0 0];

%Get educated guesses for parameters
xScalingInitial = abs((x1(end)-x1(1))/(x2(end)-x2(1)));
yScalingInitial = mean(abs(y2))/mean(abs(y1));
xOffsetInitial = (x2(1) - x1(1));

%determine how much these educated guesses can vary (here 10% for the
%Offset in x and 20% for the scaling)
x2diff = peak2peak(x2);
lowerBounds = [xScalingInitial - xScalingInitial*0.2, ...
               yScalingInitial - yScalingInitial*0.2, ...
               xOffsetInitial-x2diff*0.1];
upperBounds = [xScalingInitial + xScalingInitial*0.2, ...
               yScalingInitial + yScalingInitial*0.2, ...
               xOffsetInitial+x2diff*0.1];

%fminsearch performs an unconstrained search, fminsearchbnd (from
%MatlabFileExchange) also handles constrained functions. First input is the
%function to be minimized (see below, function to_minimize()). The output
%of this function (dist) will be minimized by fminsearch. @(x) tells matlab
%which variable is varied (see anonymous functions for reference). x1, x2,
%y1, y2 are parameters that are not changed by fminsearch. x is an array
%containing the three variable values for xScaling, yScaling and xOffSet
%which are varied. 
optim = fminsearchbnd(@(x) ...
                    to_minimize(x, x1, y1, x2, y2), ...
                    [xScalingInitial, yScalingInitial, xOffsetInitial],...
                    lowerBounds, upperBounds);

function dist = to_minimize(x, x1, y1, x2, y2)
%Assign variables from input array
xScaling = x(1);
yScaling = x(2);
xOffSet  = x(3);
%Get the scaled version of arrays y2 and x2.
y2Scaled = y2*yScaling;
x2Scaled = x2*xScaling-xOffSet;
%Linspace() creates 100 (default, this can be set to more or less if you
%want more or less precision) points between x1(1) and x1(end) (same for x2), 
%linearly spaced
x2Interp = linspace(x2Scaled(1), x2Scaled(end));
x1Interp = linspace(x1(1), x1(end));
%Interpolate y values at these points
y2Interp = interp1(x2Scaled, y2Scaled, x2Interp);
y1Interp = interp1(x1, y1, x1Interp);
%Now that we have two arrays of same size and scale we can compare by
%taking the point to point distance and minimizing it.
dist = sum((x1Interp-x2Interp).^2+(y1Interp-y2Interp).^2);
%Plot, comment out for performance!
clf
hold on;
plot(x1, y1);
plot(x2Interp, y2Interp);
hold off;
pause(0.01);

Output: Method 2: If the shape of the curves is virtually the same you are better off comparing the minima and maxima, and adjusting the indices/scaling of x2/y2 accordingly. Could be done as follows:

%find minima and maxima (and their indices) in original graphs
[~, minIdxY1] = min(y1);
[y1max, maxIdxY1] = max(y1);
[~, minIdxY2] = min(y2);
[y2max, maxIdxY2] = max(y2);

%Compare maxima, to get scaling factor for y
yScaling = y1max/y2max;
y2Scaled = y2*yScaling;

%Get x distance between minimum and maximum for both graphs
x1diff = abs(x1(minIdxY1) - x1(maxIdxY1));
x2Diff = abs(x2(minIdxY2) - x2(maxIdxY2));

%Stretch x2 to the shape of x1 by multiplying with the ratio of the two
%above distances
stretchFactor = x1diff/x2Diff;
x2stretched = x2* stretchFactor;

%Get new offset by comparing the midpoint between maximum and minimum
%of graph1 (x1/y1) and x2stretched/y2
midX1 = (x1(minIdxY1)+x1(maxIdxY1))/2;
midX2 = (x2stretched(minIdxY2)+x2stretched(maxIdxY2))/2;
x2stretchedOffset = x2stretched + (midX1-midX2);

figure;
hold on;
plot(x1, y1);
plot(x2stretchedOffset, y2Scaled);
legend('x1/y1', 'x2stretched/y2');
hold off;

Output: