数据评估
Data evaluation
我在MATLAB中有一些数据,想在这些数据超过指定的阈值(例如-50
)时区分起点和终点,并保存起来,然后计算下该部分的大致面积-50
如果它低于某个定义值,则忽略这些点并检查接下来的两个点。见下图:
图中左边的两点用红色标注x
,需要的区域用绿色标注。我想对整个图做这个。
有什么想法吗?谢谢
关于绘图,评论中提到了可能的方法,而我通常会使用 patch
to plot filled polygonal regions. For the area approximation, you can use the 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)), thr+10, '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,1,2))
yint = T
% Plot
plot(x,y);
hold on
plot(xint,yint,'o','markersize',2,'color','r')
然后只需将这些新的插值添加到您的原始向量中即可:
[x,pos] = sort([x xint]);
y = [y yint];
y = y(pos);
% Now apply the @HansHirse's solution.
% ...
我在MATLAB中有一些数据,想在这些数据超过指定的阈值(例如-50
)时区分起点和终点,并保存起来,然后计算下该部分的大致面积-50
如果它低于某个定义值,则忽略这些点并检查接下来的两个点。见下图:
图中左边的两点用红色标注x
,需要的区域用绿色标注。我想对整个图做这个。
有什么想法吗?谢谢
关于绘图,评论中提到了可能的方法,而我通常会使用 patch
to plot filled polygonal regions. For the area approximation, you can use the 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)), thr+10, '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,1,2))
yint = T
% Plot
plot(x,y);
hold on
plot(xint,yint,'o','markersize',2,'color','r')
然后只需将这些新的插值添加到您的原始向量中即可:
[x,pos] = sort([x xint]);
y = [y yint];
y = y(pos);
% Now apply the @HansHirse's solution.
% ...