面積を求めたい

7 views (last 30 days)
yuta
yuta on 6 Aug 2022
Commented: yuta on 11 Aug 2022
曲線を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以上の曲線で行いたいので、時間がかかりすぎてしまいます。
何か効率よく面積を求める方法がありましたらご教示いただければ幸いです。
よろしくお願いいたします。

Accepted Answer

Hernia Baby
Hernia Baby on 6 Aug 2022
傾きは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 Comments
yuta
yuta on 11 Aug 2022
関数を作成していただきありがとうございました。
関数を作成すればこんなにも可読性がよくなるのですね。
勉強になりました。ありがとうございます。
また、いつもわかりやすくご回答いただき感謝いたします。
今後ともよろしくお願いいたします。

Sign in to comment.

More Answers (0)

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!