2Hz(2つの波)の2つ目の最小値の抽出方法について
8 views (last 30 days)
Show older comments
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-
0 Comments
Accepted Answer
Takumi
on 18 Aug 2020
①および②を検出するプログラムを作成してみました.以下順を追って説明してみます.
まず添付していただいたファイルをテーブルとして読み込みます.
%% ファイル読み込み
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 Comments
Takumi
on 18 Aug 2020
Edited: Takumi
on 18 Aug 2020
提供していただいたサンプルデータの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'); % 最小値
More Answers (2)
Shogo
on 19 Aug 2020
1 Comment
Takumi
on 19 Aug 2020
例ではテーブル変数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
See Also
Categories
Find more on 記述統計 in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!