Вычислить площадь выше и ниже установленного уровня, используя trapz MATLAB в диапазоне

Question

Data Plot

 %% Data 
% Imports the array data (x and y) to the workspace by loading Excel data 
% that has been imported and converted to *.mat format. 
    load('900day_r') % y data. 
    load('x_degreesb'); % x axis data 

%% Remove non linear trends 
    opol = 0;% Degree of filtering on original profile 
    [p,s,mu] = polyfit(x_degreesb,x900day_r,opol); 
    f_y = polyval(p,x_degreesb,[],mu); 
    x900day_r = x900day_r - f_y; 

    max_x = max(x900day_r);% Find maximum in array 
    x900day_r = x900day_r/max_x;% Normalize height to max 
    min_x = min(x900day_r);% Find minimum in array 
    x900day_r = x900day_r - min_x;% Shift profile (lowest value in array = 0) 

%% Find Peaks & Valleys 
    [pks, locs] = findpeaks(x900day_r); % returns peaks & locations 
    x900day_r_Inv = max(x900day_r)-x900day_r; % invert y data 
    vlys = max(vlys)-vlys; % invert data for valley markers 

%% Plot Profile 
% Plot profile and markers for peaks 
    plot(x_degreesb,x900day_r,'b',x_degreesb(locs),pks+0.04,'v','markersize',5,'markerfacecolor','b','linewidth',1); 
    hold on 
% Plot profile and markers for valleys 
    plot(x_degreesb(min_locs),vlys-0.04,'^','markersize',5,'markerfacecolor','b','linewidth',1); 

% Plot characteristics 
    axis('tight') % Makes the graph fill the figure. 
    grid on % Turns the grid on (major not minor). 

% Sets up the figure (fig) for display 
    fig = gca; 
    set(fig,'fontname','Bodoni 72','fontsize',20); 

% Set y limits to auto. Set x limits and x tick marks 
    ylim('auto'); % Sets the y limits 
    xlim([0, 360]); % Sets the x limits from 0 to 360 
    set (fig, 'XTick', [0, 45, 90, 135, 180, 225, 270, 315, 360]) % Sets x axis tick marks 

% Set fig for presentation appearance 
    x = xlabel('$Location On Cylinder Perimiter {[degrees]}$'); % x label. 
    y = ylabel('$Curve Height {[mm]}$'); % y label. 
    t = title('$900 Days Cylinder$ r'); % Title (presentation fraction). 

% Set vertical lines at quadrant boundaries 
    hold on 
    x1 = 90; 
    x2 = 180; 
    x3 = 270; 
    x4 = 360; 
    y1 = get(gca, 'ylim'); 
    plot ([x1 x1],y1,'k')%Caudal Medial(0:90) 
    plot ([x2 x2],y1,'k')%Cranial Medial(90:180) 
    plot ([x3 x3],y1,'k')%Cranial Lateral(180:270) 
    plot ([x4 x4],y1,'k')%Cadual Lateral(270:360) 
    hold on 

% Interpretation of text characters for presentation 
    set(t,'Interpreter','Latex'); 
    set(x,'Interpreter','Latex'); 
    set(y,'Interpreter','Latex'); 

%% Isolate Cranial Medial & Lateral Section of Profile 
    x = x_degreesb;% Quadrant data 
    y = x900day_r;% Profile data 

    indx = (x >= 90) & (x <= 270);% Index section 
    [pks, locs, widths, proms] = findpeaks(y(indx));% Find peaks in section 

    Rmax = max(pks);% Find maximum peak  
    Rmin = min(y(indx));% Find minimum in section 
    CL = (Rmax + Rmin)/2;% Center line of sectioned profile  

%% Plot Center Line 
    hold on 
    x1 = 90; 
    x2 = 270; 
    y1 = CL; 
    plot ([x1 x2],[y1 y1],'r','linewidth',1.5); 

%% Plot Rmax 
    hold on 
    x1 = 90; 
    x2 = 270; 
    y1 = Rmax; 
    plot ([x1 x2],[y1 y1],'k','linewidth',1.5); 

%% Plot Rmin 
    hold on 
    x1 = 90; 
    x2 = 270; 
    y1 = Rmin; 
    plot ([x1 x2],[y1 y1],'k','linewidth', 1.5); 

%% Highlight Region of Interest 
%subplot(2,1,2) 
    x = x_degreesb;% Quadrant data 
    y = x900day_r;% Profile data 
    indx = (x >= 90) & (x <= 270);% Index section 

    plot(x(indx),y(indx),'k','linewidth',2)% Plot 90:270 curve in black 
    level = CL;% set the centerline as the level for area shading 
% Shade the area of the curve in the indexed section grey [.7 .7 .7] above the level CL 
    area(x(indx), max(y(indx), level),level, 'EdgeColor', 'none', 'FaceColor',[.7 .7 .7],'showbaseline', 'off'); 
% Shade the area of the curve in the indexed section dark grey below the level CL 
    area(x(indx), min(y(indx), level),level, 'EdgeColor', 'none', 'FaceColor',[.5 .5 .5],'showbaseline','off'); 

Кто-нибудь знает, как я могу найти площадь выше (светло-серый) и ниже (темно-серый) средняя линия (красная) для конкретного диапазона между 90: 270 с использованием MATLAB? Я пытаюсь использовать trapz, устанавливая уровень (красная линия в изображении Data Plot), но, похоже, не может получить trapz для вычисления только выделенных областей. Я разместил код, но не данные, как его довольно большой набор данных, составляющий кривую. Любая помощь будет принята с благодарностью!

RP

@Some Guy: Спасибо за быстрый ответ. Вот пример скрипта, который вы можете запустить. Вы можете изменить соотношение синей области к серой области, изменив значение уровня. В pic 1 с уровнем, установленным в 2, они должны быть равны. В pic 2 синий цвет должен быть более серого с уровнем, установленным в 1.6. То, что я пытался сделать. Есть предположения?

%% Example 
x = 0:.01:4*pi;% x data 
y = sin(x)+2;% y data 
level = 1.6;% level 
plot(x, y) 
hold on 
x_interest = 0:.01:x(length(y)); 
y_interest = sin(x_interest)+2; 
xlim ([0 x(length(y))]) 

% Shaded area above level 
area(x_interest, max(y_interest, level), level, ... 
    'EdgeColor', 'none', 'FaceColor', [.6 .7 .8], ... 
    'ShowBaseLine', 'off'); 

% Shaded area below level 
area(x_interest, min(y_interest, level), level, ... 
    'EdgeColor', 'none', 'FaceColor', [.5 .5 .5], ... 
    'ShowBaseLine', 'off'); 

y_above = max(y_interest - level,0); % Take only the part of curve above y_split 
y_below = max(-y_above,0); % Take the part below y_split 
A_above = trapz(y_above) 
A_below = trapz(y_below) 

Answer 1

Если бы я имел данные в векторе y и скалярную y_split, на котором я хотел расколоть я хотел бы сделать:

y_above = y - y_split; 
y_below = max(-y_above,0); % Take the part below y_split 
y_above = max(y_above,0); % Take the part above y_split 
A_above = trapz(y_above); 
A_below = trapz(y_below); 

Вы можете построить y_above и y_below к убедитесь, что вы интегрируетесь по своему усмотрению.

EDIT: С OPs пример сценария области с level = 2 является:

A_above = 

    399.9997 


A_below = 

    399.9976 

И level = 1.6 является

A_above = 

    683.5241 


A_below = 

    181.1221