面積を求めたい

曲線を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以上の曲線で行いたいので、時間がかかりすぎてしまいます。

何か効率よく面積を求める方法がありましたらご教示いただければ幸いです。

よろしくお願いいたします。

Answer 1

[Hernia Baby]

傾きは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

[yuta]

詳しく説明いただきありがとうございます。

こんなにもスマートに書けるのですね。助かります。

ちなみになんですが、、、

Aが既に長さ10のcell型の時(nan含む）もこのような方法が使えるのでしょうか？

load('A.mat')
A

A = 1×10 cell array
    {2582×1 double}    {1382×1 double}    {2671×1 double}    {2731×1 double}    {2425×1 double}    {3853×1 double}    {3554×1 double}    {2482×1 double}    {[NaN]}    {3913×1 double}

for ii = 1:length(A)
    if(isnan(A{ii}))
        idx1{ii} = nan;
        idx2{ii} = nan;
    else
        idx1{ii} = gradient(A{ii}) >= 0;
        idx2{ii} = A{ii} >= 0;
    end
end

idx1

idx1 = 1×10 cell array
    {2582×1 logical}    {1382×1 logical}    {2671×1 logical}    {2731×1 logical}    {2425×1 logical}    {3853×1 logical}    {3554×1 logical}    {2482×1 logical}    {[NaN]}    {3913×1 logical}

idx2

idx2 = 1×10 cell array
    {2582×1 logical}    {1382×1 logical}    {2671×1 logical}    {2731×1 logical}    {2425×1 logical}    {3853×1 logical}    {3554×1 logical}    {2482×1 logical}    {[NaN]}    {3913×1 logical}

ここまではできたのですが、この後を以前ご教示いただいた無名関数を使用してみましたが、うまくいきません。

B{1} = cellfun(@(x) x(idx1 &  idx2),A,'UniformOutput',false);
B{2} = cellfun(@(x) x(~idx1 &  idx2),A,'UniformOutput',false);
B{3} = cellfun(@(x) abs(x)(idx1 &  ~idx2),A,'UniformOutput',false);

Invalid array indexing.

Error in connector.internal.fevalMatlab

Error in connector.internal.fevalJSON

B{4} = cellfun(@(x) abs(x)(~idx1 &  ~idx2),A,'UniformOutput',false);

このようなことは可能なのでしょうか？

重ねての質問になり、申し訳ありませんがよろしくお願いいたします。

[Hernia Baby]

そうですね…少し複雑になるのでわからなければ、また聞いてください

まずはNaNしかないcellは邪魔なので省きたいと思います

clear,clc;
load('A.mat')

ここでlogicalとして抽出したいと思います

cell2mat(cellfun(@(x) all(isnan(x),'all'),A,'UniformOutput',false))

ans = 1×10 logical array
   0   0   0   0   0   0   0   0   1   0

ここでNaNの部分を消します

A(cell2mat(cellfun(@(x) all(isnan(x),'all'),A,'UniformOutput',false))) = []

A = 1×9 cell array
    {2582×1 double}    {1382×1 double}    {2671×1 double}    {2731×1 double}    {2425×1 double}    {3853×1 double}    {3554×1 double}    {2482×1 double}    {3913×1 double}

同じ操作を行います

idx1 = cellfun(@(x) gradient(x) >= 0, A, 'UniformOutput', false);
idx2 = cellfun(@(x) x >= 0, A, 'UniformOutput', false);
B = cellfun(@(x,y,z) x(y & z),A,idx1,idx2,'UniformOutput',false);
C = cellfun(@(x,y,z) x(~y & z),A,idx1,idx2,'UniformOutput',false);
D = cellfun(@(x,y,z) abs(x(y & ~z)),A,idx1,idx2,'UniformOutput',false);
E = cellfun(@(x,y,z) abs(x(~y & ~z)),A,idx1,idx2,'UniformOutput',false);

cellを縦に連結させて各々の合計値を出します

F = vertcat(B,C,D,E);
SUMs = cell2mat(cellfun(@(x) sum(x,'omitnan'), F, 'UniformOutput', false))

SUMs = 4×9

1.0e+04 * 6.9582 7.1652 8.2969 8.5283 7.4717 8.8279 9.6455 6.2463 7.5476 7.6429 6.4132 8.4827 8.4332 7.1449 9.6814 9.2576 7.7036 6.9997 2.4473 1.1328 3.7146 3.8082 1.3255 3.3651 3.6856 1.8693 2.3406 3.5207 1.9725 2.9688 3.7890 1.4947 4.1428 4.7756 2.6868 2.8294

[yuta]

返答が遅れすみません。

説明が不十分で申し訳ありません。私の取り扱っているデータの都合上、NANを省きたくなく、最後のSUMsのところはNANを含む10列で表示したいのですが、可能でしょうか？

お時間ある際に、ご教示いただければ幸いに存じます。よろしくお願いいたします。

[Hernia Baby]

返信ありがとうございます。

そのような場合は関数を作って、最後に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]

関数を作成していただきありがとうございました。

関数を作成すればこんなにも可読性がよくなるのですね。

勉強になりました。ありがとうございます。

また、いつもわかりやすくご回答いただき感謝いたします。

今後ともよろしくお願いいたします。