Поиск точек вдоль графика в Matlab

Question

У меня есть следующий сюжет и файл данных, который создает этот сюжет. Я хотел бы иметь Matlab найти следующие моменты для меня:

1. [у, х] для пика, отмеченного 100% линии
2. [х] для того, где участок пересекает у = 0 линии
3. [х] для где у 50% и 20% от пика находится в части 1.

существует ли дополнение инструментов или пакеты, которые люди знают о которых может помочь мне это сделать? Мне нужно сделать это для сбора сюжетов, чтобы что-то разумно автоматизированное было бы идеальным.

Я могу, конечно, программировать и вычислять детали в Matlab, это просто вопрос загрузки файла данных, сопоставления его с кривой или функцией и поиска различных [x, y] координат ,

Answer 1

ОК, здесь идет. Насколько я знаю, в 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% пересечений - все неправильно к концу. Итак, перед копированием, проверьте все немного подробнее :)

Как я уже говорил, и, как вы можете видеть, интерполирование одним полиномом часто приводит к множеству проблем, если базовые данные не подходят для этого. Кроме того, как вы можете ясно видеть из этих графиков, метод интерполяции очень сильно влияет на места пересечения - крайне важно, чтобы у вас было хотя бы какая-то идея того, какая модель лежит в основе ваших данных.

Для общих случаев кубические сплайны, как правило, являются лучшим способом. Тем не менее, это общий метод и даст вам (и читателям вашей публикации) ложное восприятие точности ваших данных. Используйте кубические сплайны, чтобы получить первое представление о том, что такое пересечения и как они себя ведут, но всегда возвращайтесь и пересматривайте свой анализ, когда реальная модель становится более ясной. Конечно, не публикуйте с кубическими сплайнами, когда это используется только для создания более плавных, более «визуально привлекательных» кривых через ваши данные :)

Answer 2

Это не полный ответ, но Matlab имеет встроенные функции, которые должны выполнить большую часть того, что вы хотите сделать.

• max может помочь вам найти 100% линии
• polyfit даст вам полиномом для множества точек, в смысле наименьших квадратов. Если вы хотите, чтобы он прошел ровно через n точек, я считаю, что вам нужно использовать хотя бы степень n-1.
• roots даст вам нулевые пересечения полинома, который вы только что нашли. Вы также можете использовать его для поиска переходов 20% и 50% путем вычитания константы. Там, где есть несколько пересечений, вам понадобятся те, которые ближе всего к вашему максимуму, найденному вами изначально. (Уверен, что переходы всегда будут существовать ли вы?)

Answer 3

Для поиска глобального максимума использования функции 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); 

---

Все это предполагает, что вы хотите пересечения пролива линий, как на вашем участке, не для интерполяции.