Filter löschen
Filter löschen

面積を求めたい

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

Akzeptierte Antwort

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

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu 言語の基礎 finden Sie in Help Center und File Exchange

Produkte


Version

R2022a

Community Treasure Hunt

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

Start Hunting!