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-

Answer 1

[Takumi]

①および②を検出するプログラムを作成してみました．以下順を追って説明してみます．

まず添付していただいたファイルをテーブルとして読み込みます．

%% ファイル読み込み
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)));

で確認できます．

以上参考になりましたでしょうか？一つの方法として参考にしていただけましたら幸いです．

[Shogo]

Takumi様

他のスレッドでも拝見させて頂いておりました。

コメントありがとうございます。

じっくり読ませて頂きます。

取り急ぎお礼申し上げます。

[Shogo]

Takumi様

自分のプログラムにおいても上手く動作し、idx(2)とminidx(2)の行数を特定抽出することができました。

<filenameのみ編集して以下で動作確認済>>>>>>

−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

%% ファイル読み込み

filename = '/Users/shogo/Desktop/MATLAB answers用サンプル/sample file.csv'; % ファイル名

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');

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});

−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

ご指導、誠にありがとうございます。

一つ一つ、解読し自分のプログラムに当てはめて試行錯誤中です。

重ね重ね申し訳ないのですが、この2点の行間範囲内における他の列ベクトルの最小値を抽出するためのコードを作成したいと思っており調べながら徹夜したのですが方法がわかりませんでした。

以下の手順を行う方法をアドバイス頂けないでしょうか。。

%idx(2)からminidx(2)までの行数間を2nd_Phaseとし、その範囲内で別の列の最小値を抽出する

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

[Takumi]

提供していただいたサンプルデータの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]

Takumi様

早急に解決策をお教えくださり誠にありがとうございます。

こちらのサンプルデータで確かに得たいデータが算出されていることを動作確認いたしました。

今からweb検索しながら、解読して理解し、自分のプログラムの作成に移りたいと思います。

重ねて感謝申し上げます。

早朝からありがとうございました。

Answer 2

[Shogo]

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行目,１列目の変数名を変更

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)));

%% ２回目に最小になるインデックス取得

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

−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

[Takumi]

例ではテーブル変数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

Answer 3

[Shogo]

Takumi 様

大変丁寧なご説明ありがとうございます。

仰るとおりです。そういうことだったのだと大変勉強になりました。

お送りさせて頂きました3つのファイル（trial 3回分）をインポートし、一人のsubject、一つのconditionでプログラムがエラーが無く動作致しました。重ねて感謝申し上げます。

残る疑問点としては以下の２点があります。

もしまだお時間許すようであればご回答頂ければ大変助かります。

どうぞよろしくお願い申し上げます。

①インポートファイル数がエラーが返ってくる件

「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部分のコードを削除すれば良いかと考えております。引き続き色々、試行錯誤していじってみます。

②抽出した２つの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」

以上、長々と大変申し訳ありません。

よろしくお願い申し上げます。

[Takumi]

①はプログラムのどの行でエラーが出るのですか？

②はもちろん可能ですが，まず魚田さまでコードを作成してみて，わからないところがあった場合に聞いていただければと思います．

変数を保存しておきたい場合ですが，配列にすると良いと思います．例えば

[minKVA,minKVAIdx] = min(KVA(idx(2):minidx(2))); % 最小値とそのインデックス

という箇所を

[minKVA(k),minKVAIdx] = min(KVA(idx(2):minidx(2))); % 最小値とそのインデックス

とするとminKVAにすべての値を格納できると思います．

[Shogo]

早々にご回答頂き誠にありがとうございます。 先ほど質問させて頂頂き①「Index exceeds the number of array elements (1).」というエラーは、かなり多くのファイルを計算した後に、最後の行に現れます。 明日、プログラムを改めていじって課題解決に取り組みます。 ご意見に深く感謝いたします。

[Shogo]

Takumi様

昨日エラーが現れた①「Index exceeds the number of array elements (1).」というエラーはインポートファイルを少なくしたりしている間に現れなくなりました。

その代わりに別のエラー「Unable to perform assignment because the left and right sides have a different number of elements.」が現れましたが、以下の理由で解決できました。

原因：いくつかのファイルでminidx（波形2回目で最小になるインデックス取得）において、初期にマイナスの波形ができていたため、idx(2)からminidx(2)までの間と、その後定義しているのに、minidx(2)がidx(2)よりも先に現れエラーとなっていた。インポートファイルの修正で解決しました。

②3つのファイルを平均し表に書き出すコードは作成中ですのでまたアドバイス頂ければ幸いです。

取り急ぎ御報告と御礼申し上げます。

[Shogo]

Takumi様

大変お世話になっております。

教えて頂いたプログラムを自分なりに解読し、編集させて頂きました。

しかしながら、抽出データを表に書き出す部分が上手く書けずエラーが続いております。色々とAnswersなどを調べさせて頂いたのですが、どう自分のプログラムに当てはめていけばわからず苦戦しております。以下の「%% 表に結果書き出し」の部分についてアドバイス頂けないでしょうか。ご検討のほどよろしくお願い致します。

インポートフォルダ

https://www.dropbox.com/sh/menwos2rgjmzjub/AAD-oD0IrrL7L6PvSYYtISoha?dl=0

clear
close all
clc

%フォルダからファイルをインポート
cd /Users/shogo/Desktop/CMDJexportCSV
Files = dir('*.csv');
numfiles = length(Files);
mydata = cell(1, numfiles);

%一つのファイルを開けて抽出を繰り返す（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(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(k) =  (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(k),KVAminframe] = min(KVA(FC:COMmin)); % 2回目の最小値特定とそのインデックス
    KVApkTime(k) =  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');   % 最小値


 %% 表に結果書き出し　%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% 右順に列のラベルとする（ファイル名,FC,COMmin,COMminTime,KVAminframe,KVApkTime,KVAmin）
    T = table(XXXXX,FC,COMmin,COMminTime,KVAminframe,KVApkTime,KVAmin);
    T(:,7);

    %% 3つのTrialファイルごと(上から3行ごと)に平均する(移動平均で可能？)
    M = movmean(T,3) % 3試技ごとに移動平均算出
  
    %% C1,C2,C3で別々のExcelシートへ書き出し
    filename = 'MATLAB抽出後データ.xlsx';
    writetable(T,Files(k).name,'Sheet','MATLAB抽出後データ','WriteVariableNames',false);


end
fprintf('<<<<DONE!>>')

[Takumi]

何をどのように出力したいのか，どのようなエラーが出るのか具体的に教えていただけますか？

最終的に出力したい結果を手打ちでもいいのでExcelファイルで再現していただけると指摘しやすいです．

あと，平均のところも何の平均かわからないので教えて下さい．

[Shogo]

大変お世話になります。

添付ファイルのように、まず一枚目のシートに抽出した値を書き出し、フォルダ内のインポートファイル名がT1,2,3と続きますので、こちらを平均し、シート2,3,4,5枚目にインポートファイルのC1,2,3,4ごとに整理、書き出しを行いたいと考えています。

まず躓いているのは、Excelの一行目にファイル名を行ベクトルで書き出すために、XXXXXの部分へ「Files」や「myfilename」と代入してみましたが上手く動きませんでした。

%% 表に結果書き出し　%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% 右順に列のラベルとする（ファイル名,FC,COMmin,COMminTime,KVAminframe,KVApkTime,KVAmin）
    T = table(XXXXX,FC,COMmin,COMminTime,KVAminframe,KVApkTime,KVAmin);
    T(:,7);
    %% 3つのTrialファイルごと(上から3行ごと)に平均する(移動平均で可能？)
    M = movmean(T,3) % 3試技ごとに移動平均算出
  
    %% C1,C2,C3で別々のExcelシートへ書き出し
    filename = 'MATLAB抽出後データ.xlsx';
    writetable(T,Files(k).name,'Sheet','MATLAB抽出後データ','WriteVariableNames',false);

その後、Excelの一枚目シートをもとに、T1,2,3を移動平均などで、値を平均化し、シート2,3,4,5枚目に、それぞれC1,2,3,4ごとに整理できればと考えています。

インポートフォルダ

https://www.dropbox.com/sh/menwos2rgjmzjub/AAD-oD0IrrL7L6PvSYYtISoha?dl=0

なんとか仕組みを理解したいと思うのですが、知識不足でアイディアがなく前進できていません。色々と無理を言い大変申し訳ありませんがアドバイス頂ければ幸いです。

何卒よろしくお願い致します。

[Takumi]

丁寧にありがとうございます．

まず，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]

エクセルの出力方法についてサンプルまで作成頂きご説明ありがとうございます。

自分なりに編集して、3つ目の方法で出力を試みました。

以下のエラーが消えないのですが、{}この中括弧が問題なのでしょうか？

Brace indexing is not supported

for variables of this type.

https://jp.mathworks.com/matlabcentral/answers/401438-error-brace-indexing-is-not-supported-for-variables-of-this-type

なにかご存知であればご教授頂ければ幸いです。

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

clear
close all
clc

%フォルダからファイルをインポート
cd /Users/shogo/Desktop/CMDJexportCSV
Files = dir('*.csv');
numfiles = length(Files);
mydata = cell(1, numfiles);

%一つのファイルを開けて抽出を繰り返す（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(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(k) =  (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(k),KVAminframe] = min(KVA(FC:COMmin)); % 2回目の最小値特定とそのインデックス
    KVApkTime(k) =  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');   % 最小値


end



 %% 表に結果書き出し　%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>> % セル配列で書き込む場合
header = {'Files','FC','COMmin','COMminTime','KVAminframe','KVApkTime','KVAmin'};

>> % ヘッダだけ先に書き込んで順番に変数を書き込んで行きたい場合
writecell(header,'myData.xlsx'); % ヘッダを書き込む
nfile = length(Files);
for i=1:nfile
    fname = Files{i};
    a = FC(i);
    b = COMmin(i);
    c = COMminTime(i);
    d = KVAminframe(i);
    e = KVApkTime(i);
    f = KVAmin(i);
    str = [fname,num2cell([a,b,c,d,e,f])];
    writecell(str,'myData.xls','Range',['FC',num2str(i+1,'%d')]);    % ２行目以降にデータを書き込む
end


fprintf('<<<<DONE!>>')

[Takumi]

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]

Takumi様

プログラムコードを修正してくださり、誠にありがとうございます。

こんなにシンプルにできるものなのですね。本当に凄いです。

残る作業もなんとか取り組んでいきます。

ここまで大変親切にして頂き、ぜひ御礼させて頂きたいのでFacebookなどで繋がることは可能でしょうか。

[Takumi]

私もまだまだ未熟でスマートな書き方ではないかもしれませんが、少しでも参考にしていただけたなら回答者としても大変嬉しく思います。 御礼なんて気を使っていただく必要はないのですが、私のプロフィールにFacebookのリンクを設けておきましたので、またなにかお困りでしたらそちらから連絡いただいても結構です。 よろしくお願いいたします。

[Shogo]

とんでもありません。とても助かりました。

Facebookリンクなのですが、トップページに到達してしまいますので、マイプロフィールなどに繋げて頂ければと思います。

今後ともどうぞよろしくお願い致します。

[Takumi]

失礼しました。リンク編集しました