关于绘图，注释中提到了一些可能的方法，而我通常会使用patch来绘制填充的多边形区域。对于面积近似，可以使用trapz函数进行梯形数值积分。

那将是我的解决方案，包括检测间隔，并忽略面积不足的间隔（这有点冗长，并且充满了绘制所有间隔的循环；可以肯定地对其进行优化）：

% Set up function,and parameter(s)
x = linspace(-0.125*pi,4.125*pi,10001);
y = linspace(60,100,10001) .* sin(x);
thr = -50;
thr_area = 30;

% Find y values lower than threshold
y_idx = find(y <= thr);

% Get start and end of intervals
idx_int = find(diff(y_idx) > 1);
n_int = numel(idx_int)+1;
s = zeros(n_int,1);
e = zeros(n_int,1);
s(1) = y_idx(1);
e(end) = y_idx(end);
for k = 1:n_int-1
  e(k) = y_idx(idx_int(k));
  s(k+1) = y_idx(idx_int(k)+1);
end

% Calculate areas
Q = zeros(n_int,1);
for k = 1:n_int
  Q(k) = abs(trapz(x(s(k):e(k)),y(s(k):e(k))-thr));
end

% Visualization
figure(1);
hold on;
plot(x,y);
xlim([x(1),x(end)]);
ylim([min(y)-10,max(y)+10]);
plot([x(1),x(end)],[thr thr],'k');
for k = 1:n_int
  patch(x(s(k):e(k)),y(s(k):e(k)),'k');
  plot([x(s(k)),x(e(k))],[y(s(k)),y(e(k))],'r.','MarkerSize',15);
  text(x(s(k)),thr+20,num2str(Q(k)));
  if (Q(k) < thr_area)
    text(x(s(k)),thr+10,'Area too low');
  else
    text(x(s(k)),'Area OK');
  end
end
hold off;

结果如下：

您现在应该已经掌握了所有信息，可以进行进一步的计算，分析等。

希望有帮助！

免责声明：我使用Octave 5.1.0测试了代码，但是我很确定，它应该与MATLAB完全兼容。如果没有，请发表评论，我将尝试解决可能的问题。

除了@HansHirse的答案外，对数据进行插值以找到阈值交叉点可能很有用。

例如，如果您的数据如下所示：

x = [ 1  2  3  4];
y = [47 49 51 53];

y不包含确切的阈值（50），因此我们可以对这些数据进行插值，以猜测根据x，我们将到达y = 50的位置。

x_interp = [ 1  2 2.5  3  4];
y_interp = [47 49 50  51 53];

没有插值的交叉点：

% Dummy data
x = 0:0.2:5*pi;
y = sin(x)*10;

% Threshold
T = 5;

% Crossing points
ind = find(abs(diff(sign(y-T)))==2)+1
xind = x(ind)
yind = y(ind)

% Plot
plot(x,y);
hold on
plot(xind,yind,'o','markersize',2,'color','r')

带有插值的交叉点：

% Dummy data
x = 0:0.2:5*pi;
y = sin(x)*10;

% Threshold
T = 5;

%% Crossing points interpolation
% Index where intersection occurs
ind = [find(abs(diff(sign(y-T)))==2)+1].'+[-1,0]

% For example we could obtain:
%        [5;              [4,5;  %We cross the threshold between y(4) and y(5)
% ind =  10;   + [-1,0] =  9,10;  %We cross the threshold between y(9) and y(10)
%        18]              17,18]  %...

xind = x(ind)
yind = y(ind)-T
% Linear interpolation
xint = xind(:,1)-yind(:,1)./(diff(yind,1,2)./diff(xind,2))
yint = T

% Plot
plot(x,y);
hold on
plot(xint,yint,'r')

然后只需将这些新的内插值添加到原始向量中即可

[x,pos] = sort([x xint]);
y       = [y yint];
y       = y(pos);

% Now apply the @HansHirse's solution.
% ...