在 Matlab 中寻找最佳波形重叠

Finding optimal waveform overlap in Matlab

我正在尝试查看一些有关如何找到两个波形之间的最佳重叠的示例。这是一些示例数据。

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];

请注意,波形具有不同的 x 值和 x 值的范围。具体来说,我只想修改 x2y2x 值必须保持它们的相对间距(即任意两个连续 x 值之间的距离(在上面的 x2 示例中,即 1)必须相同,但在查找最佳重叠,新距离可以与原始距离不同)。 也就是说,第二个波形的"shape"必须保持不变。

一般需要进行两步:

  1. 分配新的 x2
  2. 通过一些全局因素缩放 y2

我想将最佳重叠定义为两个波形之间的最小减法。也就是说,从第一个 x 值(x1x2 中的最低值)到最后一个 x 值(x1 和 [=13 中的最高值) =]), 两个波形相减是最小值。请注意,即使波形具有不同的 x 值,也可以通过在点之间进行插值来减去波形。如果数据不存在,减法基本上应该只包括一个波或另一个波的 0

有什么想法吗?提前致谢!

方法 1: 如果整体形状相同但没有可靠比较的绝对特征(如极值),则效果很好。

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);

输出: 方法2:如果曲线的形状几乎相同,则最好比较最小值和最大值,并相应地调整x2/y2 的indices/scaling。可按如下方式进行:

%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;

输出: