如何使用正交或任何其他方法查找八度音程中的点所包围的区域
How to find area enclosed by points in octave using Quadrature or any other method
我有两组坐标(正值和负值,不一定按升序排列,很多情况下不同的y值对应相同的x值 ) 我可以将其加载到两个大小相等的行向量中。
我想计算曲线围成的面积。
八度如何做到?
我试过这个 answer 但它不起作用,因为打印的区域 (204.64) 似乎太高了(见图)。
我试过代码:
function showdata(fName)
M = dlmread(fName);
H = M(2:end, 1); % starting row number is 2
B = M(2:end, 2);
aux = figure();
plot(H, B,'linewidth',2);
xlabel ("Auxilary field H (A/m)");
ylabel ("Magnetic Field B (Tesla)");
area = polyarea(H,B)
axis([min(H), max(H), min(B), max(B)]);
grid on;
grid minor on;
title (area,"fontsize",20);
然后我在 Octave 中调用 showdata('data.txt')
。
数据点图片:
This是我使用的数据文件
Octave中有一个计算凸包的函数叫做“convhull”。 return是构成凸包数据的点的索引。
M = dlmread("data.txt"); #I removed the header in data.txt
x = M(:,1);
y = M(:,2);
k = convhull(x,y);
plot (x(k), y(k), "r-", x, y, "b+");
n = rows(k);
x_prime = vertcat(x(k(n)), x(k(1:n-1)));
y_prime = vertcat(y(k(n)), y(k(1:n-1)));
A = .5*abs(x_prime'*y(k)-y_prime'*x(k)); #80.248
polyarea(x(k), y(k)) == A and true
也许凸包不是很好的面积估计,因为左上角和右下角的线离这些点有点远。还有其他方法可以从数据
形成多边形
,其中之一可能是 alpha 形状。但是alpha shape比较复杂,Octave中没有对应的预建函数。
更新:
每个 x 对应至少一个 y 坐标。我在同一个 x 上标记了最高点和最低点,并再次估计了面积。
有代码:
[uni, ~] = sort(unique(x));
n = rows(uni);
outline = [];
for i = 1:n
y_list = y(x==uni(i));
[y_max, ~] = max(y_list);
outline(i, :)= [uni(i), y_max];
[y_min, ~] = min(y_list);
outline(2*n-i+1,:)= [uni(i), y_min];
endfor
figure;
plot (x(k), y(k), "r-", x, y, "b+", outline(:,1), outline(:,2), "g-", "linewidth", 3);
polyarea(outline(:,1), outline(:,2)) #74.856
顺便说一下,如果函数 polyarea 的参数没有形成闭合曲线函数 polyarea 将 return 错误的区域。
单位正方形上的四个点:
[(0,0), (1,0), (1,1), (0,1)], [(0,0), (1,1), (1,0), (0,1) ]
polyarea([0,1,1,0],[0,0,1,1])!==polyarea([0,1,1,0],[0,1,0,1])
.
我有两组坐标(正值和负值,不一定按升序排列,很多情况下不同的y值对应相同的x值 ) 我可以将其加载到两个大小相等的行向量中。
我想计算曲线围成的面积。 八度如何做到?
我试过这个 answer 但它不起作用,因为打印的区域 (204.64) 似乎太高了(见图)。
我试过代码:
function showdata(fName)
M = dlmread(fName);
H = M(2:end, 1); % starting row number is 2
B = M(2:end, 2);
aux = figure();
plot(H, B,'linewidth',2);
xlabel ("Auxilary field H (A/m)");
ylabel ("Magnetic Field B (Tesla)");
area = polyarea(H,B)
axis([min(H), max(H), min(B), max(B)]);
grid on;
grid minor on;
title (area,"fontsize",20);
然后我在 Octave 中调用 showdata('data.txt')
。
数据点图片:
This是我使用的数据文件
Octave中有一个计算凸包的函数叫做“convhull”。 return是构成凸包数据的点的索引。
M = dlmread("data.txt"); #I removed the header in data.txt
x = M(:,1);
y = M(:,2);
k = convhull(x,y);
plot (x(k), y(k), "r-", x, y, "b+");
n = rows(k);
x_prime = vertcat(x(k(n)), x(k(1:n-1)));
y_prime = vertcat(y(k(n)), y(k(1:n-1)));
A = .5*abs(x_prime'*y(k)-y_prime'*x(k)); #80.248
polyarea(x(k), y(k)) == A and true
形成多边形
,其中之一可能是 alpha 形状。但是alpha shape比较复杂,Octave中没有对应的预建函数。
更新: 每个 x 对应至少一个 y 坐标。我在同一个 x 上标记了最高点和最低点,并再次估计了面积。 有代码:
[uni, ~] = sort(unique(x));
n = rows(uni);
outline = [];
for i = 1:n
y_list = y(x==uni(i));
[y_max, ~] = max(y_list);
outline(i, :)= [uni(i), y_max];
[y_min, ~] = min(y_list);
outline(2*n-i+1,:)= [uni(i), y_min];
endfor
figure;
plot (x(k), y(k), "r-", x, y, "b+", outline(:,1), outline(:,2), "g-", "linewidth", 3);
polyarea(outline(:,1), outline(:,2)) #74.856
polyarea([0,1,1,0],[0,0,1,1])!==polyarea([0,1,1,0],[0,1,0,1])
.