面積を求めたい

13 ビュー (過去 30 日間)
yuta
yuta 2022 年 8 月 6 日
コメント済み: yuta 2022 年 8 月 11 日
曲線をyが0以上の時の加速成分と減速成分、yが0以下の時の加速成分と減速成分の4つに分け、
それぞれ面積の合計を求めたいと考えております
% dataの準備
x = linspace(0,1,1000);
base = 4*cos(2*pi*x);
Pos = [1 2 3 5 7 8]/10;
Hgt = [3 7 5 5 4 5];
Wdt = [1 3 3 4 2 3]/100;
for n = 1:length(Pos)
Gauss(n,:) = Hgt(n)*exp(-((x - Pos(n))/Wdt(n)).^2);
end
PeakSig = sum(Gauss)+base;
% 使用するデータ
A = PeakSig(1:600);
plot(A)
% 局所最大値を求める
[pks,locs] = findpeaks(A)
pks = 1×4
6.2344 8.2539 3.7888 0.9992
locs = 1×4
101 199 299 500
% 同様に局所最小値を求める
[pks1,locs1] = findpeaks(-A)
pks1 = 1×5
-3.5528 -2.6766 -0.6439 3.4179 3.4179
locs1 = 1×5
79 142 256 423 578
% 先ほど求めた局所最大値のインデックスをもとに区分けをする
plot(A)
yline(0)
xline(locs(1)),xline(locs(2)),xline(locs(3)),xline(locs(4)),xline(locs1(1)),xline(locs1(2)),xline(locs1(3)),xline(locs1(4)),xline(locs1(5))
以下は求めたいものになります。
yが0以上である青色の面積と赤色の面積
yが0以下である緑色とオレンジ色の面積
4つの区画に分け、それぞれの合計値を算出したいと考えております。
% 正から負に変わるタイミング(yが0)がわからないため絶対値に直す
B = abs(A);
plot(B)
再度色をつけると上記のような区分けになります。
% インデックスに格納する
idx = locs1(1) >= x & x <= locs(1);
x1 = x(idx)
x1 = 1×1000
0 0.0010 0.0020 0.0030 0.0040 0.0050 0.0060 0.0070 0.0080 0.0090 0.0100 0.0110 0.0120 0.0130 0.0140 0.0150 0.0160 0.0170 0.0180 0.0190 0.0200 0.0210 0.0220 0.0230 0.0240 0.0250 0.0260 0.0270 0.0280 0.0290
% 面積を求める
trapz(0.001,x1)
ans = 0.4995
この作業を全区域で行い、それぞれ足し合わせる
現状、私ができる方法としては以上になります。
上記のように1つ1つインデックスに格納していけば面積の合計値を求めることは可能かと思いますが、
この作業を異なる1000以上の曲線で行いたいので、時間がかかりすぎてしまいます。
何か効率よく面積を求める方法がありましたらご教示いただければ幸いです。
よろしくお願いいたします。

採用された回答

Hernia Baby
Hernia Baby 2022 年 8 月 6 日
傾きはgradient関数を使えば大丈夫だと思います。
実際にやってみましょう
■前準備:同じデータです
x = linspace(0,1,1000);
base = 4*cos(2*pi*x);
Pos = [1 2 3 5 7 8]/10;
Hgt = [3 7 5 5 4 5];
Wdt = [1 3 3 4 2 3]/100;
for n = 1:length(Pos)
Gauss(n,:) = Hgt(n)*exp(-((x - Pos(n))/Wdt(n)).^2);
end
PeakSig = sum(Gauss)+base;
% 使用するデータ
A = PeakSig(1:600);
plot(A)
■条件分離
idx1 傾きが0以上
idx2 値が0以上
idx1 = gradient(A)>=0;
idx2 = A >= 0;
値が負のものB{3}とB{4}は絶対値をとっています
B{1} = A( idx1 & idx2);
B{2} = A(~idx1 & idx2);
B{3} = abs(A( idx1 & ~idx2));
B{4} = abs(A(~idx1 & ~idx2));
可視化してみます
figure
hold on
cellfun(@plot,B)
legend({'傾き:正 & 値:正','傾き:負 & 値:正','傾き:正 & 値:負','傾き:負 & 値:負'})
hold off
面積を求めます
B_sum = cellfun(@sum,B,'UniformOutput',false)
B_sum = 1×4 cell array
{[486.8008]} {[767.1826]} {[209.6346]} {[371.8635]}
cell型からdouble型に変換する場合はcell2matを使用します
B_sum = cell2mat(B_sum)
B_sum = 1×4
486.8008 767.1826 209.6346 371.8635
  5 件のコメント
Hernia Baby
Hernia Baby 2022 年 8 月 10 日
返信ありがとうございます。
そのような場合は関数を作って、最後にcellfunを適用するのがおすすめです。
面倒くさい場合も無名関数よりも関数を一度自作し、一気にやらせる方が可読性もよくなります。
load('A.mat')
SUMS = cell2mat(cellfun(@(x) MyFcn1(x),A,'UniformOutput',false))
SUMS = 4×10
1.0e+04 * 6.9582 7.1652 8.2969 8.5283 7.4717 8.8279 9.6455 6.2463 NaN 7.5476 7.6429 6.4132 8.4827 8.4332 7.1449 9.6814 9.2576 7.7036 NaN 6.9997 -2.4473 -1.1328 -3.7146 -3.8082 -1.3255 -3.3651 -3.6856 -1.8693 NaN -2.3406 -3.5207 -1.9725 -2.9688 -3.7890 -1.4947 -4.1428 -4.7756 -2.6868 NaN -2.8294
関数1
・cell内のすべてがNaNの場合、NaNで表示
・違う場合は関数2へ移行
function y = MyFcn1(x)
idx = all(isnan(x),'all');
if idx
y = repmat(x,[4 1]);
else
y = MyFcn2(x);
end
end
関数2
・条件分離し合計したものを縦に連結
function y = MyFcn2(x)
idx1 = gradient(x) >= 0;
idx2 = x >= 0;
A = sum(x( idx1 & idx2),'omitnan');
B = sum(x(~idx1 & idx2),'omitnan');
C = sum(x( idx1 &~idx2),'omitnan');
D = sum(x(~idx1 &~idx2),'omitnan');
y = vertcat(A,B,C,D);
end
yuta
yuta 2022 年 8 月 11 日
関数を作成していただきありがとうございました。
関数を作成すればこんなにも可読性がよくなるのですね。
勉強になりました。ありがとうございます。
また、いつもわかりやすくご回答いただき感謝いたします。
今後ともよろしくお願いいたします。

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File Exchange言語の基礎 についてさらに検索

製品


リリース

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!