matlab - 在 Matlab 中沿绘图查找点

Question

我有以下图和创建该图的数据文件。我想让 Matlab 为我找到以下几点：

[y,x] 表示 100% 线指出的峰值
[x] 表示绘图与 y=0 线相交的位置
[x] 其中 y 是第 1 部分中发现的峰值的 50% 和 20%。

是否有任何人们知道的附加工具或软件包可以帮助我完成此任务？我需要为一组图执行此操作，因此合理自动化的东西将是理想的。

我当然可以在 Matlab 中进行编程和计算部分，只需能够加载数据文件，将其与曲线或函数匹配，并找到各种 [x,y] 坐标。

数据图

score 6 · Accepted Answer

好的，这就去。据我所知，Matlab 中没有内置的例程来做你想做的事。你必须自己做一个。需要注意的几点：

应该很明显，线性插值数据是最容易做到的，应该没有问题
使用单个多项式插值也不是太难，尽管还有一些细节需要处理。找到峰值应该是首要任务，这涉及找到导数的根（roots例如，使用）。找到峰值后，通过将多项式偏移该数量，在所有所需级别（0%、20%、50%）处找到多项式的根。
使用三次样条 ( spline) 是最复杂的。对于具有完整三次的所有子区间，应重复上面概述的一般多项式例程，考虑到最大值也可能位于子区间边界上的可能性，并且找到的任何根和极值都可能位于区间之外三次方有效（也不要忘记x使用的 -offsets spline）。

这是我对所有 3 种方法的实现：

%% Initialize
% ---------------------------

clc
clear all

% Create some bogus data
n = 25;

f = @(x) cos(x) .* sin(4*x/pi) + 0.5*rand(size(x));
x = sort( 2*pi * rand(n,1));
y = f(x);


%% Linear interpolation
% ---------------------------

% y peak
[y_peak, ind] = max(y);
x_peak = x(ind);

% y == 0%, 20%, 50%
lims = [0 20 50];
X = cell(size(lims));
Y = cell(size(lims));
for p = 1:numel(lims)

    % the current level line to solve for
    lim = y_peak*lims(p)/100;

    % points before and after passing through the current limit
    after    = (circshift(y<lim,1) & y>lim) | (circshift(y>lim,1) & y<lim);
    after(1) = false;
    before   = circshift(after,-1);

    xx = [x(before) x(after)];
    yy = [y(before) y(after)];

    % interpolate and insert new data
    new_X = x(before) - (y(before)-lim).*diff(xx,[],2)./diff(yy,[],2);
    X{p} = new_X;
    Y{p} = lim * ones(size(new_X));

end

% make a plot to verify
figure(1), clf, hold on
plot(x,y, 'r') % (this also plots the interpolation in this case)

plot(x_peak,y_peak, 'k.') % the peak

plot(X{1},Y{1}, 'r.') % the 0%  intersects
plot(X{2},Y{2}, 'g.') % the 20% intersects
plot(X{3},Y{3}, 'b.') % the 50% intersects

% finish plot
xlabel('X'), ylabel('Y'), title('Linear interpolation')
legend(...
    'Real data / interpolation',...
    'peak',...
    '0% intersects',...
    '20% intersects',...
    '50% intersects',...
    'location', 'southeast')



%% Cubic splines
% ---------------------------

% Find cubic splines interpolation
pp = spline(x,y);

% Finding the peak requires finding the maxima of all cubics in all
% intervals. This means evaluating the value of the interpolation on 
% the bounds of each interval, finding the roots of the derivative and
% evaluating the interpolation on those roots: 

coefs = pp.coefs;
derivCoefs = bsxfun(@times, [3 2 1], coefs(:,1:3));
LB = pp.breaks(1:end-1).'; % lower bounds of all intervals
UB = pp.breaks(2:end).';   % upper bounds of all intervals

% rename for clarity
a = derivCoefs(:,1);
b = derivCoefs(:,2);  
c = derivCoefs(:,3); 

% collect and limits x-data
x_extrema = [...
    LB, UB,...     
    LB + (-b + sqrt(b.*b - 4.*a.*c))./2./a,... % NOTE: data is offset by LB
    LB + (-b - sqrt(b.*b - 4.*a.*c))./2./a,... % NOTE: data is offset by LB
    ];

x_extrema = x_extrema(imag(x_extrema) == 0);
x_extrema = x_extrema( x_extrema >= min(x(:)) & x_extrema <= max(x(:)) );

% NOW find the peak
[y_peak, ind] = max(ppval(pp, x_extrema(:)));
x_peak = x_extrema(ind);

% y == 0%, 20% and 50%
lims = [0 20 50];
X = cell(size(lims));
Y = cell(size(lims));    
for p = 1:numel(lims)

    % the current level line to solve for
    lim = y_peak * lims(p)/100;

    % find all 3 roots of all cubics
    R = NaN(size(coefs,1), 3); 
    for ii = 1:size(coefs,1) 

        % offset coefficients to find the right intersects
        C = coefs(ii,:);
        C(end) = C(end)-lim;

        % NOTE: data is offset by LB
        Rr = roots(C) + LB(ii); 

        % prune roots
        Rr( imag(Rr)~=0 ) = NaN;
        Rr( Rr <= LB(ii) | Rr >= UB(ii) ) = NaN;
        % insert results
        R(ii,:) = Rr;
    end

    % now evaluate and save all valid points    
    X{p} = R(~isnan(R));
    Y{p} = ppval(pp, X{p});

end

% as a sanity check, plot everything 
xx = linspace(min(x(:)), max(x(:)), 20*numel(x));
yy = ppval(pp, xx);

figure(2), clf, hold on

plot(x,y, 'r') % the actual data
plot(xx,yy) % the cubic-splines interpolation 

plot(x_peak,y_peak, 'k.') % the peak

plot(X{1},Y{1}, 'r.') % the 0%  intersects
plot(X{2},Y{2}, 'g.') % the 20% intersects
plot(X{3},Y{3}, 'b.') % the 50% intersects

% finish plot
xlabel('X'), ylabel('Y'), title('Cubic splines interpolation')
legend(...
    'Real data',...
    'interpolation',...
    'peak',...
    '0% intersects',...
    '20% intersects',...
    '50% intersects',...
    'location', 'southeast')


%% (N-1)th degree polynomial
% ---------------------------

% Find best interpolating polynomial
coefs = bsxfun(@power, x, n-1:-1:0) \ y;
% (alternatively, you can use polyfit() to do this, but this is faster)

% To find the peak, we'll have to find the roots of the derivative: 
derivCoefs = (n-1:-1:1).' .* coefs(1:end-1);
Rderiv = roots(derivCoefs);
Rderiv = Rderiv(imag(Rderiv) == 0);
Rderiv = Rderiv(Rderiv >= min(x(:)) &  Rderiv <= max(x(:)));

[y_peak, ind] = max(polyval(coefs, Rderiv));
x_peak = Rderiv(ind);

% y == 0%, 20%, 50%
lims = [0 20 50];
X = cell(size(lims));
Y = cell(size(lims));
for p = 1:numel(lims)

    % the current level line to solve for
    lim = y_peak * lims(p)/100;

    % offset coefficients as to find the right intersects
    C = coefs;
    C(end) = C(end)-lim;

    % find and prune roots
    R = roots(C);
    R = R(imag(R) == 0);
    R = R(R>min(x(:)) & R<max(x(:)));

    % evaluate polynomial at these roots to get actual data
    X{p} = R;
    Y{p} = polyval(coefs, R);

end

% as a sanity check, plot everything 
xx = linspace(min(x(:)), max(x(:)), 20*numel(x));
yy = polyval(coefs, xx);

figure(3), clf, hold on

plot(x,y, 'r') % the actual data
plot(xx,yy) % the cubic-splines interpolation 

plot(x_peak,y_peak, 'k.') % the peak

plot(X{1},Y{1}, 'r.') % the 0%  intersects
plot(X{2},Y{2}, 'g.') % the 20% intersects
plot(X{3},Y{3}, 'b.') % the 50% intersects

% finish plot
xlabel('X'), ylabel('Y'), title('(N-1)th degree polynomial')
legend(...
    'Real data',...
    'interpolation',...
    'peak',...
    '0% intersects',...
    '20% intersects',...
    '50% intersects',...
    'location', 'southeast')

这导致这三个图：

线性插值多项式插值三次样条插值

（请注意，第 (N-1) 次多项式出现问题；20% 的交叉点在最后都是错误的。因此，在复制粘贴之前，请更彻底地检查所有内容：)

正如我之前所说，并且您可以清楚地看到，如果基础数据不适合它，使用单个多项式进行插值通常会引入很多问题。此外，正如您从这些图中可以清楚地看到的那样，插值方法非常强烈地影响交叉点的位置 - 至关重要的是，您至少对数据的基础模型有所了解。

对于一般情况，三次样条通常是最好的方法。但是，这是一种通用方法，会让您（以及您出版物的读者）对数据的准确性产生错误的认识。使用三次样条来初步了解相交是什么以及它们的行为方式，但一旦真实模型变得更加清晰，总是回来重新审视您的分析。当仅用于通过数据创建更平滑、更“视觉吸引力”的曲线时，当然不要使用三次样条发布 :)

score 1 · Accepted Answer

这不是一个完整的答案，但 Matlab 具有应该完成大部分您想做的事情的内置函数。

max可以帮你找到100%的线
polyfit将以最小二乘的方式为您提供一组点的多项式拟合。如果您希望它准确地通过 n 个点，我相信您至少需要使用 n-1 度。
roots会给你刚刚找到的多项式的零交叉点。您还可以使用它通过减去一个常数来找到 20% 和 50% 的交叉点。如果有多个交叉口，您将需要最接近最初找到的最大值的交叉口。（你确定过境点会一直存在吗？）

score 1 · Accepted Answer

求全局最大值使用 MAX 函数：

[ymax, imax] = max(y);
xmax = x(imax);
line(xlim,[ymax ymax],'Color','r')
line(xmax,ymax,'Color','r','LineStyle','o')

对于其余部分，您可以使用出色的 FileExchange 提交 - “Fast and Robust Curve Intersections”。

y=0 处的线可以用xlim和定义yline0 = [0 0];。然后你可以做

[x0, y0] = intersections(x,y,xlim,yline0); % function from FileExchange
x0close(1) = xmax - min(xmax-x0(x0<xmax));
x0close(2) = xmax + min(x0(x0>xmax)-xmax);
y0close = y0(ismember(x0,x0close));
line(xlim,yline0,'Color','r')
line(x0close,y0close,'Color','r','LineStyle','o')

20% 和 50% 也可以这样做，除了

yline20 = repmat((ymax - y0(1))*0.2,1,2);

所有这一切都假设您希望在绘图上出现海峡线的交叉点，而不是用于插值。

matlab - 在 Matlab 中沿绘图查找点

3 回答 3

Related

Reference