2Hz(2つの波)の​2つ目の最小値の抽出​方法について

4 ビュー (過去 30 日間)
Shogo
Shogo 2020 年 8 月 17 日
コメント済み: Takumi 2020 年 8 月 24 日
2Hz(2つの波)の2つ目の最小値の抽出方法について
大変お世話になります。
現在、サンプルデータのB列およびC列のような2つの波があるデータの2つ目の波における以下2点の行数を抽出する方法で悩んでいます。
①C列の値が0から数値が高まり、再度0になります。その後、2回目の波が現れる始まりの行数
具体的には0から10以上になった瞬間の行数
②B列の2回目の波が来たときの最小値が現れた行数
Minとfind関数で奮闘しましたが、できませんでしたので、他に方法があればご教授頂ければ幸いです。
よろしくお願いします。
Reference answers
①時間-周波数解析の実践的基礎
https://jp.mathworks.com/help/signal/examples/practical-introduction-to-time-frequency-analysis.html
②配列の最小要素
https://jp.mathworks.com/help/matlab/ref/min.html?s_tid=srchtitle
③行列から条件を指定して値を取り出すには?
https://jp.mathworks.com/matlabcentral/answers/487795-

採用された回答

Takumi
Takumi 2020 年 8 月 18 日
①および②を検出するプログラムを作成してみました.以下順を追って説明してみます.
まず添付していただいたファイルをテーブルとして読み込みます.
%% ファイル読み込み
filename = 'sample file.xlsx'; % ファイル名
T = readtable(filename); % ファイル読み込み
T.Properties.VariableNames{1} = 'Frame'; % 1列目の変数名を変更
エクセルの第一列目の変数名だけ’Frame’に変更しています.(#が変数名に使えないので)
次に,質問にある①ですが,0からの立ち上がりが2回あるので,両方を検出して二番目の立ち上がりの行数を求めることを考えました.また立ち上がりはデータ(GRF)が10以上になるところを検出したいということでしたので,まず10をしきい値として二値化してみます.
%% 2回目の立ち上がり位置のインデックス取得
mask = T.GRF>=10; % しきい値10で2値化
立ち上がり位置の検出をどのように行うか迷ったのですが,ここでは二値化した矩形波(mask)の勾配を計算し,関数findpeaksでピークを検出することで求めました.
(find関数を用いれば,idx = find(T.GRF>=10,1,'first') というようにして最初にデータが10以上になるインデックスを取得することができますが,2回目にデータが10以上になるインデックスを取得する方法がわかりませんでした.他によい方法を御存知の方がおられましたら教えて下さい)
dmask = gradient(mask); % 勾配を計算
[~,idx] = findpeaks(dmask); % 立ち上がり位置検出
idx = idx+1; % 立ち上がり位置は+1進んだところ
条件を満たす立ち上がり位置はfindpeaksによって検出されたピーク位置に+1したところです.また2回目の立ち上がり位置はidx(2)に対応します.
図で表すと以下のようになります.
figure;
subplot(3,1,1);
plot(T.Frame,T.GRF,'-r',T.Frame(idx(2)),T.GRF(idx(2)),'*b');
ylabel(T.Properties.VariableNames{3});
subplot(3,1,2);
plot(T.Frame,mask,'-r',T.Frame(idx(2)),mask(idx(2)),'*b');
ylabel('binarization')
subplot(3,1,3);
plot(T.Frame,dmask,'-r',T.Frame(idx(2)),dmask(idx(2)),'*b');
ylabel('gradient')
xlabel(T.Properties.VariableNames{1});
青いところが二回目の立ち上がり位置(10を超える)です.対応するFrameは
fprintf('2回目の波が現れる始まりのFrame(具体的には0から10以上になった瞬間のFrame)は%dです\n',T.Frame(idx(2)));
で確認してください.
次に②についてですが,こちらはfindpeaks関数で負のピーク位置を求めると良いと思います.
%% 2回目に最小になるインデックス取得
peakTH = -0.8; % ピーク検出用のしきい値
[~,minidx] = findpeaks(-T.VerticalPstn,'MinPeakHeight',peakTH); % 2つの負のピーク検出
figure;
plot(T.Frame,T.VerticalPstn,'-r');
hold on
plot(T.Frame(minidx(2)),T.VerticalPstn(minidx(2)),'*b');
xlabel(T.Properties.VariableNames{1});
ylabel(T.Properties.VariableNames{2});
具体的には,データの正負を反転させて,振幅のしきい値を指定して求めています.
2回目に最小値を取るインデックスはminidx(2)に対応します.
図で表すと以下のようになります.
青の*が2回目に最小になる位置です.
対応するFrameは
fprintf('VerticalPstnの2回目の波が来たときの最小値が現れるFrameは%dです\n',T.Frame(minidx(2)));
で確認できます.
以上参考になりましたでしょうか?一つの方法として参考にしていただけましたら幸いです.
  4 件のコメント
Takumi
Takumi 2020 年 8 月 18 日
編集済み: Takumi 2020 年 8 月 18 日
提供していただいたサンプルデータの4列目に,testデータを追加して,その列から最小値を抽出してみました.
clear
close all
clc
%% ファイル読み込み
filename = 'sample file.xlsx'; % ファイル名
T = readtable(filename); % ファイル読み込み
T.Properties.VariableNames{1} = 'Frame'; % 1列目の変数名を変更
%% 2回目の立ち上がり位置のインデックス取得
mask = T.GRF>=10; % しきい値10で2値化
dmask = gradient(mask); % 勾配を計算
[~,idx] = findpeaks(dmask); % 立ち上がり位置検出
idx = idx+1; % 立ち上がり位置は+1進んだところ
figure;
subplot(3,1,1);
plot(T.Frame,T.GRF,'-r',T.Frame(idx(2)),T.GRF(idx(2)),'*b');hold on
ylabel(T.Properties.VariableNames{3});
subplot(3,1,2);
plot(T.Frame,mask,'-r',T.Frame(idx(2)),mask(idx(2)),'*b');
ylabel('binarization')
subplot(3,1,3);
plot(T.Frame,dmask,'-r',T.Frame(idx(2)),dmask(idx(2)),'*b');
ylabel('gradient')
xlabel(T.Properties.VariableNames{1});
fprintf('2回目の波が現れる始まりのFrame(具体的には0から10以上になった瞬間のFrame)は%dです\n',T.Frame(idx(2)));
%% 2回目に最小になるインデックス取得
peakTH = -0.8; % ピーク検出用のしきい値
[~,minidx] = findpeaks(-T.VerticalPstn,'MinPeakHeight',peakTH); % 2つの負のピーク検出
figure;
plot(T.Frame,T.VerticalPstn,'-r');
hold on
plot(T.Frame(minidx(2)),T.VerticalPstn(minidx(2)),'*b');
xlabel(T.Properties.VariableNames{1});
ylabel(T.Properties.VariableNames{2});
fprintf('VerticalPstnの2回目の波が来たときの最小値が現れるFrameは%dです\n',T.Frame(minidx(2)));
%% idx(2)からminidx(2)までの間のtestの最小値
[minTest,minTestIdx] = min(T.test(idx(2):minidx(2))); % 最小値とそのインデックス
figure;
plot(T.Frame,T.test,'-r'); % testデータ全体
hold on
plot(T.Frame(idx(2):minidx(2)),T.test(idx(2):minidx(2)),'-g'); % idx(2)からminidx(2)までの間のtestデータ
plot(T.Frame(minTestIdx+idx(2)-1),minTest,'*b'); % 最小値
Shogo
Shogo 2020 年 8 月 18 日
Takumi様
早急に解決策をお教えくださり誠にありがとうございます。
こちらのサンプルデータで確かに得たいデータが算出されていることを動作確認いたしました。
今からweb検索しながら、解読して理解し、自分のプログラムの作成に移りたいと思います。
重ねて感謝申し上げます。
早朝からありがとうございました。

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

その他の回答 (2 件)

Shogo
Shogo 2020 年 8 月 19 日
Takumi様
大変お世話になります。
自分のプログラムを修正し
以下のエラーが返りますが、これは"Frame"の定義ができていないということでしょうか。
"Reference to non-existent field 'Frame'."
−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
%まずフォルダからファイルをインポートし必要な行列データをworkspaceへ
cd /Users/shogo/Desktop/Export
Files = dir('*.csv');
numfiles = length(Files);
mydata = cell(1, numfiles);
for k = 1:numfiles
myfilename = Files(k).name;
mydata{k} = readtable(myfilename,'HeaderLines',9);
T.Properties.VariableNames{9,1} = 'Frame'; %9行目,1列目の変数名を変更
GRFy = mydata{k}.Var64;
COMVPstn = mydata{k}.Var59;
KVA = mydata{k}.Var26;
%%%以下に望む処理を追加%%%
%% 2回目の立ち上がり位置のインデックス取得
mask = GRFy>=10; %しきい値10で2値化
dmask = gradient(mask); % 勾配を計算
[~,idx] = findpeaks(dmask); % 立ち上がり位置検出
idx = idx+1; % 立ち上がり位置は+1進んだところ
figure;
subplot(3,1,1);
plot(T.Frame,GRFy,'-r',T.Frame(idx(2)),GRFy(idx(2)),'*b');
ylabel(T.Properties.VariableNames{3});
subplot(3,1,2);
plot(T.Frame,mask,'-r',T.Frame(idx(2)),mask(idx(2)),'*b');
ylabel('binarization')
subplot(3,1,3);
plot(T.Frame,dmask,'-r',T.Frame(idx(2)),dmask(idx(2)),'*b');
ylabel('gradient')
xlabel(T.Properties.VariableNames{1});
fprintf('2回目の波が現れる始まりのFrame(具体的には0から10以上になった瞬間のFrame)は%dです\n',T.Frame(idx(2)));
%% 2回目に最小になるインデックス取得
peakTH = -0.8; % ピーク検出用のしきい値
[~,minidx] = findpeaks(-T.COMVPstn,'MinPeakHeight',peakTH); % 2つの負のピーク検出
figure;
plot(T.Frame,T.COMVPstn,'-r');
hold on
plot(T.Frame(minidx(2)),T.COMVPstn(minidx(2)),'*b');
xlabel(T.Properties.VariableNames{1});
ylabel(T.Properties.VariableNames{2});
fprintf('COMVPstnの2回目の波が来たときの最小値が現れるFrameは%dです\n',T.Frame(minidx(2)));
%% idx(2)からminidx(2)までの間のKVAの最小値
[minKVA,minKVAIdx] = min(T.KVA(idx(2):minidx(2))); % 最小値とそのインデックス
figure;
plot(T.Frame,T.KVA,'-r'); % KVAデータ全体
hold on
plot(T.Frame(idx(2):minidx(2)),T.KVA(idx(2):minidx(2)),'-g'); % idx(2)からminidx(2)までの間のKVAデータ
plot(T.Frame(minKVAIdx+idx(2)-1),minKVA,'*b'); % 最小値
end
−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
  1 件のコメント
Takumi
Takumi 2020 年 8 月 19 日
例ではテーブル変数Tに読み込みましたが,魚田さんはmydataに読み込んでいるので,T.*と書かれているところをmydata{k}.*に変更する必要があります.例えば
T.Properties.VariableNames{9,1} = 'Frame'; %9行目,1列目の変数名を変更
は,正しくは
mydata{k}.Properties.VariableNames{1} = 'Frame';
です.ただ,mydata{k}.Var1 mydata{k}.Var2....とやるのは面倒ですから,すでにやられているように
Frame = mydata{k}.Var1;
とするのが一番楽だと思います.
したがって最終的なプログラムはこのようにすると良いのではないでしょうか.
%まずフォルダからファイルをインポートし必要な行列データをworkspaceへ
cd /Users/shogo/Desktop/Export
Files = dir('*.csv');
numfiles = length(Files);
mydata = cell(1, numfiles);
for k = 1:numfiles
myfilename = Files(k).name;
mydata{k} = readtable(myfilename,'HeaderLines',9);
% mydata{k}.Properties.VariableNames{1} = 'Frame'; %9行目,1列目の変数名を変更
Frame = mydata{k}.Var1;
GRFy = mydata{k}.Var64;
COMVPstn = mydata{k}.Var59;
KVA = mydata{k}.Var26;
%%%以下に望む処理を追加%%%
%% 2回目の立ち上がり位置のインデックス取得
mask = GRFy>=10; %しきい値10で2値化
dmask = gradient(mask); % 勾配を計算
[~,idx] = findpeaks(dmask); % 立ち上がり位置検出
idx = idx+1; % 立ち上がり位置は+1進んだところ
figure;
subplot(3,1,1);
plot(Frame,GRFy,'-r',Frame(idx(2)),GRFy(idx(2)),'*b');
ylabel('GRFy');
subplot(3,1,2);
plot(Frame,mask,'-r',Frame(idx(2)),mask(idx(2)),'*b');
ylabel('binarization')
subplot(3,1,3);
plot(Frame,dmask,'-r',Frame(idx(2)),dmask(idx(2)),'*b');
ylabel('gradient')
xlabel('Frame');
fprintf('2回目の波が現れる始まりのFrame(具体的には0から10以上になった瞬間のFrame)は%dです\n',Frame(idx(2)));
%% 2回目に最小になるインデックス取得
peakTH = -0.8; % ピーク検出用のしきい値
[~,minidx] = findpeaks(-COMVPstn,'MinPeakHeight',peakTH); % 2つの負のピーク検出
figure;
plot(Frame,COMVPstn,'-r');
hold on
plot(Frame(minidx(2)),COMVPstn(minidx(2)),'*b');
xlabel('Frame');
ylabel('COMVPstn');
fprintf('COMVPstnの2回目の波が来たときの最小値が現れるFrameは%dです\n',Frame(minidx(2)));
%% idx(2)からminidx(2)までの間のKVAの最小値
[minKVA,minKVAIdx] = min(KVA(idx(2):minidx(2))); % 最小値とそのインデックス
figure;
plot(Frame,KVA,'-r'); % KVAデータ全体
hold on
plot(Frame(idx(2):minidx(2)),KVA(idx(2):minidx(2)),'-g'); % idx(2)からminidx(2)までの間のKVAデータ
plot(Frame(minKVAIdx+idx(2)-1),minKVA,'*b'); % 最小値
end

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


Shogo
Shogo 2020 年 8 月 19 日
Takumi 様
大変丁寧なご説明ありがとうございます。
仰るとおりです。そういうことだったのだと大変勉強になりました。
お送りさせて頂きました3つのファイル(trial 3回分)をインポートし、一人のsubject、一つのconditionでプログラムがエラーが無く動作致しました。重ねて感謝申し上げます。
残る疑問点としては以下の2点があります。
もしまだお時間許すようであればご回答頂ければ大変助かります。
どうぞよろしくお願い申し上げます。
①インポートファイル数がエラーが返ってくる件
「Index exceeds the number of array elements (1).」
こちらは,合計1529というインポートファイル数が多すぎるためでしょうか。3つのファイルであればエラーが返ってきません。
https://jp.mathworks.com/matlabcentral/answers/478489-index-exceeds-the-number-of-array-elements-1
動作確認としてグラフで確認できれば、すべてのデータ分析時には、plot部分のコードを削除すれば良いかと考えております。引き続き色々、試行錯誤していじってみます。
②抽出した2つのFrame"idx(2)", "minidx(2)",と "minKVA" の値は、「writetable」というコードで
Workspaceのtable に並べ、ファイル3つ(T1,T2,T3)で平均し、書き出し可能でしょうか。
https://jp.mathworks.com/help/matlab/import_export/exporting-to-excel-spreadsheets.html
3つのファイルでプログラムがエラーが無く完了しても、workspaceのminKVAのtableには最終の値のみしか表示されていません。
例えば、以下の表をworkspaceに書き出し、xlsxなどで別名保存することを考えております。
※インポートファイルの名前の意味は以下となっており、Trialは3回分で1つに平均したいと考えています。
「S4C1T1」→Subject 4, Condition1, Trial 1
列1「Subject Num.」(インポートファイルのSの後の値)
列2「idx(2) frame@Condition 1」
列3「minidx(2) frame@Condition 1」
列4「minKVA@Condition 1」
列5「idx(2) frame@Condition 2」
列6「minidx(2) frame@Condition 2」
列7「minKVA@Condition 2」
列8「idx(2) frame@Condition 3」
列9「minidx(2) frame@Condition 3」
列10「minKVA@Condition 3」
以上、長々と大変申し訳ありません。
よろしくお願い申し上げます。
  13 件のコメント
Shogo
Shogo 2020 年 8 月 24 日
とんでもありません。とても助かりました。
Facebookリンクなのですが、トップページに到達してしまいますので、マイプロフィールなどに繋げて頂ければと思います。
今後ともどうぞよろしくお願い致します。
Takumi
Takumi 2020 年 8 月 24 日
失礼しました。リンク編集しました

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

タグ

製品


リリース

R2020a

Community Treasure Hunt

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

Start Hunting!