2Hz（2つの波）の2つ目の最小値の抽出方法について
显示 更早的评论
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
2020-8-18
①および②を検出するプログラムを作成してみました．以下順を追って説明してみます．
まず添付していただいたファイルをテーブルとして読み込みます．
%% ファイル読み込み
filename = 'sample file.xlsx'; % ファイル名
T = readtable(filename); % ファイル読み込み
T.Properties.VariableNames{1} = 'Frame'; % １列目の変数名を変更
エクセルの第一列目の変数名だけ’Frame’に変更しています．（＃が変数名に使えないので）
次に，質問にある①ですが，０からの立ち上がりが２回あるので，両方を検出して二番目の立ち上がりの行数を求めることを考えました．また立ち上がりはデータ(GRF)が10以上になるところを検出したいということでしたので，まず10をしきい値として二値化してみます．
%% 2回目の立ち上がり位置のインデックス取得
mask = T.GRF>=10; % しきい値10で2値化
立ち上がり位置の検出をどのように行うか迷ったのですが，ここでは二値化した矩形波(mask)の勾配を計算し，関数findpeaksでピークを検出することで求めました．
（find関数を用いれば，idx = find(T.GRF>=10,1,'first') というようにして最初にデータが10以上になるインデックスを取得することができますが，２回目にデータが10以上になるインデックスを取得する方法がわかりませんでした．他によい方法を御存知の方がおられましたら教えて下さい）
dmask = gradient(mask); % 勾配を計算
[~,idx] = findpeaks(dmask); % 立ち上がり位置検出
idx = idx+1; % 立ち上がり位置は+1進んだところ
条件を満たす立ち上がり位置はfindpeaksによって検出されたピーク位置に＋１したところです．また２回目の立ち上がり位置は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関数で負のピーク位置を求めると良いと思います．
%% ２回目に最小になるインデックス取得
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});
具体的には，データの正負を反転させて，振幅のしきい値を指定して求めています．
２回目に最小値を取るインデックスはminidx(2)に対応します．
図で表すと以下のようになります．
青の＊が２回目に最小になる位置です．
対応するFrameは
fprintf('VerticalPstnの2回目の波が来たときの最小値が現れるFrameは%dです\n',T.Frame(minidx(2)));
で確認できます．
以上参考になりましたでしょうか？一つの方法として参考にしていただけましたら幸いです．
4 个评论
Shogo
2020-8-18
Shogo
2020-8-18
提供していただいたサンプルデータの4列目に，testデータを追加して，その列から最小値を抽出してみました．
clear
close all
clc
%% ファイル読み込み
filename = 'sample file.xlsx'; % ファイル名
T = readtable(filename); % ファイル読み込み
T.Properties.VariableNames{1} = 'Frame'; % １列目の変数名を変更
%% 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)));
%% ２回目に最小になるインデックス取得
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
2020-8-18
更多回答（2 个）
Shogo
2020-8-19
0 个投票
1 个评论
Takumi
2020-8-19
例ではテーブル変数Tに読み込みましたが，魚田さんはmydataに読み込んでいるので，T.*と書かれているところをmydata{k}.*に変更する必要があります．例えば
T.Properties.VariableNames{9,1} = 'Frame'; %9行目,１列目の変数名を変更
は，正しくは
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行目,１列目の変数名を変更
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)));
%% ２回目に最小になるインデックス取得
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
2020-8-19
0 个投票
13 个评论
Takumi
2020-8-19
①はプログラムのどの行でエラーが出るのですか？
②はもちろん可能ですが，まず魚田さまでコードを作成してみて，わからないところがあった場合に聞いていただければと思います．
変数を保存しておきたい場合ですが，配列にすると良いと思います．例えば
[minKVA,minKVAIdx] = min(KVA(idx(2):minidx(2))); % 最小値とそのインデックス
という箇所を
[minKVA(k),minKVAIdx] = min(KVA(idx(2):minidx(2))); % 最小値とそのインデックス
とするとminKVAにすべての値を格納できると思います．
Shogo
2020-8-19
Shogo
2020-8-20
Takumi
2020-8-23
丁寧にありがとうございます．
まず，XXXXXの部分を例えばmyfilenameに変更して実行してみるとtableを作成するときに次のようなエラーが出ますね．
エラー:
T = table(myfilename,FC,COMmin,COMminTime,KVAminframe,KVApkTime,KVAmin);
原因:
文字ベクトル 'S4C1T1.csv' を変数の 1 つとしてもつ 1 行の table を作成しようとした可能性があります。テキスト データを table に保存するには、文字配列ではなく、string 配列または文字ベクトルの cell 配列を使用してください。または、1 行の cell 配列を作成してからCELL2TABLE を使用して table に変換してください。
これは，myfilenameが一つの要素しかなくてもcell配列にしてくださいねということなので，次のように1×1のcell配列にするとエラーを回避できます．また，table変数を作成するときはすべての変数の行数が揃っている必要があることに注意してください．いまはすべての変数が1行ですね．
fname = {myfilename};
T = table(fname,FC,COMmin,COMminTime,KVAminframe,KVApkTime,KVAmin);
iが1のときのTの出力は以下のとおりです．
T =
1×7 table
fname FC COMmin COMminTime KVAminframe KVApkTime KVAmin
______________ ____ ______ __________ ___________ _________ _______
{'S4C1T1.csv'} 6672 7186 0.36715 205 0.14643 -16.284
これをエクセルに書き込みたい場合は
writetable(T,'MATLAB抽出後データ.xlsx');
とすればOKです．ただ，for文の中にこのまま書くと毎回上書きされてしまうことが容易にわかりますね．
したがって，エクセル出力は最後に一気にやるといいと思うのですが，FC，COMmin等の変数を配列にしたり大変だと思いますし，イメージが沸かないと思いますので，練習としてエクセル出力のかんたんな例を作成しました．
まず，エクセル出力にはいくつか方法があります．一つはwritetable関数を使う方法ですが，もう一つおすすめの方法はwritecell関数を使う方法です．（あとwritematrixもあります，他にもいくつかあります．．．）
２つの関数を使って出力するサンプルプログラムを作りました．
まずテストデータを作ります．
% テストデータ
filename = {'test1';'test2';'test3'};
A = [1;2;3];
B = [4;5;6];
C = [7;8;9];
それぞれ3×1の配列です．テーブル変数を作成して一気に書き込む場合はwritetable関数を使ってこのようにかけます．
% テーブルで書き込む場合
T = table(filename,A,B,C);
writetable(T,'myDataTable.xls')
次に同じことをwritecell関数でやってみます．writetableではテーブル変数名が自動的にエクセルの1行目に追加されますが，writecellでは自分で書きます．まずヘッダの文字列を作成します．
% セル配列で書き込む場合
header = {'filename','A','B','C'};
次に書き込むデータをまとめましょう．
data = [filename,num2cell([A,B,C])]; % cell配列に変換
writecell関数ではすべてcell配列にする必要があるので，数値データはcell配列に変換してからfilename配列と結合していることに注意してください．
最後にヘッダーとデータを結合して一気に書き込みます．
writecell([header;data],'myDataCell.xls')
myDataTable.xlsとmyDataCell.xlsが全く同じ出力になることを確認してください．
最後にデータを順々に書き込んでいきたい場合の例を示します．
まずエクセルの1行目に変数名だけ先に書き込みます．
% ヘッダだけ先に書き込んで順番に変数を書き込んで行きたい場合
writecell(header,'myData.xls'); % ヘッダを書き込む
次に2行目以降に追加で書き込んでいきましょう．
nfile = length(filename);
for i=1:nfile
fname = filename{i};
a = A(i);
b = B(i);
c = C(i);
str = [fname,num2cell([a,b,c])];
writecell(str,'myData.xls','Range',['A',num2str(i+1,'%d')]); % ２行目以降にデータを書き込む
end
ここではせっかく配列になっているfilename,A,B,Cをあえてバラバラにして書き込んでいっています．
そして書き込む行数をRangeで指定しています．
以上3通りの方法を説明しましたが，最後の方法がやりやすいんじゃないかと思います．
これでエクセル出力については解決しましたでしょうか．．．？
Shogo
2020-8-23
Takumi
2020-8-24
fname = Files{i};
ではなく
fname = Files(i).name;
です．
あと，いくつか間違っているので編集しました．
clear
close all
clc
%フォルダからファイルをインポート
cd /Users/shogo/Desktop/CMDJexportCSV
Files = dir('*.csv');
numfiles = length(Files);
mydata = cell(1, numfiles);
% ヘッダーだけ先に書き込む
header = {'Files','FC','COMmin','COMminTime','KVAminframe','KVApkTime','KVAmin'};
writecell(header,'myData.xlsx'); % ヘッダを書き込む
%一つのファイルを開けて抽出を繰り返す（for~end）
for k = 1:numfiles
myfilename = Files(k).name
mydata{k} = readtable(myfilename,'HeaderLines',9);
Frame = mydata{k}.Var1;
GRFy = mydata{k}.Var64;
COMVPstn = mydata{k}.Var59;
KVA = mydata{k}.Var26;
%%%以下に望む処理を追加%%%
%% 2nd_Landing_Foot_Contact特定（波形2回目のGRFy立ち上がり位置のindexA取得）
mask = GRFy>=10; %しきい値10で2値化
dmask = gradient(double(mask)); % 勾配を計算
[~,idxA] = findpeaks(dmask); % 立ち上がり位置検出
idxA = idxA+1; % 立ち上がり位置は+1進んだところ
figure;
subplot(3,1,1);
plot(Frame,GRFy,'-r',Frame(idxA(2)),GRFy(idxA(2)),'*b');
ylabel('GRFy');
subplot(3,1,2);
plot(Frame,mask,'-r',Frame(idxA(2)),mask(idxA(2)),'*b');
ylabel('binarization')
subplot(3,1,3);
plot(Frame,dmask,'-r',Frame(idxA(2)),dmask(idxA(2)),'*b');
ylabel('gradient')
xlabel('Frame');
fprintf('2回目FootContactFrameは%dです\n',Frame(idxA(2)));
%% 2nd_LandingCOMmin特定（波形2回目でCOMVPstnが最小になるindexB取得）
peakTH = -1; % ピーク検出用のしきい値（COMが1m未満で検索できるように）
[~,idxB] = findpeaks(-COMVPstn,'MinPeakHeight',peakTH); % 2つの負のピーク検出
COMminTime = (idxB(2) - idxA(2))*0.00071429; % 接地からCOMminまでの秒数(1400Hz)
figure;
plot(Frame,COMVPstn,'-r');
hold on
plot(Frame(idxB(2)),COMVPstn(idxB(2)),'*b');
xlabel('Frame');
ylabel('COMVPstn');
fprintf('2回目COMminFrameは%dです\n',Frame(idxB(2)));
%% FCからCOMminまでの間（2nd_Landing）のKVAの最小値
FC= idxA(2);
COMmin= idxB(2);
[KVAmin,KVAminframe] = min(KVA(FC:COMmin)); % 2回目の最小値特定とそのインデックス
KVApkTime = KVAminframe*0.00071429; % 接地からKVApkまでの秒数(1400Hz)
figure;
plot(Frame,KVA,'-r'); % KVAデータ全体
hold on
plot(Frame(FC:COMmin),KVA(FC:COMmin),'-g'); % FCからCOMminの間のKVA最小値
plot(Frame(KVAminframe+FC-1),KVAmin,'*b'); % 最小値
%% 表に結果書き出し %%
str = [myfilename,num2cell([FC,COMmin,COMminTime,KVAminframe,KVApkTime,KVAmin])];
writecell(str,'myData.xlsx','Range',['A',num2str(k+1,'%d')]); % ２行目以降にデータを書き込む
end
fprintf('<<<<DONE!>>')
これでいかかでしょうか？
Shogo
2020-8-24
Takumi
2020-8-24
私もまだまだ未熟でスマートな書き方ではないかもしれませんが、少しでも参考にしていただけたなら回答者としても大変嬉しく思います。 御礼なんて気を使っていただく必要はないのですが、私のプロフィールにFacebookのリンクを設けておきましたので、またなにかお困りでしたらそちらから連絡いただいても結構です。 よろしくお願いいたします。
Shogo
2020-8-24
Takumi
2020-8-24
失礼しました。リンク編集しました
类别
在 帮助中心 和 File Exchange 中查找有关 フーリエ解析とフィルター処理 的更多信息
